In silico analysis of single nucleotide polymorphism (SNP) in human TNF-α gene.

The TNF-α gene mutations are seen in many diseases especially inflammatory diseases. Hence, before planning a larger population study, it is advisable to sort out the possible functional SNPs. To accomplish this goal, data available in the dbSNP database and different computer programs can be used. Therefore, this study was undertaken to find the functional nsSNPs (non-synonymous single nucleotide polymorphisms) in TNF-α. Out of the total 169 SNPs, 48 were nsSNPs (non-synonymous single nucleotide polymorphisms), 23 occurred in the mRNA 3' UTR, 10 occurred in 5' UTR region, 41 occurred in intronic regions and the rest were other types of SNPs. SIFT and PolyPhen predicted 2 out of 48 nsSNPs as damaging. Among the predicted nsSNPs, rs4645843 and rs1800620 were identified as deleterious and damaging by the SIFT (Sorting Intolerant from Tolerant) and PolyPhen programs. Additionally, I-Mutant and nsSNPAnalyzer showed a decrease in stability for these nsSNPs upon mutation. Protein structural analysis with these amino acid variants was performed by using I-Mutant, Swiss PDB viewer, ANOLEA (Atomic Non-Local Environment Assessment), MUSTER (MUlti-Sources ThreadER) and NOMAD-Ref servers to check their molecular dynamics and energy minimization calculations. This study suggested that P84L and A94T variants of TNF-α could directly or indirectly destabilize the amino acid interactions and hydrogen bond networks thus explaining the functional deviations of protein to some extent.

Introduction

A single nucleotide polymorphism (SNP) is a source variance in a genome. A SNP is a single base mutation in DNA. SNPs are the simplest form and most common source of genetic polymorphism in the human genome (90% of all human DNA polymorphisms). There are several publically available databases for SNPs, such as dbSNP, GWAS Central and SwissVar. Only the non-synonymous SNPs (nsSNPs), also called as missense variants are particularly important as they result in to changes in the translated amino acid residue sequence. It is likely that nsSNPs play a major role in the functional diversity of coded proteins in human populations and have been linked with many diseases. nsSNPs may affect the protein function by reducing protein solubility or by destabilizing protein structure and they may affect gene regulation by altering transcription and translation (Barroso et al., 1999, Chasman and Adams, 2001, Lander, 1996, Smith et al., 1994).

Tumor necrosis factor (TNF) is a cytokine with pleomorphic actions. TNF-α is pivotal in host defense against infections and has a major role in autoimmune diseases as well. It is also a crucial cytokine for granuloma formation. The level of TNF-α varies from individual to individual and is genetically determined (Wilson et al., 1992a, Wilson et al., 1992b). The gene for TNF-α is located within the major histocompatibility complex (MHC) region on chromosome 6p21.3 which is a highly polymorphic region. There are many biallelic single nucleotide polymorphisms (SNPs) in and around the TNF-α gene. TNF-α is a pleiotropic cytokine produced mainly by monocytes/macrophages that are involved in systemic inflammation (Sugimoto et al., 1999). TNF-α exerts cytotoxic effects on renal cells (McCarthy et al., 1998, Min et al., 2009) and it has been shown to participate in renal damage development in experimental models of renal disease including lupus nephritis, glomerulonephritis, nephropathy, hypertension, and diabetes (Elmarakby and Sullivan, 2012, McCarthy et al., 1998). Several investigators have implicated tumor necrosis factor-alpha (TNF-α) cytokine in modulating the progression of ESRD (end stage renal disease) (Klahr, 1999, Sankaran et al., 1999).

Several polymorphisms have been identified in the upstream regulatory region of TNF-α. Among these are common biallelic G (TNF1 allele) transition to A (TNF2 allele) 308 (rs1800629) nucleotide upstream from the transcription initiation site in the TNF-α promoter. This transition is considered to be an important enhancer of transcriptional activation associated with elevated levels of TNF-α (Wilson et al., 1997a, Wilson et al., 1997b), which have been shown to be involved in increased susceptibility to different eye diseases including diabetic retinopathy and glaucoma (Huang et al., 2006, Limb et al., 1999, Yoshioka et al., 2006). TNF-α-308 promoter gene polymorphism has been reported to be associated with several autoimmune disorders including systemic lupus erythematosus, rheumatoid arthritis and infections such as tuberculosis (Hajeer and Hutchinson, 2000). Taking into account all these considerations and the central role played by TNF-α in many diseases, in the present study we aimed to determine the influence of various polymorphisms in TNF-α on its protein structure that may have an important role in disease susceptibility.

Materials and methods

The data on human TNF-α gene was collected from Entrez Gene on National Center for Biological Information (NCBI) web site. The SNP information (protein accession number and SNP ID) of the TNF-α gene was retrieved from the NCBI dbSNP (http://www.ncbi.nlm.nih.gov/snp/) and SwissProt databases (http://expasy.org/).

Evaluation of the functional impact of coding nsSNPs using a sequence homology tool (SIFT)

SIFT takes a query sequence and uses multiple alignment information to predict tolerated and deleterious substitutions for every position of the query sequence (Ng and Henikoff, 2003). It is a multistep procedure that, given a protein sequence, (1) searches for similar sequences, (2) chooses closely related sequences that may share similar function, (3) obtains the multiple alignment of these chosen sequences, and (4) calculates normalized probabilities for all possible substitutions at each position from the alignment. Substitutions at each position with normalized probabilities less than a tolerance index of 0.05 are predicted to be intolerant or deleterious; those greater than or equal to 0.05 are predicted to be tolerated (Ng and Henikoff, 2001, Ng and Henikoff, 2006).

The analysis was performed by allowing the algorithm to search for homologous sequences using the default settings (UniProt-TrEMBL 39.6 database, median conservation of sequences of 3.00, and allowance to remove sequences more than 90% identical to query sequence). The TNF-α FASTA amino acid sequence of the NCBI protein accession id CAA26669 was used as the query sequence and a total of 169 TNF-α nsSNPs filtered from the dbSNP database were analyzed.

PolyPhen

PolyPhen (version 2) (http://genetics.bwh.harvard.edu/pph2/) predicts the influence of amino acid substitution on the structure and function of proteins by using the specific empirical rules. Protein sequence, database ID/accession number, amino acid position and amino acid variant details are the input options for PolyPhen (Ramensky et al., 2002). The tool estimates the position-specific independent count (PSIC) score for every variant and calculates the score difference between variants.

nsSNPAnalyzer

It uses a machine learning method called Random Forest to classify the nsSNPs. It was trained using a curated SNP dataset prepared from the SwissProt database. nsSNPAnalyzer calculates the three types of information from the user's input: 1) the structural environment of the SNP, including the solvent accessibility, environmental polarity and secondary structure (Bowie et al., 1991); 2) the normalized probability of the substitution in the multiple sequence alignment (Ng and Henikoff, 2001); and 3) the similarity and dissimilarity between the original amino acid and mutated amino acid.

The minimum input includes a protein sequence in FASTA format and a substitution file denoting the SNP identities to be analyzed. Substitution file format: a substitution is denoted as X#Y where X is the original amino acid in one letter, # is the position of the substitution (starting from 1) and Y is the mutated amino acid in one letter. Multiple substitutions should be separated by new line characters.

Identification of cis regulatory elements

PROSCAN version 1.7 Web Promoter Scan Service (http://bimas.dcrt.nih.gov/molbio/proscan/)

PROSCAN version 1.7 Web Promoter Scan Service predicts promoter regions based on homologies with putative eukaryotic Pol II promoter sequences. The site is serviced and maintained by Dr. Dan Prestridge at the Advanced Biosciences Computing Center, University of Minnesota.

Promoter 2.0 Prediction Server (http://www.cbs.dtu.dk/services/Promoter/)

Promoter 2.0 predicts transcription start sites of vertebrate Pol II promoters in DNA sequences. It has been developed as a frequently updated database of simulated transcription factors that interact with sequences in promoter regions. It builds on principles that are common to neural networks and genetic algorithms. The site is serviced and maintained by Steen Knudsen at The Center for Biological Sequence Analysis at the Technical University of Denmark.

TSSG

Recognition of human PolII promoter regions and transcription start sites from Softberry (http://www.softberry.com/) — TSSG is the most accurate mammalian cis element prediction program.

Modeling nsSNP locations on protein structure

I-Mutant and FOLD-X

I-Mutant (version 2.0) (http://folding.uib.es/i-mutant/i-mutant2.0.html) is a neural network based tool for the routine analysis of protein stability and alterations by taking into account the single-site mutations. The FASTA sequence of protein retrieved from UniProt is used as an input to predict the mutational effect on protein stability. I-Mutant also provides the scores for free energy alterations, calculated with the FOLD-X energy based web server. FOLD-X is a computer algorithm for quantitative estimation of interactions facilitating the stability of proteins. The FOLD-X tool provides the comparison between wild type and mutant models in the form of van der Waals clashes, which greatly influence the energy decomposition. Sometimes the mutations can cause and reduce strain in the original native structure (Abagyan and Totrov, 1994, Schymkowitz et al., 2005).

Models generated by MUSTER

MUSTER is a MUlti-Source ThreadER program, which considers six different sources: (1) sequence-derived profiles; (2) secondary structures; (3) structured-derived profiles; (4) solvent accessibility; (5) torsion angles (psi and phi angles); and (6) hydrophobic scoring matrix. The optimized threading is found by global dynamic programming. MUSTER provides the Z-score and complete full length models by using MODELLER v8.2. If the calculated Z-score is greater than 7.5, the corresponding template is considered good otherwise designated as bad (Wu and Zhang, 2008).

Results

The dbSNP database contains both validated and non-validated polymorphisms. In spite of this drawback, we opted to avail the dbSNP because the allelic frequency of most of nsSNPs of TNF-α has been recorded there and that is the most extensive SNP database (Bhagwat, 2010). It contained a total of 169 SNPs, out of which 48 were nsSNPs, 23 occurred in the mRNA 3′ UTR, 10 occurred in 5′ UTR region, 41 occurred in intronic regions and the rest were other types of SNPs as shown in Fig. 1. We selected missense nsSNPs, 3′ UTR SNPs, and intronic SNPs for our investigation.

Fig. 1

A graphical representation of distribution of nonsynonymous, 5′ UTR, 3′ UTR and intronic SNPs for TNF-α gene (based on the dbSNP database).

Deleterious nsSNPs by SIFT program

Protein sequence with mutational position and amino acid residue variants associated to 24 missense nsSNPs were submitted as input to the SIFT server, and the results are shown in Table 1, along with the corresponding heterozygosity and validation status description for each SNP, when available from dbSNP.

Table 1

List of nsSNPs that were analyzed by SIFT.

SNP	Amino acid change	Protein ID	Amino acid	Prediction	Score
rs1800620	A94T	NP_000585	A	Tolerated	1
T	Damaging	0.02
rs3179060	H52N	NP_000585	H	Tolerated	1
N	Tolerated	0.1
rs4645843	P84L	NP_000585	P	Tolerated	1
L	Tolerated	0.12
rs11574936	I194N	NP_000585	I	Tolerated	1
N	Damaging	0
rs35131721	P64L	NP_000585	P	Tolerated	1
L	Tolerated	0.06
rs104895105	H478Y	NP_000234	H	Tolerated	0.56
Y	Tolerated	0.09

Protein sequence with mutational position and amino acid residue variants associated to nsSNPs were submitted as input to the SIFT server, and the results are shown in Table 1. According to the classification proposed by Ng and Henikoff 2003 and the lower the tolerance index, the higher the functional impact a particular amino acid residue substitution is likely to have and vice versa.

Among the total nsSNPs analyzed, 2 nsSNPs were identified to be deleterious with a tolerance index score ≤ 0.05. One nsSNP (rs11574936) showed a highly deleterious tolerance index score of 0.00. The remaining nsSNPs (rs3179060, rs4645843, rs35131721 and rs104895105) were found to be tolerated with tolerance index scores of ≥ 0.05.

Damaged nsSNPs by PolyPhen server

All the 6 nsSNPs submitted to SIFT were also submitted to the PolyPhen server. A PSIC score difference of 1.5 and above is considered to be damaging. Two nsSNPs (rs4645843, rs1800620) were considered to be damaging and exhibited a range of PSIC score difference between 1.503 and 2.609. SNP rs1800620 was also found damaging by SIFT server.

nsSNPAnalyzer is a tool to predict whether a nonsynonymous single nucleotide polymorphism (nsSNP) has a phenotypic effect. nsSNPAnalyzer also provides additional useful information about the SNP to facilitate the interpretation of results, e.g., structural environment and multiple sequence alignment. nsSNPAnalyzer uses information contained in the multiple sequence alignment and information contained in the three dimensional protein structure to make predictions.

Two amino acid variants (P84L and A94T) were submitted to nsSNPAnalyzer tool along with protein sequence of TNF-α in FASTA format. We could not find any significant result for the P84L variant. The output of final result of nsSNPAnalyzer tool is shown in Table 2.

Table 2

Output for nsSNPAnalyzer.

Amino acid variant	Phenotype	Environment	Area buried	Frac polar	Secondary structure
A94T	Neutral	P1S	0.303	0.115	S

Amino acid variant was found to be neutral.

After analyzing the coding regions, we focused our attention on non-coding variants also. Specifically, we analyzed the genotypic data and haplotypic information of Phase 1 of the 1000 Genomes Project. We found 331 non-coding genetic variants. To evaluate the impact of these variants on TNF-α gene, we performed different in silico analyses like SIFT and nsSNPAnalyzer. 36 non-coding SNPs were found to be affected and among them 9 SNPs were observed deleterious by SIFT. When we enroll the same non-coding SNPs in nsSNPAnalyzer, 3 SNPs were found to be affected (Renato et al., 2014). Among these SNPs, I194N was found to be of disease type.

Analysis of cis regulatory elements

PROSCAN: version 1.7

Cis element region was predicted on forward strand at 197 bp to 447 bp.

Softberry TSSG

1 promoter is predicted.

Promoter position: 469 LDF.

Promoter 2.0 Prediction Server

Structural analysis of mutant structures

Out of 6 nsSNPs predicted to be deleterious by SIFT or PolyPhen, 2 (rs4645843, rs1800620) were mapped to the native structure by I mutant 2.0 server. Then, energy minimizations were performed by GROMACS for the native structure and the mutant modeled structures (Table 4).

Table 4

Total energy of native and mutant structures after energy minimization.

Amino acid variants	Total energy after minimization (kJ/mol)	Electrostatic constraint
Native	− 6807.950	− 4973.34
P84L	− 8366.239	− 5316.41
A94T	− 8601.871	− 5441.50

Total energy and electrostatic constraint were found to be negative.

Prediction of protein structural stability

I-Mutant is a neural network based routine tool used in the analysis of protein stability alterations by considering the single-site mutation. I-Mutant also provides the scores for free energy alterations, calculated with the FOLD-X energy based web server. By assimilating the FOLD-X estimations with those of I-Mutant, the 93% precision can achieved. The two mutations (84, P → L and 94, A → T) of TNF α gene have been selected on the basis of prediction scores of PolyPhen. These variants were given to I-Mutant web server to predict the DDG stability and reliability index (RI) upon mutation (Table 3).

Table 3

Protein structural stability based on standard free energy change.

Mutation	Position	WT	New	pH	Temperature	Stability	DDG(kcal/mol)
84: P–L	84	P	L	7.0	25 °C	Decrease	− 0.29
94: A–T	94	A	T	7.0	25 °C	Decrease	− 1.02

Where, “WT” is the amino acid in native protein, “New” is mutant amino acid and DDG is the stability (DDG < 0: decrease stability, DDG > 0: increase stability).

Modeling of mutant structure

Information about mapping the deleterious nsSNPs into protein structure was obtained from dbSNP. The available structure for TNF alpha gene has a PDB id 4G3Y. Two nsSNPs were found to be the highest deleterious nature among all the nsSNPs. Hence we selected these nsSNPs for structural analysis. The mutational position and amino acid variant associated with this nsSNP is P → L at the residue position 84 and at the residue 94 A → T were mapped by ANOLEA to get modeled structure as shown in Fig. 2.

Fig. 2

(a) Native structure showing phenyl alanine and arginine at positions 84 and 94 respectively. (b) Mutant modeled structure showing tyrosine residue at position 94; deep view of superimposed structure of wild and mutant residue at position 94. (c) Mutant modeled structure showing leucine residue at position 84; deep view of superimposed structure of wild and mutant residue at position 84.

MUSTER

If Z-score is greater than 7.5, the corresponding template is considered good otherwise designated as bad. It was found that for all the alignments of TNF-α Z score was > 7.5, which indicates that all the templates can be considered as good type (Table 5).

Table 5

Z score value of different templates analyzed by MUSTER.

Rank	Template	Align_length	Coverage	Z score	Seq_id	Type
1	3it8C	152	0.652	14.246	0.993	Good
2	1tnrA	140	0.6	12.032	0.35	Good
3	4msvA	137	0.587	11.587	0.299	Good
4	2re9A	146	0.626	11.445	0.329	Good
5	3ugnA	140	0.6	11.173	0.307	Good
6	1i9rA	140	0.6	10.783	0.264	Good
7	1s55A	141	0.605	10.496	0.227	Good
8	1d4vB	142	0.609	9.9	0.232	Good
9	4mxwY	127	0.545	9.494	0.291	Good

Different templates were found based on alignment score.

Discussion

Tumor necrosis factor (TNF) is a cytokine which plays an important role in many diseases like diabetes, arthritis, and multiple myeloma. TNF-α is pivotal in host defense against infections and has a major role in autoimmune diseases as well. The level of TNF-α varies from individual to individual and is genetically determined. The gene for TNF-α is located within the major histocompatibility complex (MHC) region on chromosome 6p21.3 which is a highly polymorphic region (Wilson et al., 1997a, Wilson et al., 1997b). It has been observed that TNF enhancer polymorphism is implicated in several diseases, and TNF rs1800629 and rs361525 SNPs are the most important in human disease susceptibility as these might influence the transcription of TNF gene. It has been found that TNF enhancer polymorphism influences the serum level of TNF in different human diseases and thus affects the susceptibility to diseases. The presence of DNA sequence variation in TNF gene causes the modification of transcriptional regulation and is thus responsible for the association of susceptibility/resistance with human diseases (Qidwai and Khan, 2011).

Many human SNPs that are now recognized (in excess of 4-million unique SNPs) (http://www.ncbi.nlm.nih.gov/SNP/index.html), along with the genome sequence and other proteome information, provide an opportunity for a much broader understanding of the association between genotype and phenotype. Therefore, an effort was made to identify SNPs that can modify the structure, function and expression of the TNF-α gene. There were 48 nsSNPs submitted to the SIFT as well as to the PolyPhen server, out of these 48 nsSNPs, 6 nsSNPs were found to be damaging by SIFT and 2 nsSNPs were found to be deleterious by PolyPhen. Associations between polymorphism (rs1799964, rs1799724, rs1800630) and immune-mediated diseases such as rheumatoid arthritis and Crohn's disease (CD) have been reported. It has been found that in patients with diabetic nephropathy, TNF variant like rs1800629 was more frequent than in non-diabetic patients with chronic renal failure (Wilson et al., 1997a, Wilson et al., 1997b). There was not any study reported which shows association between damaging nsSNPs (rs4645843 and rs1800620) and any disease. Therefore the validation of these nsSNPs in any disease is required to complement this finding. Similar approach has been used to investigate the effect of nsSNP of BARD1 gene and BRCA1 gene (Ali, et al., 2012, Rajasekaran et al., 2007). The combination of the analysis of human genetic variations of the TNF-α gene, together with the computational method to predict their possible functional impact, can facilitate the analysis of TNF-α gene variant and their effects on protein functional characteristics. Specifically, this approach allows us to estimate the probability of an amino acid change to be deleterious for TNF-α function (Rosanna et al., 2014). In conclusion, the presence of a non-complete concordance among clinical evidences, experimental investigations, and in silico analyses strongly suggested that improvements are necessary in all these fields, in order to enhance our understanding of TNF-α nsSNPs and their role in the pathogenesis of the TNF-α related disease.

Molecular dynamics (MD) simulation was performed to study the explicit solvent behavior of native and mutant structures to examine the difference in dynamics and stability of native and fetal mutations. We compared RMSD value and total energy values (kcal/mol) of native structure and mutated modeled structure for TNF-α gene variant. Two variants P84L and A84T were found to have decreased the stability of protein structure. This may be due to mutant residue which is bigger than the wild type and cannot fit within the available space. Due to carrying the less rigid behavior, the mutant residue can possibly disturb the original core structure of native protein. The two identified TNF variants i.e. P84L and A84T could de-regulate the inflammatory signaling cascade, thus may participate in human inflammatory diseases. Hence, P84L and A84T variants constitute a unique resource of genetic markers that may considerably increase the power of TNF-α gene mutation-screening in disease epidemiological studies.

Conclusion

Functional and structural impact of SNPs in the TNF-α was found out gene using computational prediction tools. Out of a total of 169 SNPs in the TNF-α gene, 48 were nsSNPs, 23 occurred in the mRNA 3′ UTR, 10 occurred in 5′ UTR region, 41 occurred in intronic regions and the rest were other types of SNPs. Out of 48 missense nsSNPs, 2 were found to be deleterious by SIFT, and 2 were found to be damaging by the PolyPhen tool. Structural analysis results showed that the amino acid residue substitutions which had the greatest impact on the stability of the TNF-α protein were mutations P84L (rs4645843) and A94T (rs1800620). Based on our results, we conclude that these SNPs should be considered important candidates in causing diseases related to TNF-α gene malfunction.

Position (bp)	Score	Likelihood
700	1.077	Highly likely prediction
1500	0.665	Marginal prediction
3200	0.577	Marginal prediction