The hZP2 gene data was obtained from Entrez Gene on National Center for Biological Information (NCBI) website and Ensembl genome database. The SNP information for hZP2 (rsIDs, chromosomal position, residue change, and global minor allele frequency (MAF)) was retrieved from the NCBI dbSNP database. Among all the SNPs present in hZP2, the validated ones were cataloged into coding and non-coding. The coding SNPs were further categorized into non-synonymous, synonymous, nonsense, and frameshift. The non-synonymous SNPs were then subjected to a variety of in silico tools as shown in Figure 1.

To filter out the deleterious non-synonymous SNPs (nsSNPs), six different bioinformatics tools were employed. These include SIFT (Sorting Intolerant From Tolerant; https://sift.bii.a-star.edu.sg/) (Ng and Henikoff, 2001; Kumar et al., 2009), PolyPhen-2 (Polymorphism phenotyping v2; http://genetics.bwh.harvard.edu/pph2/) (Adzhubei et al., 2010, 2013), PROVEAN (Protein Variation Effect Analyzer; http://provean.jcvi.org/index.php) (Choi et al., 2012; Choi and Chan, 2015), SNPs&GO (http://snps.biofold.org/snps-and-go/snps-and-go.html) (Calabrese et al., 2009; Capriotti et al., 2013), PhD-SNP (Predictor of human Deleterious Single Nucleotide Polymorphisms; http://snps.biofold.org/phd-snp/phd-snp.html) (Capriotti et al., 2006) and SNAP2 (Screening for Non-Acceptable Polymorphism; https://www.rostlab.org/services/snap/) (Hecht et al., 2015). All these algorithms use different approaches to classify a non-synonymous single nucleotide variation as deleterious or not. The inputs for these were given either in the form of rsIDs or amino acid substitutions (AAS) corresponding to all the nsSNPs in hZP2. The deleterious nsSNPs which were common in at least 3 or 4 algorithms were chosen for further characterization.

MutPred2 (http://mutpred.mutdb.org) was used to predict the structural and functional alterations caused by the deleterious nsSNPs. MutPred2 quantifies the pathogenicity of amino acid substitutions and categorizes them as pathogenic or benign in humans and also predicts their impact on 50 different protein properties (Pejaver et al., 2017). The protein sequence in FASTA format along with the AAS corresponding to the selected nsSNPs was submitted as input in this web server. A p-value threshold of 0.05 and a prediction score ranging between 0 and 1 was used. A higher score reflects a higher probability of pathogenicity and the possible alterations in properties were represented as gain/loss of protein structure and/or function.

I-Mutant2.0   (http://folding.biofold.org/i-mutant/i-mutant2.0.html) predicts the change in stability of a protein, upon single point mutation. It is based on the dataset derived from the ProTherm database which is the most inclusive database of thermodynamic experimental data of free energy changes of protein stability upon mutation under different conditions (Capriotti et al., 2005). hZP2 protein sequence in FASTA format and individual AAS corresponding to selected deleterious nsSNPs were given as input and the corresponding change in stability (Reliability Index; RI) and free energy (Kcal/mol; represented as DDG) was obtained. The RI ranges from 0–10 with 10 being the highest in reliability.

The 3D model of hZP2 protein was obtained using I-TASSER (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) which is the most advanced protein structure prediction server. It uses LOMETS, a multiple threading approach to identify structural templates, and generates a 3D atomic model by comparing it with structurally similar known proteins (Roy et al., 2010). The amino acid sequence of hZP2 protein in FASTA format was given as input. Out of the top five models generated in output, the one with the highest C-score was selected. The model thus selected was viewed in Chimera 1.11 (Pettersen et al., 2004) which allows interactive visualization and analysis of molecular structures and related data. For looking at the effect of SNPs, all the deleterious amino acid changes were substituted manually using its rotamer function, and any new network of contacts or clashes formed were assessed. In addition, 12 mutant models of hZP2 were also generated using I-TASSER by manually substituting the amino acids in FASTA sequence of hZP2 at positions corresponding to the 12 deleterious SNPs.

PROCHECK (https://servicesn.mbi.ucla.edu/PROCHECK/) was used to assess the quality of the 3D model that was generated above. This program gives an evaluation of the overall quality of the structure, based on various stereochemical properties (like phi-psi angles in most favored regions of Ramachandran plot, side-chain parameters, chi1-chi2 plots, etc.) by comparing them with the well-refined protein structures of the same resolution and also highlights the regions that may need further investigation (Laskowski et al., 1993). The selected model of hZP2 in PDB format was submitted as input. The 3D model was also verified by ERRAT (https://servicesn.mbi.ucla.edu/ERRAT/) (Colovos and Yeates, 1993) which compares the statistics of non-bonded interactions of different atoms of the submitted protein model with those of highly refined structures.

TM-Align (https://zhanglab.ccmb.med.umich.edu/TM-align/) was used to calculate TM-score and RMSD values of wild type and mutant models. TM-score tells about the topological similarity between wild type and mutant models and RMSD helps in measuring the average distance between alpha-carbon backbones of wild type and mutant models (Zhang and Skolnick, 2005). The TM-score varies between 0 and 1, with a value of one representing a perfect match between two structures. The RMSD value, on the other hand, represents the deviation of mutant structure from the wild type. A higher RMSD value means greater deviation.

Since the evolutionarily conserved positions in a protein are considered important in terms of its structure and function, the conservation score of all the amino acid positions corresponding to deleterious nsSNPs was calculated using ConSurf (https://consurf.tau.ac.il/) (Ashkenazy et al., 2010). This bioinformatics tool uses PSI-BLAST, CSI-BLAST, or BLAST to find the homologous sequences for the given input sequence and performs multiple sequence alignment using different programs like MAFFT, PRANK, TCOFFEE, MUSCLE, or CLUSTALW and finally gives output in the form of a score that ranges from one to nine where nine represents most conserved and one represents highly variable amino acid position.

Since the glycosylation sites on native human zona are not known, we used prediction software to determine this. Potential N- and O-glycosylation sites in the full-length hZP2 (1–745 aa) were predicted using NetNGlyc-1.0 (http://www.cbs.dtu.dk/services/NetNGlyc/) and NetOGlyc-4.0 (http://www.cbs.dtu.dk/services/NetOGlyc/) respectively. NetNGlyc-1.0 predicts N-linked glycosylation sites in human proteins using artificial neural networks that examine the sequence context of Asn-Xaa-Ser/Thr sequons (Gupta R, 2002). NetOGlyc-4.0 predicts O-GalNAc (mucin type) glycosylation sites in mammalian proteins using neural network predictions (Steentoft et al., 2013). Boja et al., 2003 have characterized the N- and O-linked glycosylation sites in mouse ZP2 (mZP2) using mass spectrometry of native mouse zona proteins. Using this data, manual assertions were made for glycosylation in hZP2 by aligning hZP2 and mZP2 protein sequences using Clustal Omega (Madeira et al., 2019). These manual assertions were compared with those predicted by NetNGlyc-1.0 and NetOGlyc-4.0. In addition, to analyze any deviation in terms of loss or gain of N- or O-glycosylation sites caused by the shortlisted 12 deleterious nsSNPs, FASTA sequences corresponding to the polymorphisms were analyzed using NetNGlyc-1.0 and NetOGlyc 4.0 respectively. This was done to ascertain if a change in N- or O-glycosylations resulted in altered interaction with sperm or, modified the zona structure.

To predict the functional effect of both coding and non-coding SNPs, FuncPred (https://snpinfo.niehs.nih.gov/snpinfo/snpfunc.html) was used. This web-based server selects the SNPs from Genome Wide Association Studies (GWAS) and uses GWAS-SNP p-value data to predict the effect of SNPs on functional characteristics like splice sites, Transcription factor binding sites (TFBS), microRNA binding sites, etc. (Xu and Taylor, 2009). The rsIDs of all the validated nsSNPs (coding and non-coding) in the hZP2 gene were used as input and predictions were obtained.

According to the dbSNP database, a total of 5,152 SNPs were reported in the human ZP2 gene (transcript ID: NM_003460 and protein ID: NP_003451). Out of these 5,152, only 1,069 were found to be validated. Further, among the validated, 229 SNPs were in the coding region and the remaining 840 were in the non-coding region (3’ & 5′ near gene region, 3’ & 5’ UTRs, and introns). Among the 229 coding SNPs, 165 were non-synonymous SNPs (missense; nsSNPs), 56 were synonymous (same-sense; sSNPs), four were nonsense, and four were frameshift (insertions and deletions).

nsSNPs produce amino acid allelic variants of the gene which may affect the structure and function of the protein. Hence, the 165 nsSNPs were selected for further investigations. Of those 22 were predicted to be deleterious by the SIFT web server, with a SIFT score of ≤0.05 (Supplementary Table S1 ). Two of the nsSNPs, rs374388107 and rs267604453, had a score of 0 which is considered the most damaging score. Another program, PolyPhen-2 predicted 64 nsSNPs as probably/possibly damaging. Out of these 64, 30 were marked as “probably damaging” (PolyPhen score between 0.950 and 1). The remaining 34 nsSNPs were designated as “possibly damaging” (PolyPhen score between 0.850 and 0.950) (Supplementary Table S2). PROVEAN web server predicted a total of 51 out of 165 nsSNPs as deleterious, with a PROVEAN score of less than the cut-off value (-2.5) (Supplementary Table S3). According to the predictions by another algorithm SNPs&GO, only 26 nsSNPs are found to be disease-associated as they had a probability value of >0.50 which predicts disease association (Supplementary Table S4). PhD-SNP prediction classified 54 nsSNPs as disease-associated (Supplementary Table S5).

To effectively select the most deleterious nsSNPs and reduce the rate of false-positive predictions, we shortlisted 18 nsSNPs which were commonly classified as deleterious in at least 3 or 4 out of the above mentioned five algorithmic tools by manual concordance (Table 1). These were classified as deleterious nsSNPs. These deleterious nsSNPs were cross-validated for having an effect or being a neutral variant using another web-server called SNAP2, which predicted the functional effects caused by these nsSNPs. All the 18 nsSNPs were predicted as “effect variants” with highly expected accuracy, making these good candidates for further investigation (Supplementary Table S6).

Of the 18 deleterious nsSNPs submitted to MutPred2 web-server, only eight were found to score more than 0.50 and thus, were predicted to result in structural and functional alterations like change in stability of protein, gain or loss of relative solvent accessibility, loss of disulfide linkage, loss of DNA binding sites, loss of strand, an altered transmembrane protein, altered metal-binding site, etc. (Table 2). Apart from MutPred2, I-Mutant predicted the effect of these mutations on the stability of hZP2 protein. A decrease in stability was observed for 16 out of the 18 nsSNPs. The resultant free energy change (kcal/mol) and the reliability index for each of the substitutions are shown in Table 3.

For predicting structural alterations caused due to nsSNPs, first, the 3D model of wild type hZP2 was generated using I-TASSER. Out of the five models generated, the model with the highest C-score of -1.76, an estimated TM-score of 0.50 ± 0.15, and an estimated RMSD of 12.5 ± 4.3Å was used for further analysis ( Figure 2A ). The stereo-chemical quality of the protein model was checked using the PROCHECK program based on various factors like overall G-factor, phi-psi angles, chi1-chi2 plots, side-chain parameters, etc. which were found to be within limits and thus the structure was found acceptable and worth investigating further (Figure 2B). Ramachandran plot showed 58.0% residues in the core region, 32.7% in the allowed region, 6.0% in the generously allowed region, and 3.4% in the disallowed region ( Figure 2C ). ERRAT2 program which verifies the quality of the protein model based on non-bonded interactions predicted the overall quality factor to be 70.21. The generally accepted range for a high-quality model is > 50.

The model generated above was used to study the effect of the mutations on the protein’s 3D structure using Chimera 1.11. Each of the 18 nsSNPs shortlisted above was checked to identify the formation of any new network of contacts and/clashes (Table 4). Out of the 18 nsSNPs, 12 were found to form new networks of clashes and contacts and these were used for further analysis. An example of this is shown in Figure 3 in which, Alanine⁵⁴⁷ formed no network of clashes/contacts in wild type hZP2. However, when it was replaced by Valine, five new pseudo bonds with Val⁶¹¹, Met^595, and Trp⁵⁴⁶ were formed.

Structural models for the 12 nsSNPs which formed new network of clashes and contacts were also generated using I-TASSER. Out of five models obtained in output for each of the 12 mutant proteins, models with the highest C-score were selected for further analysis. Finally, the wild-type and mutant models were compared using TM-align. The TM-Score and RMSD values for each mutant model are shown in Table 5. Mutant models for C155Y (rs761335280) and S384I (rs768663589) were found to have the lowest TM-score i.e. 0.22707 and 0.20968 and a high RMSD value i.e. 8.27 and 8.44 respectively, thus showing a greater deviation from wild type protein. When the conservation score of 12 amino acid positions (i.e. corresponding to 12 deleterious nsSNPs) on the protein was calculated using ConSurf, it was found that six out of 12 nsSNPs (C155Y, G425R, N439K, E440K, A547V, and S627Y) are in highly conserved positions, thus, showing their functional significance (Table 5 ).

NetNGlyc predicted 6 amino acid positions in the wild-type hZP2 protein (N⁸⁷, N¹⁰⁵, N¹²², N²²³, N²⁶⁹, and N⁴⁰⁰). When these were compared with the N glycosylation sites characterized in native mZP2, four of these (N⁸⁷, N²²³, N²⁶⁹, and N⁴⁰⁰) were present in mouse too. The other two positions (N¹⁰⁵ and N¹²²) were specific to hZP2 (Figure 4). Out of the 12 deleterious nsSNPs, none of them was present at a predicted N-glycosylation site or caused any change in the N-glycosylation pattern due to this polymorphism.

NetOGlyc predicted 19 potential O-glycosylation sites (S⁹, S¹¹, S¹³, S⁶³¹, T⁶³³, S⁶³⁷, S⁶³⁸, T⁶⁴⁴, T⁶⁴⁷, S⁶⁵⁵, S⁶⁷⁴, S⁶⁷⁵, S⁶⁸⁰, S⁶⁸², S⁶⁸³, S⁶⁸⁷, S⁶⁸⁹, T⁶⁹¹, and S⁶⁹⁷). In addition, based on manual assertion after comparing with mouse mass spectrophotometric data, Thr⁴⁶² was also assigned to be potentially glycosylated in hZP2 (Figure 4). O-linked glycosylation occurs by transferring oligosaccharides to serine and threonine residues. Among the 12 nsSNPs, only two SNPs (S384I and S627Y) translated into the amino acid substitutions where serine was being replaced. However, none including the above two were present at the predicted O-glycosylation sites. When the corresponding mutant sequence for these SNPs was analyzed by NetOGlyc-4.0, loss of glycosylation sites was observed only for C155Y (S⁶⁵⁵), V382F (S⁶⁷⁴, S⁶⁹¹), N439K (S⁶⁵⁵), E440K (S⁶⁹¹), L531Q (S⁶⁵⁵), A547V (S⁶⁷⁴, S⁶⁸²), P553L (S⁶⁵⁵, S⁶⁷⁴) and S627Y (S⁶³¹, S⁶⁷⁴, S⁶⁹¹). No loss/gain of O-glycosylation was observed for P47H, Q374H, S384I, and G425R.

The 1,069 validated SNPs from both coding and non-coding region were also submitted to FuncPred to detect their role in the regulation of gene expression. Only 31 of these were found to have an effect. Five coding SNPs (rs16971234, rs2075520, rs2075526, rs34159042, and rs35162028) were found to affect splicing [exonic splicing enhancers (ESE) and exonic splicing silencers (ESS)] and the remaining 26 SNPs from the non-coding region were found to affect transcription factor binding sites (TFBS). Also, no SNP in the 3′UTR region was found to create or abolish the miRNA binding site. The detailed results are shown in Table 6. It is interesting to find that most of these regulatory SNPs have high global MAF values (Supplementary Figure S1).