Genetic variation of the Slc11a2gene was analyzed from data obtained by sequencing of exonic regions of this innate immune gene which was studied to screen SNPs. Ten Bos taurus animals of Gir breed and ten Bubalus bubalis animals of Murrah breed were used for genomic DNA extraction (unpublished data). GATK software tools (version 2.8; http://www.broadinstitute.org) were used for genotype calling with recommended parameters. Genotypes were called by the GATK Unifiedgenotyper tool, and variants were filtered by depth 60.

SIFT performs multiple alignments of a number of peptide sequences until a median conservation for the sequence is reached at the default of 3.0 and then it predicts whether substitution with any of the other amino acids is tolerated or deleterious for every position in the submitted sequence [20]. The SIFT prediction was given as a tolerance index (TI) score ranging from 0.0 to 1.0, which was the normalized probability that the amino acid change was tolerated. A nsSNP with a TI score of V0.05 was considered to be deleterious i.e. an amino acids with probabilities < 0.05 were predicted to be deleterious. We submitted the amino acid sequence of Slc11a2 along with nsSNPs with corresponding amino acid positions.

The functional validation of nsSNPs predicted by SIFT was analyzed by PANTHER (Protein Analysis Through Evolutionary Relationships;www.pantherdb.org/tools/csnp). This tool estimates the likelihood of a particular non-synonymous coding SNP to cause a functional impact on the protein using Hidden Markov Models (HMM) based modeling and evolutionary relationship. It calculates the subPSEC (substitution position-specific evolutionary conservation) score based on an alignment of evolutionarily related proteins [5]. The score of subPESC ≥− 3 was predicted as a less deleterious, while ≤− 3 was predicted as the deleterious effect. Amino acid sequence in FASTA format was uploaded.

PROVEAN (Protein Variation Effect Analyzer) is a tool which predicts the impact of an amino acid substitution or indel on the biological function of a protein (http://provean.jcvi.org/index.php).This algorithm allows for the best balanced separation between the deleterious and neutral amino acids, based on a threshold. The score <− 2.5 indicates that the variant is deleterious and >− 2.5 score is considered as a neutral variant [21]. A query peptide sequence of Slc11a2 was provided in FASTA format to the PROVEAN server for predicting the functional impact of the SNPs.

I-Mutant2.0 (http://folding.biofold.org/cgi-bin/i-mutant2.0) is a Support Vector Machine-based web server for the automatic prediction of protein stability changes upon single-site mutations. The input FASTA sequence of protein along with the residues change was provided for analysis of DDG value (kcal/mol) [22]. Also the RI value (reliability index) was computed.

RaptorX, a protein structure prediction server,predicts 3D structures for protein sequences lacking close homologs in the Protein Data Bank (PDB).For given FASTA sequence RaptorX predicted its secondary and tertiary structures as well as solvent accessibility and disordered regions. RaptorX also calculates p-value for the relative global quality, GDT (global distance test) and uGDT (un-normalized GDT) for the absolute global quality, and RMSD for the absolute local quality of each residue in the model. The 3D structures were visualized by PyMOL (http://www.pymol.org/) which is an open source molecular isualization tool. Mutant model was also constructed using PyMoL tool.

Model quality was checked both of native and altered protein by Ramachandran plot using software RAMPAGE (mordred.bioc.cam.ac.uk/~rapper/rampage.php) which analyzed residue-by-residue geometry and overall structure geometry. PyMOL was used to locate nsSNPs on protein structure and for analyzing RMS deviation by superimposing both native and mutant structures. Amino acids at the position of SNPs were checked for polar interactions with other amino acids in the protein using PyMOL. In addition, total energy after energy minimization was calculated for each altered model using Swiss PDB viewer.

To find whether these identified nsSNPs are present on any epitope region or any protein binding region, we performed binding site prediction using RaptorX Binding and FT site server which predicted binding site regions in Slc11a2 protein.

Upon variant calling, total 124 SNPs were observed in Slc11a2 gene (Supplementary Table 1). Among these SNPs, 22 (17.74%) and 74 (59.67%) were found to be non-synonymous and synonymous respectively. The remaining 28 (22.58%) were found to be in the non-coding region, 3 in UTR 5′region and 25 in UTR 3′ region.

The SIFT identified 8 nsSNPs viz. I114T, G369R, Y374C, A377V, Q385H, M389V, N492S and V497M to be deleterious. Non-synonymous SNPs with SIFT prediction and SIFT score are shown in Table 1.

The results of SIFT were further confirmed by investigating the effect of nsSNPs on protein function using HMM based PANTHER tool. The analysis of 22 non-synonymous mutations revealed that 5 SNPs (G369R, Y374C, A377V, Q385H and N492S) reflected a subPSEC score >−3, thus PANTHER classified them as deleterious. Remaining SNPs of Slc11a2 had a score <− 3 and were classified as tolerated. Non-synonymous SNPs along with PANTHER score are given in Table 1.

Further confirmation of effect of nsSNPs on protein was done using PROVEAN tool which revealed 8 from 22 nsSNPs (I114T, G369R, Y374C, A377V, Q385H, M389V, S490F and N492S) to be deleterious. The higher the tolerance index is, the less functional impact a particular amino acid substitution is likely to have, and vice versa. Among the 22 nsSNPs, 8 (36.36%) were found to be deleterious, having a tolerance index score of ≤− 2.5 using PROVEAN tool (Table 1).

To add another layer of confirmation, we also analyzed effect of these nsSNPs using I-Mutant 2.0. It gave result in the form of effect of mutants on stability of protein with reliability index at pH 7.0 and temperature 25 °C. Here in our case, in Slc11a2, all 22 non-synonymous SNPs showed resulting decrease in stability of the protein. All 22 SNPs with reliability index and DDG value are given in Table2.

The 3D structure of native model generated through RaptorX was visualized using PyMoL. Slc11a2 is having 568 amino acid residues (Supplementary Fig. 1). From these 480(85%) residues were modeled and 65(11%) positions predicted as disordered. Secondary structures revealed 69% helix, 0% beta sheet and 30% loop structures. The solvent accessibility was divided into three states by 2 cut-off values: 10% and 42%. Value less than 10% was identified as buried, larger than 42% value was identified as exposed and if value was between 10% and 42% was identified as medium. Proportions of buried, medium and exposed regions in our protein were 62%, 21% and 15% respectively. Overall uGDT (GDT) value was 143(25%). The uGDT is the unnormalized GDT (global distance test) score. For a protein with > 100 residues, uGDT > 50 is a good indicator. For a protein with < 100 residues, GDT > 50 is a good indicator. If a model has acceptable uGDT (> 50) but lower GDT (< 50), it indicates that only a small portion of the model may be good. P-value is the likelihood of a predicted model being worse than the best of a set of randomly-generated models for this protein (or domain), so P-value evaluates the relative quality of a model. The smaller the p-value, the higher is the quality of the model. For alpha proteins, p-value less than 10^− 3 is a good indicator. For manly beta proteins, p-value less than 10^− 4 is a good indicator. For this model of Slc11a2, RaptorX predicted p-value of 4.55e − 07. Twenty two mutant models were generated in PyMOL (Supplementary Fig. 2, Supplementary Fig. 3, Supplementary Fig. 4, Supplementary Fig. 5, Supplementary Fig. 6, Supplementary Fig. 7, Supplementary Fig. 8, Supplementary Fig. 9, Supplementary Fig. 10, Supplementary Fig. 11, Supplementary Fig. 12, Supplementary Fig. 13, Supplementary Fig. 14, Supplementary Fig. 15, Supplementary Fig. 16, Supplementary Fig. 17, Supplementary Fig. 18, Supplementary Fig. 19, Supplementary Fig. 20, Supplementary Fig. 21, Supplementary Fig. 22, Supplementary Fig. 23).

Ramachandran plot of native protein showed 442 residues (88.8%) in favored region, 42 residues (8.4%) in allowed region and 14 residues (2.8%) in outlier region (Supplementary Fig. 24). While in case of 22 altered proteins, in case of 21 nsSNPs, Ramachandran plot showed similar pattern as native protein but for nsSNP G369R, one residue from favored region was shifted to outlier region. So in this case, 441 residues (88.6%) were in favored region and 15 residues (3.0%) were in outlier region (Supplementary Fig. 25). While checking RMSD value, it was observed that there was not much deviation from native protein. The higher the RMSD value, the more the deviation between the two structures which in turn changes their functional activity [23]. V108I, R465Q, W477L, and V517I showed somewhat higher RMSD values of 0.053, 0.057, 0.033, and 0.036 respectively which are given in Table 3. While tested for polar interactions, in case of some nsSNPs, there is a change in interaction patterns compared to native protein. SNPs V108I, I114T, T336K, T343A, Y374C, Q385H, N492S and A512V showed different interaction patterns than native protein's amino acids. In V108I, V108 formed two polar interactions with K104 and A112 while I108 in altered protein formed three interactions. Two were with K104, I108 with altered bond lengths and third extra interaction with L105 (Supplementary Fig. 26 & B). I114 in case of I114T having interaction with L110 with a bond length of 2.6 Å, but T114 had additional interaction to same amino acid with 3.3 Å bond length due to change of R group from non-polar to polar (Supplementary Fig. 27A & B). K336 which is having positively charged R group in altered protein forming 3.6 Å long interaction with V334, while T334 proposed polar R group not forming any interaction (Supplementary Fig. 28. & B). SNP T343A showed alteration of R group from polar to non-polar, which changed interactions (Supplementary Fig. 29. & B). Y374 in native protein forming one polar interaction with P370 and one with V378 of 2.9 Å and 2.8 Å length respectively. While altered amino acid C374 forming two polar interactions with P370 of 2.9 Å, 3.0 Å lengths and with V378 of 2.8 Å, 4.8 Å lengths (Fig. 1 & B). In Q385C, Q385 formed three polar interactions with T343 of 2.8 Å, 2.9 Å, and 3.2 Å lengths, one with L381 of 2.9 Å length, one with A382 of 3.1 Å length and one with M389 of 3.0 Å length but altered residue H385 had interaction with L381 of 2.9 Å length and one with V389 of 3.0 Å length because of the change in R group (Fig. 2 & B). Native residue N492 showed one interaction with I488 of 3.0 Å length and one with I491 of 3.2 Å length and altered residue S492 forming two polar interactions with I488 of 2.5 Å and 3.0 Å lengths (Fig. 3 & B).Similarly, A512V showed the change in bond length (Fig. 4 & B). While verified further for energy change, T336K, Y374C, Q385H, M389V, R465Q, W477L, L484V, S490F, N492S, D502G, V510M, A512V and V517I showed higher total energy after energy minimization than native protein which are given in Table 3.

Further when analyzed for binding site regions using several tools, results revealed ligands and ligand binding sites which are shown in Supplementary Tables 2 and 3. However, none of above nsSNPs resided in the above identified binding sites.