In silico analysis on the functional and structural impact of Rad50 mutations involved in DNA strand break repair

Multiple sequence alignment analysis and conserved domain analyses

Human Rad50 protein sequence was obtained from National Center for Biotechnology Information (NCBI). The sequence similarity search tool, BLASTP from the NCBI server (http://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to find homologs for Rad50. To investigate the similarity between Rad50 protein in human and other organisms such as Danio rerio, Mus musculus, Rattus norvegicus, Pan troglodytes, Oryctolagus cuniculus, Drosophila melanogaster, Saccharomyces cerevisiae and Schizosaccharomyces pombe, a multiple sequence analysis (MSA) was conducted using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) with default settings to determine consensus and conserved regions between the multiple sequences of different organisms (Sievers & Higgins, 2018). Meanwhile, InterPro (http://www.ebi.ac.uk/interpro/) was used to identify the domains and motifs using human sequence (Finn et al., 2017). InterPro results are classified into several types (families, domains, motif or sites) depending on the biological entity they represent (Finn et al., 2017). Using this tool, Rad50 protein sequence was classified into families and the presence of domains and important sites were predicted. ClustalX software (Thompson, Gibson & Higgins, 2002) was used to view and analyze the conserved regions within the domains and motifs in the selected proteins.

Data mining of Rad50 mutation from literature and SNPs database

Rad50 mutations were identified from previous published manuscripts using PubMed database and their functional impacts were extracted for comparison. Besides that, naturally occurring single nucleotide polymorphisms (SNPs) in Rad50 were retrieved from SNPeffect 4.0 database (http://snpeffect.switchlab.org/about) (De Baets et al., 2012) (date of access: 7 April 2018). SNPeffect 4.0 database currently contains more than 60,000 human SNPs gathered from human avariance list available at UniProt website (https://www.uniprot.org/). It specifically focuses on the molecular characterization, annotation of diseases as well as polymorphism variants in human proteins (De Baets et al., 2012). All these available Rad50 protein mutations (obtained from both literature and databases) have been aligned using pairwise alignment through Clustal Omega between human sequence and other organisms’ sequence, individually. From this analysis, we identified similar mutation sites in human. All the identified equivalent mutations in human were manually refined, for example removing the same residues and mutations that has been studied by several different researchers to identify the non-redundant mutations (Table S1). Identified mutations (after converting to equivalent residues in human) were then mapped into Fig. S1.

Secondary structure prediction and analysis of 3D modeling

The Rad50 templates identified from the BLAST analysis also were used to develop secondary structure and 3D model. The PSIPRED program (http://bioinf.cs.ucl.ac.uk/psipred/) has been utilized for secondary protein structure prediction (Buchan et al., 2013). Secondary structure prediction has revealed a clear distribution of alpha helix, beta sheet and coil in H. sapiens (Helix: 74.69%, coil; 18.29 and beta sheet; 7.01%) (Fig. S2). Databases such as UniProt (https://www.uniprot.org/) and Protein Data Bank (PDB) (https://www.rcsb.org/) were used to identify structural information regarding Rad50 protein in human. Rad50 protein sequence also has been BLAST searched against Protein Data Bank (PDB) sequence in Network Protein Sequence @nalysis (NPS@) (https://npsa-prabi.ibcp.fr/) to identify the most identical structure. The incomplete structure has been further predicted using fold recognition method using Protein Homology/analogY Recognition Engine Version 2.0 (Phyre2) (http://www.sbg.bio.ic.ac.uk/phyre2) (Kelley et al., 2015). Phyre2 is an online tool to predict and analyze protein structure, function and mutations which uses advanced remote homology detection methods to build 3D models, predict ligand binding sites and analyze the effect of amino acid variants (e.g., nonsynonymous SNPs (nsSNPs)) for a protein sequence (Kelley et al., 2015). Rad50 sequence was submitted to the webserver to interpret the secondary and tertiary structures of the model, domain composition and quality. 3D model of Rad50 was run under ‘intensive’ mode that generates a complete full-length model of a protein sequence by using multiple template modeling and simplified ab initio folding simulation (Kelley et al., 2015). UCSF Chimera software was used to view and to analyze the 3D structure (Pettersen et al., 2004).

Prediction of deleterious effects of Rad50 mutations using in silico tools

The Rad50 mutations were in silico predicted using PredictSNP to determine their possible molecular impacts in human (https://loschmidt.chemi.muni.cz/predictsnp1/) (Bendl et al., 2014). Its benchmark dataset contains over 43,000 mutations obtained from the Protein Mutant Database and the UniProt database (Bendl et al., 2014). This tool incorporated six established prediction tools; such as Multivariate Analysis of Protein Polymorphism (MAPP) (Stone & Sidow, 2005), Predictor of human Deleterious Single Nucleotide Polymorphisms (PhD-SNP) (Capriotti & Fariselli, 2017), PolyPhen-1 (Ramensky, 2002), PolyPhen-2 (Adzhubei, Jordan & Sunyaev, 2013), Sorting Intolerant from Tolerant (SIFT) (Sim et al., 2012) and Single-Nucleotide Amplified Polymorphisms (SNAP) (Bromberg & Rost, 2007) to provide a more accurate and robust comparison. We classified the mutations as deleterious if five to seven analyses performed were identified as damaging in PredictSNP. For instance, an in silico prediction was considered accurate when a given mutation predicted to be deleterious (as performed in this study) was also found experimentally deleterious (either in vitro or in vivo with phenotypes such as embryonic lethality, growth defect and/or cancer predisposition) based on previous cited studies. Conversely, the prediction is inaccurate if such deleterious mutations was predicted as neutral or tolerant.

Molecular mechanism of amino acid substitutions

To determine the molecular mechanism based on pathogenicity of amino acid substitutions in Rad50, MutPred2 (Pejaver et al., 2017) (http://mutpred2.mutdb.org/index.html) analysis was carried out. This program predicts the pathogenicity and molecular impacts of amino acid substitutions potentially affecting the phenotype. It is trained on a set of 53,180 pathogenic and 206,946 unlabeled (putatively neutral) variants obtained from the Human Gene Mutation Database (HGMD) (Stenson et al., 2017), SwissVar (Mottaz et al., 2010), dbSNP (Sherry et al., 2001) and inter-species pairwise alignment (Pejaver et al., 2017). The output of MutPred contains a general probability that the amino acid substitution is deleterious/disease-associated, and a list of rank of specific molecular alterations potentially affecting the phenotype with its p-value (<0.05).

Prediction of molecular and structural effects of protein coding variants in Rad50 mutation

Prediction of molecular and structural effects of protein coding variants in Rad50 mutations was performed using SNPeffect4.0 (De Baets et al., 2012) (http://snpeffect.switchlab.org/about). The analysis includes predictions of the aggregation prone regions in a protein sequence (TANGO), amyloid-forming regions (WALTZ) and chaperone binding site (LIMBO). The range of prediction score differences outside −50 to 50 for mutants are considered significant (De Baets et al., 2012). SNPeffect also uses FoldX (Schymkowitz et al., 2005) to analyze the effect of mutations on the structural stability. However, as structure quality is important for the accuracy of delta G predictions for stability, model structures with less than 90% sequence identity to the modeling template structure will not be modeled (De Baets et al., 2012).

Analysis of protein stability

The stability of Rad50 upon single amino acid residue mutations were predicted using MUpro (http://mupro.proteomics.ics.uci.edu/) (Cheng, Randall & Baldi, 2006) and I-Mutant 3.0 (http://gpcr2.biocomp.unibo.it/cgi/predictors/I-Mutant3.0/I-Mutant3.0.cgi) (Capriotti, Fariselli & Casadio, 2005) using default setting, for instance temperature was set at 25 °C and pH 7. Mu-Pro and I-Mutant 3.0 are valuable tools for protein stability prediction and analysis, even when the protein structure is not yet known with atomic resolution. Both use support vector machines (SVM)-based tool to predict protein stability changes for single amino acid mutations either from both sequence or structural information which correctly predicts with over 80% accuracy using cross validation methods (datasets and experimental) (Capriotti, Fariselli & Casadio, 2005; Cheng, Randall & Baldi, 2006). Rad50 protein sequence was searched against the web server and energy changes (ΔΔG) were recorded. Negative value for ΔΔG represents a decrease in protein stability whereas positive value for ΔΔG represents an increase in stability.

Rad50 data acquisition and MSA analysis

Human Rad50 sequence from NCBI database contains 1312aa with the accession number of AAB07119.1. Sequence homology search of the human Rad50 protein was performed against NCBI nonredundant protein databases (E-value ≤ 1E−05) and the result was downloaded for further analysis. Out of 500 sequences, six sequences were choosen for MSA analysis from diverse organisms such as D. rerio, M. musculus, R. norvegicus, P. troglodytes, O. cuniculus, and D. melanogaster. Two sequences, Saccharomyces cerevisiae and Schizosaccharomyces pombe were also included due to widely being used as models in previous Rad50 studies (Table S1).

Analysis of protein domains

Domain identification analysis showed that Rad50 contains three P-loop containing nucleoside triphosphate hydrolase (P-loop NTPase) domains which belong to ATP Binding Cassette (ABC) protein superfamily (De La Rosa & Nelson, 2011). It is located near the N- and C-terminal, at the residue number of 25-103, 130-227 and 1196-1279 (Fig. 1A). Residue annotation showed that Rad50 has six specific motifs including Walker A and Q-loop that are located at the N-terminal whereas Rad50 signature motif, Walker B, D-loop and H-loop/switch region are located at C terminal (Fig. 1A) (De La Rosa & Nelson, 2011). It also has a domain called zinc hook (635-734aa) located at C-terminal region (Fig. 1A) (Hopfner et al., 2002). Multiple sequence alignment (MSA) analysis between human Rad50 and its homologous genes (D. rerio, M. musculus, R. norvegicus, P. troglodytes, O. cuniculus, D. melanogaster, S. cerevisiae and S. pombe) also revealed that these specific motifs are highly conserved (Fig. 1B).

Mutation datasets from the literature and database searches

In order to identify the Rad50 mutations, literature pertaining to the topic was exhaustively searched and 18 articles over the period of 1990 to 2017 were identified. All these mutations from different organisms were listed in Table S1. There are 103 mutations identified which mostly occurred in the protein domains and motifs with various biological effects (Table S2). In order to obtain equivalent mutations in human, pairwise alignment was performed individually between each organism (D. rerio, M. musculus, R. norvegicus, P. troglodytes, O. cuniculus, D. melanogaster, S. cerevisiae and S. pombe) and the Rad50 human sequence as a reference (Table S1). Then, MSA analysis was carried out between these sequences from different organisms (including human) to identify consensus regions (Fig. 1; Fig. S1). Further refinement such as integrating similar mutations that occurred at the same positions (for examples; S1202R, K42R, S679R, P682E, V683R, R1214E, K6E, and R81I) (Table S2) from different organisms of which a total of 80 different mutations or non-redundant mutation were identified. All these mutations have been mapped based on equivalent residues in human (Fig. S1). From SNPeffect 4.0 database, another 13 SNP mutations were also identified (Table S3). However, from the total of 103 mutations obtained from literature, only 42 residues of the Rad50 protein mutations were known to contribute to the most damaging effects in vitro and in vivo such as embryonic lethality (Bender et al., 2002; Roset et al., 2014) and growth defect (Table 1; Table S2) (Alani, Padmore & Kleckner, 1990; Bhaskara et al., 2007; Waltes et al., 2009; He et al., 2012; Barfoot et al., 2015; Hohl et al., 2015). Most of these deleterious mutations reside at the specific motifs such as Walker A, Q-loop, zinc hook, Rad50 signature motif, Walker B and D-loop (Fig. 1B) that become our primary research focus (Fig. 1B).

3D structure modeling of Rad50

Currently there is no complete structure of human Rad50 available. Nonetheless, a crystal structure of Rad50 hook and coil–coil domain (HCC) that contains 182 residues (residue 585-766) has been determined (PDB ID: 5GOX) (Park et al., 2017), which represents 13% of the Rad50 structure in human. We attempted to predict a more complete human Rad50 3D structure model using homology modeling. Homology modeling program Phyre2 successfully predict the N-terminal of 276 residues (2-278) and C-terminal of 155 residues (1153-1306) of human Rad50 with 100% confidence level using the template model (PDB ID: 5DAC) from Chaetomium thermophilum (Seifert, Lammens & Hopfner, 2015) that share 66% sequence identity (Fig. S2). As a result, half of human Rad50 protein structure was obtained. The regions with no 3D structure information available are residues 279-584 and residues 767-1152 (Fig. S2), which were mainly predicted to consist of alpha helices secondary structure (Fig. S2). Pairwise alignment of Rad50 sequence between C. thermophilum (1315aa) and human (1312aa) showed about 30% sequence identity (Fig. S2). The result showed that the partial structure of C. thermophilum that has been determined (black line) are highly conserved with the human sequence (Fig. S2) suggesting that the human structure should also share high structure similarity to C. thermophilum at these regions. In agreement to this, results from Phyre2 prediction showed that the N-terminal and C-terminal of Rad50 form a globular and coil–coil domain, similar to the structure of C. thermophilum (Fig. 2A). With the generated model, six motifs of Rad50 namely Q-loop, Walker B, signature motif, D-loop, Walker A and H-loop were identified and marked in the 3D structure (Fig. 2A). All identified deleterious residues found in the domain were also marked as shown in Fig. 2B. To correlate the deleterious residues in the Rad50 HCC domain with zinc hook motif that was not found in the model, the structure of 181 residues (residue 585-766) that has been determined (PDB ID: 5GOX) independently was employed (Figs. 2C and 2D) for functional analysis.

Analyses of Rad50 mutation deleterious effects

All of these 42 mutations (based on mutations from other organisms mapped to human) (Figs. 1B and 2B) were then analyzed using bioinformatics analyses such as impact of amino acid substitutions (PredictSNP), molecular mechanism (MutPred), structural phenotyping (protein and amyloid aggregation) (SNPeffect 4.0) and protein stability (MuPro and I-Mutant 3.0) (Table 2). All raw data from each analysis has been supplied as Supplemental Data (Table S3 for PredictSNP analysis, Table S4 for MutPred analysis, Table S5 for SNPeffect analysis and finally Table S6 for I-Mutant and MuPro analysis). The results showed that most of the deleterious effects fall into specific motifs such as Walker A, Q-loop, Rad50 signature motif, Walker B and D-loop (Table 2; Fig. 2A). Previous analysis also revealed that mutations at these motifs contributed to a number of biological defects such as growth defect (Alani, Padmore & Kleckner, 1990; He et al., 2012; Hohl et al., 2015) embryonic lethality (Bender et al., 2002; Roset et al., 2014), cancer predisposition (Bender et al., 2002; Roset et al., 2014), hematopoietic and spermatogenic depletion (Bender et al., 2002; Roset et al., 2014) (Table 1). Several mutations at the zinc hook region (C681A, C681S, P682E, C684R and C684S) and ATPase/coiled-coil domain (K6E and K132E) also showed to be deleterious (Table 2; Figs. 2C and 2D).

Furthermore, mutations located at Walker A (Fig. 2A) were predicted to affect catalytic and allosteric site, loss or gain of methylation, alteration of DNA binding, metal binding, ordered interface and the loss of relative solvent accessibility, which all are depending on the types of amino acid substitutions (Table 2). These mutations were predicted to affect ATP binding site motif, N-myristolylation, casein kinase II (CK2), protein kinase A (PKA) phosphorylation site and Forkhead-associated (FHA) functional sites. Mutations at the Walker A region also might led to the decrement of protein stability as predicted by I-Mutant and MuPro. Mutation at the Q-loop region (Q159H) (Figs. 2A and 2B) also predicted to have significant deleterious effect and decreased protein stability, but no effects have been identified on its molecular mechanism and structural phenotyping as predicted by MutPred and SNPeffect 4.0, respectively (Table 2).

Mild deleterious effect was predicted at the mutated zinc hook domain (Table 2; Fig. 2E). Subsequent analysis using MutPred also revealed that any mutation at zinc hook might affect several important functional sites that involved in DNA damage repair signaling response and cell cycle checkpoints such as phosphatidylinositol 3-kinase-related kinases (PIKK) phosphorylation site, protein kinase C (PKC) phosphorylation site and BRCA1 C-terminus (BRCT) phosphopeptide ligands binding sites (Table 2). Moreover, deleterious mutation was predicted at the conserved cysteine residue located at the zinc hook motif (CXXC). For example, amino acid substitutions of alanine (A) and serine (S) at the cysteine residue position 681; (C681(A/S)) (Fig. 2D) may affect N-glycosylation, proline-directed phosphorylation and mitogen-activated protein kinases (MAPK) phosphorylation site, which possibly due to the affected zinc binding domain (Table 2). Another deleterious mutation, C684(R/S) was also predicted to not affect its molecular mechanism but might disrupt diarginine retention/retrieving signal, PKC and PIKK phosphorylation site (Table 2). Whilst P682E (Fig. 2D) mutation may lead to gain of helix, altered coiled coil domain, loss of N-linked glycosylation and CK2 phosphorylation site (Table 2).

Rad50 signature motif (Fig. 2A) is a critical site which could lead to deleterious effects if mutated as suggested by PredictSNP analysis (Table 2). All mutations in this motif (S1202A/R/M, Q1205E and K1206M/A/E) or located near this motif (G1198E, L1211W and R1214A/E) (Figs. 1B and 2B) were predicted to affect the protein allosteric and catalytic sites (Table 2), except for R1198E. Mutations at residue S1202A/R/M (Figs 1B and 2B) might affect PKA phosphorylation sites and glycosaminoglycan attachment site (Table 2). Furthermore, R1214A (Figs. 1B and 2B) mutation might affect ATP-binding cassette, ABC transporter-type, signature and profile functional sites (Table 2). We have also predicted several mutations in Rad50 signature motif such as Q1205E, L1211W and R1214A that contributed to the total defect in the structural phenotyping such as the increment in protein and amyloid aggregation and the decrement of protein stability (Table 2).

We have also predicted K6E, K132E and K105E mutations occurred at the coil-coiled domain or ATPase domain to be deleterious (Figs. 1B and 2B; Table 2). Specifically, the mutations at K6E and K132E might lead to loss of strand or loss of helix, respectively. Additionally mutation at K132E also predicted to affect casein kinase 1 (CK1) and PKC phosphorylation sites (Table 2). Even though K22M and R83I (Fig. 2B) were predicted to be neutral in PredictSNP analysis, both of these mutations have also been predicted to increase protein aggregation tendency (Table 2). The mutation at R83I might contributed to the alteration of coiled coil structure domain, DNA binding and ordered interface, that might affect the functional site such as protein–protein interactions (PPI)-docking motif (Table 2). Another neutral mutation predicted were T191E, C221E and S106E (Figs. 1B and 2B), where T191E mutation might be responsible in altering the coiled coil domain and may affect tumor necrosis factor receptor-associated factor (TRAF), serine/threonine-protein kinase (NEK2) and PKC phosphorylation site (Table 2). On the other hand, C221E and S106E (Figs. 1B and 2B) were predicted to not affect any molecular mechanism or protein aggregation (Table 2).

Consistency between bioinformatics prediction and experimental evidence

PredictSNP was used in this study to provide a more accurate prediction of disease-related mutations as it combines six best performing prediction tools for a consensus classifier (Bendl et al., 2014). Evidently, this in silico analysis was consistent with the results from the previous experimental studies where mutations at the Walker A, D-loop, signature motif, Q-loop and Walker B have shown damaging effects (Tables 1 and 2 ; Fig. 2A).

G41D and K40E (Figs. 1B and 2B) mutations at the Walker A motif (Fig. 2A) and C681A and C684R (Figs. 1B and 2D) mutations at the cysteine residue (CXXC) in the zinc hook motif (Fig. 2C) conferred an identical phenotype with the Rad50 null mutation characterized by total defect in the formation of viable spore in S. cerevisiae experiment (Table 1) (Alani, Padmore & Kleckner, 1990; He et al., 2012). This analysis also identified that mutations at Q-loop (Q159H) and D-loop (D1238N and D1238A) (Figs. 1B, 2A and 2B) were also predicted deleterious (Table 2) and were experimentally shown to interrupt all ATP-dependent activities of the complex in different organisms such as P. furiosus and T bacteriophage respectively (Table 1) (Moncalian et al., 2004; De La Rosa & Nelson, 2011). Furthermore, a E1232Q (Figs. 1B and 2B) mutation at the Walker B motif (Fig. 2A) was also predicted to be deleterious (Table 2). Similarly the mutation of Walker B at residue E798Q in Thermotoga maritima showed low ability to respond to DNA damage (Table 1) (Rojowska et al., 2014). This suggests that this motif is important for a molecular repair process, specifically during DNA binding process, which if mutated will affect the viability of an organism. Our analysis using PredictSNP has identified three mutations, which were N28A (Figs. 1B and 2B) (De La Rosa & Nelson, 2011), D1238H (Figs. 1B and 2B) (De La Rosa & Nelson, 2011) and S1202R (Figs. 1B and 2B) (Kerem et al., 1989; Moncalian et al., 2004) located at the Walker A, D-loop and Rad50 signature motif, respectively (Fig. 2A) (Kerem et al., 1989; Moncalian et al., 2004; De La Rosa & Nelson, 2011).

Mutations at the Walker A domain and Rad50 signature motif (Fig. 2A) may also affect important functional sites such as ATP binding site (Table 2). For example, K42R/M/E/A mutation at the Walker A (Figs. 1B and 2B) in S. cerevisiae and D. radiodurans has been identified experimentally to cause defective in ATPase (Table 1) (Chen et al., 2005; Koroleva et al., 2007) and S793R mutation in Pyrococcus furiosus showed the inhibition of ATP binding and disrupted communication between ATP loops (Table 1) (Moncalian et al., 2004). This mutation further distorted the surface of the C-terminal domain and thus altered the interaction between Rad50 monomers to prevent dimerization (Table 1) (Moncalian et al., 2004). We have also identified mutations at several motifs such as Walker A (G41D, K42M/R/E/A) and Walker B (E1232Q) (Figs. 1B, 2A and 2B) that might affect the binding of FHA phosphopeptide ligands that plays a critical role in DNA damage repair mechanism and cell cycle (Table 2). Many FHA domain–containing proteins localized to the nucleus showed to play a critical role in establishing or maintaining DNA repair, cell cycle checkpoints or transcriptional regulation (Durocher et al., 2000). When mutated, diseases such as Nijmegen breakage syndrome (NBS) and the hereditary cancer syndrome variant Li-Fraumeni (CHK2) will be developed (Matsuura et al., 1998; Varon et al., 1998; Carney et al., 1998; Featherstone & Jackson, 1998; Bell et al., 1999) suggesting the importance of these conserved residues within Rad50 for DNA repair and maintenance.

Mutations at or near the Rad50 signature motif (Figs. 1B and 2A) were also known to be damaging (Table 2), particularly the S1202R (Figs. 1B and 2B) mutation which has been studied the most due to its numerous biological defects in vivo. The same residue mutations of the Rad50 signature motif in yeast (S12025R) and human (S1202R) also generated complexes that were significantly diminished in adenylate kinase (AK) activity that was important for DNA tethering (Bhaskara et al., 2007). Previously, AK deficiency was found to be associated with anemia and several cases of mental retardation and psychomotor impairment (Abrusci et al., 2007), which may explain why disruption of the MRN complex also causes this phenotype on patients (Waltes et al., 2009). In addition, such deleterious mutation also contributed to inviable spores and significant telomere shortening in S. cerevisiae (Bhaskara et al., 2007). Defects in telomere length in human have been known to cause the pathology of several age-related diseases and premature aging syndrome, as well as cancer and other human diseases such as Hoyeraal-Hreidarsson syndrome, Coats plus syndrome, pulmonary fibrosis, dyskeratosis congenita, liver fibrosis and aplastic anemia (Blasco, 2005).

Additionally, most mutations such as G1199E, S1202A/R/M, and Q1205E (Figs. 1B and 2B) at the Rad50 signature motif (Fig. 2A) were identified to affect PKA phosphorylation site (Table 2) suggesting that this site is dependent upon the function of Rad50 signature motif. Phosphorylation is one of the most ubiquitous and important post translational modifications of proteins, and implicated in almost all kinds of cellular processes and pathways (Ptacek & Snyder, 2006). In neurons, enhanced PKA signaling promotes neuronal development, enhances synaptic plasticity, and elevates dopamine synthesis (Dagda & Das Banerjee, 2015). However a deterioration in PKA signaling has contributed to the etiology of several neurodegenerative diseases, such as Alzheimer and Parkinson (Dagda & Das Banerjee, 2015). We hypothesized that the defective PKA functional sites may also lead to Nijmegen breakage syndrome associated with neurological phenotype in Rad50 mutations (Waltes et al., 2009), however this potential phosphorylation sites remain to be validated.

C681A and C684R mutations (Figs. 1B and 2D) at the zinc hook motif (Fig. 2C) were identified deleterious from our analysis (Table 2) and these mutations were known to lead severe defects in various DNA damage response (DDR) such as ataxia-telangiectasia mutated (ATM) protein activation, homologous recombinant, irradiation sensitivity and ataxia telangiectasia and Rad3 related (ATR) protein activation (He et al., 2012). These findings were consistent with our bioinformatics analysis where C684S deleterious mutation at zinc hook (Figs. 1B and 2D) might affect a protein kinase called ataxia telangiectasia mutated (ATM) that belongs to the phosphatidylinositol 3-kinase-related kinase (PIKK) family (Table 2). The ATM protein was known to cause devastating ataxia-telangiectasia syndrome which is characterized by progressive neurological disorder, impaired organ maturation and immunodeficiency (Shiloh & Ziv, 2013). Rad50 phosphorylated ATM at S635 site (Figs. 1B and 2D) of which the mutation on this site showed its importance for cell cycle control signaling and DNA repair mechanism (Gatei et al., 2011).

P682E mutation at the zinc hook motif (Figs. 1B and 2D) was shown to be deleterious (Table 2), where previous study has reported that the double mutation P682E and S679R at the zinc hook motif have reduced zinc affinity and dimerization efficiency leading to mice lethality (Roset et al., 2014). In addition, crossbreeding P682E and S679R mutant mice with wildtype mice produce offsprings with hydrocephalus (accumulation of cerebrospinal fluid within the brain), defects in hematopoietic stem cells and gametogenic cells. This suggests that the hook motif has strong influence on the MRN complex associated with DDR signaling, tissue homeostasis and tumorigenesis, as well as fertility of the organism (Roset et al., 2014). This is consistent with the mutations in the yeast hook domain that has increased chromosomal fragmentation (Cahill & Carney, 2007), suggesting its presence is required for the binding or tethering of chromosomal ends.

Limitations of in silico prediction

Several mutations were functionally predicted to be neutral, in contrast with the previous experimental findings. For example, a few mutations that is, S635G (H. sapiens), S679R, C680N, P682A, V683I (S. cerevisiae), V683R (M. musculus) and Q685S (S. cerevisiae) (Figs. 1B and 2D) located at the zinc hook domain (Table 2) and mutations on K22M (S. cerevisiae and M. musculus), R83I (S. cerevisiae and M. musculus), T191E, C221E and S106E (T. maritima) in the ATPase domain (Figs. 1B and 2B; Table 2) were experimentally validated to be deleterious; some causing embryonic lethality, growth defect, cancer predisposition, as well as hematopoietic and spermatogenic depletion in vivo (Bender et al., 2002). A few previous studies have also shown discrepancies between computer prediction and experimental data. For example, an extensive in silico analysis using PolyPhen2 and MutPred tools of the ATP-binding cassette transporter ABCA1, an important target in anti-atherosclerosis treatment predicted that several nsSNPs can be neutral, contradicting with previous experimental data findings (Marín-Martín et al., 2014). Furthermore, another in silico analysis performed using PolyPhen and SIFT on proteins related to several hereditary diseases such as glucose-6-phosphate dehydrogenase deficiency (G6PD), the receptor 1 for tumor necrosis factor-(TNFRSF1A), and familial mediterranean fever (MEFV) has concluded that some nsSNPs impact may also not be predicted deleterious to correspond to previous phenotypic effect (Tchernitchko, Goossens & Wajcman, 2004). Moreover, in silico identification of PmrAB virulence targets in Salmonella typhimurium also demonstrated false positive prediction when validated experimentally (Marchal et al., 2004) suggesting that more work has to be done to develop a more accurate bioinformatics prediction platforms in the future. In contrast, various SNP prediction software have predicted that these mutations were not damaging (Table 2).

Such discrepancy between the computational prediction and experimental results may be due to several limitations in the bioinformatic tools used in our analysis. Several web-based prediction tools may supply conflicting results (Wan et al., 2008) and even with an integrated predictor, PredictSNP (Bendl et al., 2014), it is also limited by the differences in algorithms, principles, training datasets and information used. For example, MAPP, PANTHER and SIFT in the PredictSNP used alignment scores for functional prediction whereas SNAP, PoplyPhen-1 and PolyPhen-2 used neural network, support vector machine and Naïve Bayes algorithm, respectively (Bendl et al., 2014). Interestingly, we identified that the software predicts most accurately (in agreement with experimental results) for the motifs or sites located at the highly conserved position (Fig. 1B). Conversely, most residues that were predicted to be neutral are located at non-conserved positions in the Rad50 protein (Fig. S3; Table S7), suggesting that these prediction software may have only been trained and preferentially biased towards conserved regions (Gardner et al., 2017). This suggests that computer prediction should also consider and take into account the effect of non-conserved regions outside the motifs/domains too for future improvement in their algorithms. Furthermore, any subsequent prediction studies should also be aware of this limitation (whether located in conserved or non-conserved regions) to carefully deduce the function of their protein mutation of interest.

Nonetheless, we cannot rule out the possibility that these mutations derived from other organisms may not be readily affected or transferred to other organisms including human. This is because certain organisms may possess gene compensation to compensate or mask the effect of such mutations and that the different proteins from different organisms may not have perfectly superimposable function. Hence, future experiments should focus on their validation especially in human cell line studies to better understand the roles of these mutated residues in Rad50 function.