Forecasting most deleterious nsSNPs in human TLR9 gene and their cumulative impact on biophysical features of the protein using in silico approaches

Abstract

In women, the uterine cervix and corpus uteri are two main suspects, playing a major role in cancer-associated-mortality. Immunologically, Toll-like receptors (TLRs) associated with the innate immune system, can recognize pathogens and induce immune responses against pathogens. Cellularly, TLR9 expression occurs in immune system cells including macrophages, natural killer cells, dendritic cells, and other antigen-presenting cells. TLR9 recognizes and interacts with viral and bacterial DNA comprising cytosine-phosphate-guanine (CpG) dideoxynucleotide motif. The current study is designed to identify the most deleterious nonsynonymous single nucleotide polymorphisms (nsSNPs) in the TLR9 gene and to delineate their deleterious effect on the structural and functional features of proteins at the molecular level. Based on the implementation of various computational tools and algorithms eight most deleterious nsSNPs (P139H, R257C, C265Y, L283P, G514D, L544Q, H566Y, and W670R) have been identified in the human TLR9 gene as potentially damaging SNPs. Further, our study suggests highly conserved patterns at deleterious nsSNPs sites could influence protein stability and its functional features. Additionally, this study identifies two nsSNPs (G514D and W670R) associated with the severity of Uterine corpus endometrial carcinoma. In support of our computational findings, the validation of key results using polymerase chain reaction and other experimental methods is warranted in the Indian population. In general, this study might be able to delineate the guideline for identifying the most damaging SNPs and enhances the understating of the risk factors for cancer and disease susceptibilities.

Keywords:

• TLR9
• missense SNPs
• uterine cervix
• Corpus uteri
• phenotype

Introduction

Worldwide, in gynecological cancers, uterine cervix and corpus uteri (endometrial cancer) play a major role in cancer-associated mortality. In India, enhancement in cervical cancer and corpus uteri cancer cases has been seen with time (Ferlay et al. 2021). Recent studies indicate that the involvement of immunological responses plays a crucial role in cancer development (Singhal et al. 2015; Bandil et al. 2017). Innate immunity works as the first line of defense against antigens i.e., viruses and bacteria. Innate immune cells need to be activated in order to perform effective immune responses. Subsequently, innate immune cells rely on various types of pattern recognition receptors (PRRs) expressed intracellularly or on the cellular surface. Functionally, different PRRs recognize unique structures known as pattern-associated molecular patterns (PAMPs) located on antigen’s surfaces (Pradere et al. 2014). One such PRR is Toll-like receptors (TLRs) that comprise a wide range of recognition sensors to initiate the immunological response against the pathogens for quick production of inflammatory cytokines and activation of adaptive immunity (Hemmi et al. 2000).

In homo sapiens, based on primary sequence, TLRs are classified into 10 subclasses likewise TLR1–TLR10 (Imler and Hoffmann 2001). Physiologically, TLR1, TLR2, TLR4, TLR5, and TLR6 are localized at the plasma membrane, crucially identifying bacterial and mycobacterial populations. Cellularly, TLR3, and TLR7 are allocated in endolysosomes and play a key role in the detection of dsRNA and ssRNA viruses respectively (Kawai and Akira 2009). Additionally, TLR9 is also localized into the endolysosome and detects cytosine phosphate-guanine oligodeoxynucleotide (CpG-ODN) of bacteria and viruses (Bauer et al. 2001; Kumar 2018).

In humans, TLR9 is expressed in mast cells, platelets, T cells, innate lymphoid cells (ILCs), NK cells, neutrophils, epithelial and endothelial cells (Bauer et al. 2001; Krug et al. 2004; Akira et al. 2006; Pradere et al. 2014). The CpG motifs comprise a hexamer (RRCGYY) with centrally unmethylated CpG, where R represents purine and Y represents pyrimidine (Park et al. 2008). Biologically, the immunostimulatory activity of the CpG-DNA binding motif is influenced by the density of CpG motifs, its length, and amino acid sequences flanking the CpG motif (Krieg et al. 1995). TLR9 also comprises long insertion loop (Z-loop), where proteolytic processing occurs to perform its functional activity (Ewald et al. 2008; Onji et al. 2013). Recent studies show an association of the TLR9 in the progression of various cancer types like interactions between HPV and host immune response facilitated by TLR9 which influences the transformation to cervical cancer malignancy (Janeway and Medzhitov 2002; Werner et al. 2012; Sharma et al. 2019; Moura et al. 2022). Clinically, TLR9 is highly polymorphic and represents a crucial drug target to explicate the critical role in the development of endometrial cancer (Krieg 2008; Luquez et al. 2021).

The presence of single nucleotide polymorphism (SNP) in the TLR9 gene is responsible for amino acid substitutions causing functional changes or deleterious effects caused by respective SNP. The improper function of TLR9 could lead to inappropriate activation of the downstream immune response pathway (Lamphier et al. 2006). The deleterious nsSNPs in TLR9 could also change its structural features or the mutations in the promoter region disturb gene regulation (Oubounyt et al. 2019). The alterations in gene sequence subsequently affect the alteration of protein function such as protein charge, hydrophobicity, translation, stability, and inter/intra protein interactions leading to serious disease (Kamaraj and Purohit 2013; Emadi et al. 2020; Rozario et al. 2021). Accordingly, the current study is designed to identify the most deleterious nsSNPs in TLR9 and delineate their deleterious effect on the structural and functional features of proteins at the molecular level. In respect, several computational methods and tools have been used to achieve the goal of this study. Subsequently, based on the scientific background and medical importance of TLR9, the key findings of the current study are discussed with the vital role of deleterious nsSNPs in respect of cervical and endometrial cancer.

Results and discussion

Identification of most deleterious nsSNPs in human TLR9 gene sequence using various computational methods

In the TCGA database, of the total 263 SNPs associated with the human TLR9 gene sequence, 21 SNPs are 5′-UTR (8.0%), 7 SNPs are 3′-UTR (2.7%), 156 SNPs are missense (59.3%), and 79 SNPs are others (30.0%) (Figure 1). In others, 67 SNPs are synonymous, 7 SNPs are stop gained, and 5 SNPs are frameshift (Figure 1). Here, we have performed a consensus method with SIFT, SNAP2, Align GVGD, PolyPhen-2, PROVEAN, PhD-SNP, and p-MUT to identify the most deleterious nsSNPs from the data set of missense SNPs. Accordingly, the reported 156 missense SNPs data set has been used with every individual SNPs prediction method, and their altogether results have been manually analyzed and compared (Supplementary Table S1).

Figure 1. Three-dimensional pie-chart analysis of human TLR9. The figure shows the association of percentage missense SNPs, 5′-UTR, 3′-UTR, and other types of SNPs. The data set of correlated TLR9 SNPs is retrieved from the TCGA database. The total of retrieved 263 SNPs are comprising the human TLR9 gene sequence, where 21 SNPs are 5′-UTR (8.0%), 7 SNPs are 3′-UTR (2.7%), 156 SNPs are missense (59.3%), and 79 SNPs are others (30.0%).

Table 1. Analyzing the most deleterious nsSNPs of human TLR9 gene using various bioinformatics tools.

SNP ID	Substitutions	Codons	Sift		SNAP2		Align-GVGD	PolyPhen		PROVEAN		Phd-SNP		P-Mut
			PREDICTION	SCORE	PREDICTION	SCORE	PREDICTION	PREDICTION	SCORE	PREDICTION	SCORE	PREDICTION	SCORE	PREDICTION
I	C265F	TGT→TTT	Deleterious	1	Effect	67	C65	Probably damaging	1	Deleterious	−10.45	Disease	1	Disease
		TGC→TTC
II	G514D	GGT→GAT	Deleterious	1	Effect	54	C65	Probably damaging	1	Deleterious	−5.84	Disease	6	Disease
		GGC→GAC
III	H566Y	CAT→TAT	Deleterious	1	Effect	76	C65	Probably damaging	1	Deleterious	−5.44	Disease	1	Disease
		CAC→TAC
IV	L283P	CTT→CCT	Deleterious	1	Effect	70	C65	Probably damaging	1	Deleterious	−4.39	Disease	1	Disease
		CTG→CCG
V	L544Q	CTA→CAA	Deleterious	1	Effect	59	C65	Probably damaging	1	Deleterious	−4.86	Disease	3	Disease
		CTG→CAG
VI	P139H	CCC→CAC	Deleterious	1	Effect	56	C65	Probably damaging	1	Deleterious	−8.22	Disease	7	Disease
		CCT→CAT
VII	R257C	CGT→TGT	Deleterious	1	Effect	69	C65	Probably damaging	1	Deleterious	−7.61	Disease	2	Disease
		CAC→TAC
VIII	W670R	CAC→TAC	Deleterious	1	Effect	87	C65	Probably damaging	1	Deleterious	−12.72	Disease	1	Disease

Here, we have shown only most deleterious nsSNPs which have been identified as common by every individual tool. Detailed mutation IDs have been shown in Supplementary Table S1.

Sift probability score: <0.05 =t deleterious and score ≥0.05 = tolerated.

SNAP2 probability score: neutral: −100 ≤ score ≤ 0 and effect: 0< SNAP2 score≤ 100.

Align GVGD probability score: C15 = least likely and C65= most likely, effect protein function.

PolyPhen-2 probability score: 0= Neutral, 0 < Score > 0.9= possibly damaging or 1.0= probably damaging.

PROVEAN probability score: ≤ 2.5= deleterious and otherwise = neutral.

PhD-SNP probability score: >0.5= Pathogenic and otherwise = neutral.

PMut probability score: 0–0.5= neutral and 0.5–1= pathological.

Of the total 156 missense SNPs, SIFT tool identifies 62 missense SNPs as the most deleterious nsSNPs. Similarly, 40, 22, 79, 61, 63, and 22 nsSNPs have been predicted as the most deleterious nsSNPs by using SNAP2, p-LIGN-GVGD, PLYPHEN, PROVEAN, PhD-SNP, and p-MUT techniques respectively (Supplementary Table S1). Our comparative outcomes suggest that of the total 156 missense SNPs, only 8 are identified as common nsSNPs by all the techniques (Table 1). Therefore, we could argue that only eight nsSNPs could be the most deleterious and might be associated with the establishment of various cancer types and diseases (Table 1). With respect to the TLR9 protein sequence, the eight identified deleterious nsSNPs have been labeled and are named as P139H, R257C, C265F, L283P, G514D, L544Q, H566Y, and W670R (Table 1). Since the identified most deleterious sites have been replaced by the distinct nature of amino acid residues the changes at these sites might influence the structural and functional features of the TLR9 protein.

Analyzing the effect of most deleterious nsSNPs on the structural stability of TLR9 protein

The present study is aimed to identify the most deleterious nsSNPs of the TLR9 gene sequence and delineate their cumulative role in the structural and functional parameters of the protein. Accordingly, the effect of the most deleterious nsSNPs on TLR9 protein stability has been analyzed by the DynaMut2, SDM, and I-MUTANT Suite servers. The findings of the DynaMut2 and I-MUTANT server revealed DDG < 0 for seven deleterious nsSNPs (P139H, R257C, C265Y, L283P, G514D, L544Q, and W670R) suggesting destabilizing the TLR9 protein (Table 2) (Rodrigues et al. 2018). However, the reverse has been seen for H566Y SNP with positive DDG = 0.41 and 0.37 respectively by both the servers, suggesting an enhancement in the protein stability (Table 2). Further, SDM server analysis for all deleterious SNPs revealed negative DDG < 0 values and a reduction in the protein stability (Table 2). Here, the comparative outcomes of DynaMut2, SDM, and I-MUTANT servers analysis suggest that point mutations might affect the conformational changes that lead to altering the protein stability.

Table 2. Identifying the effect of deleterious nsSNPs on protein stability and functional features of TLR9.

Amino acid residues	I - Mutant suite		SDM		DynaMut2
	Stability	DDG value (kcal /mol)	Stability	DDG value (kcal /mol)	Stability	DDG value (kcal /mol)
P139H	Decrease	−1.73	Decrease	−1	Destabilizing	−1.9
R257C	Decrease	−0.98	Decrease	−1.65	Destabilizing	−1.03
C265Y	Decrease	−0.14	Decrease	−1.03	Destabilizing	−1.16
L283P	Decrease	−1.37	Decrease	−2.81	Destabilizing	−1.74
G514D	Decrease	−0.87	Decrease	−2.45	Destabilizing	−1.53
L544Q	Decrease	−1.83	Decrease	−1.55	Destabilizing	−2.3
H566Y	Increase	0.37	Decrease	−0.31	Stabilizing	0.41
W670R	Decrease	−1.01	Decrease	−0.5	Destabilizing	−2.51

DDG value = free energy change

>−0.5 kcal/mol = largely unstable

<0.5 kcal/mol = largely stable

if value in between = neutral.

Conservation profile analysis of most deleterious nsSNPs in TLR9 by ConSurf server

Looking into insight, an evolutionarily conserved amino acid profile of most deleterious nsSNPs provides essential information to understand the possible impact of deleterious nsSNPs on the structural and functional features of the protein. The ConSurf web server identifies whether the position of every individual amino acid is exposed on the protein surface or buried inside the core protein. The outcomes of the ConSurf web server revealed six deleterious nsSNPs (P139H, R257C, L283P, G514D, L544Q, and H566Y) as evolutionary conserved (Table 3 and Figure 2). However, the other two deleterious nsSNPs C265Y and W670R are predicted to be as moderately conserved (Table 3). The two variants (G514D and H566Y) are structurally exposed whereas the other six nsSNPS (P139H, R257C, L283P, L544Q, and W670R) are buried inside the core protein (Table 3). Evolutionarily, conserved amino acid residues are known to play a vital role in crucial structural and functional features of the protein. Therefore, it is argued that mutations at highly conserved sites have a direct impact on changes in the structural and functional features of TLR9.

Figure 2. Conservation profile analysis of most deleterious nsSNPs. Figure showing the amino acid conservation profile of human TLR9 protein. Highly conserved residues are directly proportional to the density of dark red color (score varies between 6 and 9), whereas the density of blue color (score varies between 6 and 9) shows variability in respective amino acids. Identified, eight most deleterious nsSNPs sites are highlighted by red arrows.

Table 3. Evolutionary conservation profile analysis of high-risk nsSNPs in TLR9.

Amino acid residues	Consurf			Netsurf
	Class assignment	Conservation score	Conserved/ Variable	Class assignment	RSA	ASA (Å^2)
P139H	Buried	8	Conserved	Buried	0.19	27.67
R257C	Buried	8	Conserved	Exposed	0.46	105.98
C265Y	Buried	7	Conserved	Buried	0.13	18.77
L283P	Buried	8	Conserved	Buried	0.063	11.55
G514D	Exposed	8	Conserved	Exposed	0.49	39.03
L544Q	Buried	8	Conserved	Buried	0.098	18.03
H566Y	Exposed	9	Conserved	Exposed	0.38	45.47
W670R	Buried	7	Conserved	Buried	0.171	41.02

RSA: relative surface accessibility = ASA / Max ASA.

ASA: absolute surface accessibility.

Conservation score: 1= most variable position and 9= most conserved position.

Analyzing the effect of most deleterious nsSNPs on protein stability and solvent accessibility by NetSurfP server

The solvent accessibility and protein stability are essential physiological parameters to understand the possible impact of most deleterious nsSNPs on the functional behaviors of TLR9 protein. The NetSurfP web server characterized the cumulative impact of every individual deleterious nsSNPs in terms of exposed or buried amino acid along with solvent accessibility parameters i.e., absolute surface accessibility (ASA) and relative surface accessibility (RSA). NetSurfP predicts that five nsSNPs (P139H, C265Y, L283P, L544Q, and W670R) are buried inside the backbone of core protein; however, three (R257C, G514D, and H566Y) are exposed over the protein surface (Table 3). The RSA score of most deleterious nsSNPs varies between 0.063 ≤ RSA ≥ 0.49 and ASA score varies between 11.55 ≤ ASA ≥ 105.98. The NetSurfP results showed higher value of RSA and ASA score for three SNPs and these are exposed on protein surface i.e., R257C (RSA = 0.46, ASA= 105.98), G514D (RSA= 0.49, ASA= 39.03), and H566Y (RSA= 0.38, ASA= 45.47) (Table 3). Of the total, five most deleterious nsSNPs are buried in core protein and could be responsible for protein folding and functional changes. However, the three most predicted deleterious nsSNPs exposed on the protein surface might alter the physiological parameters of TLR9 protein i.e., protein thermostability, surface charge potential, accessible surface area, and volume. Further, based on the altogether outcomes of Consurf and NetSurfP analysis, we could argue that the identified most deleterious eight nsSNPs might have a direct impact on altering the biological and functional features of TLR9 and that could be leading to the establishment of various types of cancer and other chronic diseases.

Structural assessment analysis of most deleterious nsSNPs on TLR9 by HOPE server

Project HOPE is a web-based application used to identify the effect of amino acid substitution on the native structure of TLR9 protein, including physiological and chemical properties i.e., hydrophobicity, spatial structure, charge, size changes between wild and mutant, and protein functions (Table 4). Accordingly, HOPE indicates, of the total three amino acid residues corresponding to the most deleterious nsSNPs (P139H, G514D, and L544Q) are bigger as compared to the wild-type. However, the two most deleterious nsSNPs L283P and R257C are comparatively smaller than the wild-type. Our previous Consurf analysis suggests that the mutation site G514D is exposed on the protein surface. Additionally, the inspection of P139H and L544Q on TLR9 showed mutated residues are bigger than native residues and might not perfectly fit in the core protein. Further, at site P139H, the mutated residue is positively charged whereas the native residue is neutral, which can lead to a problem in protein folding.

Table 4. Interpreting the impact of amino acid changes on TLR9 protein structure, hydrophobicity, charge, and size by HOPE.

Amino acid change	Change of structure	Change of hydrophobicity	Change of size	Change of charge	Interpretation
P139H					Might disturb the main activity of the protein. Might disturb the local structure. Disturb local structure.
R257C				Positive > Neutral	Loss of Hydrogen bonds. Might disturb correct folding.
C265Y				Neutral > Negative	Flexibility of protein may get abolished. It can cause repulsion of ligands. It might lead to bumps and disturb the torsion angles. cysteine bridge interactions might also be lost.
L283P					Might lead to loss of interactions.
G514D				Neutral > Negative	Might lead to bumps. The torsion angles for this residue are unusual. Disturb local structure.
L544Q					Might lead to bumps. Hydrophobic interactions either in the core or on the surface of the protein may be lost.
H566Y					May result in loss of hydrogen bonds. Might disturb the protein correct folding.
W670R				Positive > Neutral	Might lead to loss of interactions. Hydrophobic interactions, either in the core of the protein or on the surface, will be lost.

The HOPE server identifies the effect of point mutations based on the physiological changes in the nature of amino acids. The change in physical features of amino acid residues might affect their impact on the interacting surrounding of the protein, which could lead to altered protein conformation and protein folding.

Interestingly, of the total, six mutants (P139H, R257C, G514D, L544Q, H566Y, and W670R) showed an increased level of hydrophobicity in mutant-type as compared to the wild-type, whereas the reverse has been observed for mutant L283P (Table 4). In addition, the charge potential differences at key point mutations (R257C, C265Y, G514D, W670R, and D683K) are also observed between the wild-type and mutant-type. The inspection at these mutation sites showed that in R257C mutation residue C is neutral, whereas wild-type residue R has a positive charge. Similarly, in the charge potential inspection of other mutants i.e., C265Y, G514D, W670R, and D683K, the mutant residues have negative, negative, positive, and positive charges, respectively, whereas wild-type residues have neutral, neutral, neutral, and negative charges, respectively. Here, our comparative analysis suggests that the distribution of hydrophobicity differences between wild-type and mutant-type residues in most deleterious nsSNPs can alter the side-chain intermolecular hydrogen bonding network within two or more neighboring residues. These changes could be leading to disrupting the mechanism of protein folding and other crucial structural differences in TLR9 protein. In general, we could suggest that above discussed crucial changes might affect the catalytic behavior of TLR9 protein, which could be leading to altering the biological role of the protein.

Protein–protein interaction analysis of TLR9 by STRING database

Globally, the STRING database is a popular tool to understand protein–protein interactions. Any variation in protein structure and function can alter its ability to interact with other molecules. Thus, the interacting neighbors of mutated protein could be leading to phenotypic effects. The investigation of the STRING database revealed that TLR9 protein is related to 10 different proteins likewise MYD88, HMGB1, UNC93B1, TIRAP, TICAM1, PLSCR1, IL6, IL1B, TICAM2, and TLR2 (Figure 3). Our findings of the STRING database have been associated with some of the previously reported experimental studies where the interactions of TLR9 with other protein molecules are well known (Takeda and Akira 2015; Han et al. 2018; Kim et al. 2020; Dong et al. 2022).

Figure 3. Protein–protein interactions of TLR9 protein. The figure shows protein–protein interactions of TLR9 protein with its related proteins MYD88, HMGB1, UNC93B1, TIRAP, TICAM1, PLSCR1, IL6, IL1B, TICAM2, and TLR2. The interaction network has been generated by the STRING database. A brief introduction about all interacting partners of TLR9 has been also highlighted below in the network diagram.

Three-dimensional structural assessment analysis of TLR9 protein

Three-dimensional structural assessment analysis of human TLR9 protein is performed to understand the effect of most deleterious nsSNPs on structural features of wild-type versus mutant-type. The 3D coordinates of wild-type and mutant-type TLR9 are predicted using Modeller10.1 and TLR9 from Equus caballus (PDB ID: 3WPB) as templates (Ohto et al. 2015). Ramachandran plot analysis showed that for each structure > 99% of residues have been falling in the favorable region of the Ramachandran plot (data not shown).

Structural superimposition of wild-type and mutant-type suggested that the protein backbone of both structures is aligning perfectly with each other; however, deviations are noted among all points of mutation (Figure 4). Further, it has been suggested that mutation sites identified as buried inside the core protein, might alter intermolecular hydrogen bonding with their adjacent amino acids. This might be led to altering the protein folding mechanism of TLR9 and might alter interactions between TLR9 and ligand molecules leading to affecting its catalytic processes. Additionally, mutation sites identified at the protein surface of TLR9 affect the protein–protein interactions network with surrounding proteins leading to the phenotype effect. Therefore, our comparative structural assessment analysis of wild versus mutant TLR9 suggests that the positional deviations at mutation sites could lead to structural and functional changes among the mutant-type.

Figure 4. Structural alignment of wild and mutant-type TLR9 proteins: The C-alpha atoms structural superimposition of wild-type (sky blue color) and mutant-type (sand color) TLR9 proteins has been shown in cartoon format. All sites have been shown in green and red color sticks for wild and mutant-type, respectively. The figure shows that both structures are aligning perfectly where slight deviations are noted among all points of mutations. The overall structural alignment shows the very close C-alpha atom’s root mean square deviation (RMSD) of 0.064 Å between the two structures.

Amino acid network analysis for identifying the effect of mutation on noncovalent intra-molecular interactions of TLR9 protein using the RING 3.0 server

Comparative amino acid network (AAN) analyses for wild and mutant TLR9 proteins have been performed in order to identify the impact of point mutations on noncovalent intramolecular interactions (Figure 5). The outcomes of the RING 3.0 server suggest changes in the degree of nodes at the point of mutations (Figure 5). Accordingly, we observed the degree of nodes as 0, 4, 4, 3, 1, 6, and 6 for the amino acids P139, R257, C265, L283, G514, L544, H566, and W670 respectively of the wild TLR9 protein. Whereas the degree of nodes at the points of mutations P139H, R257C, C265Y, L283P, G514D, L544Q, H566Y, and W670R are observed as 5, 3, 7, 2, 4, 6, 7, 6, and 5 respectively for the mutant TLR9 protein (Figure 5). The altogether outcomes of AAN for wild and mutant TLR9 protein suggests the direct impact of point mutations on various types of noncovalent intramolecular interactions (Figure 5 and Table 5). Accordingly, we observed 458 and 514 intra-molecular H-bonds for wild and mutant TLR9 protein, respectively (Table 5). The higher density of π-π stacking, π-cation, and van der Waal interactions have been observed for mutant as compared to the wild TLR9 protein (Table 5). However, three ionic interactions are observed to be conserved in both cases, and eight and seven disulfide bridges have been observed for the wild and mutant TLR9 protein, respectively (Table 5). Further, the density of van der Waals forces for the wild and mutant TLR9 protein is observed as 774 and 812, respectively (Table 5). Interestingly, we observed the total of noncovalent intramolecular binding interactions as 1277 and 1372 for the wild and mutant TLR9 protein, respectively (Table 5). Since our results showed a higher density of noncovalent intramolecular binding interactions for the mutant TLR9 as compared to the wild TLR9, it is suggested that point mutations directly affect the noncovalent intramolecular interactions network of TLR9 protein. The point mutations in TLR9 protein enhance the overall noncovalent intramolecular binding interactions, which might stabilize the structural scaffold of mutant TLR9 protein as compared to the wild case. Further, based on the detailed finding of our analysis we can suggest that the variability in the degree of nodes results in the alteration of the protein conformation leading to changes in protein folding. In general, we can suggest structural changes in TLR9 proteins resulting in variabilities of biophysical and functional features led by point mutations (Yan et al. 2014).

Figure 5. Amino acid network (AAN) analysis of wild-type and mutant-type TLR9 proteins. The figure shows comparative AAN analyses for wild and mutant TLR9 proteins. The nodes of interest (point of mutations) have been labeled and highlighted by black circles. The size of solid spheres is directly proportional to the density of the degree of nodes. The various types of noncovalent bonded intramolecular contacts have been shown by lines. The colors of nodes are varying according to the nature of the amino acids (Supplementary Figure S1). (A) The figure shows AAN analysis for the wild-type TLR9 protein. For wild-case TLR9, only seven nodes have been highlighted because an amino acid P139 does not have any degree of nodes. (B) The figure shows AAN analysis for the mutant-type TLR9 protein.

Table 5. Detailed intramolecular noncovalent interactions analysis of wild and mutant TLR9 protein by using online available Ring server.

Types of interactions	Wild	Mutant
Intramolecular H-bond	458	514
π-π stacking forces	34	35
π-cation interactions	00	01
Ionic interactions	03	03
Disulfide bridges	08	07
Van der Waal forces	774	812
Total noncovalent interactions	1277	1372

Ring 3.0 server is used for the analysis of noncovalent intramolecular interactions in TLR9 wild and mutant protein. The physiological parameters such as intramolecular H-bond, π-π stacking forces, π-cation interactions, and van der Waal forces are found to be higher in the mutant protein, whereas disulfide bridges have been higher in the wild protein.

Association of the damaging nsSNPs with cancer

The cBioPortal and canSAR Black web server are used to identify the cancer-specific nsSNPs. The comparative outcomes of both tools identify two nsSNPs G514D and W670R that are associated with uterine corpus endometrial carcinoma (UCEC). However, none of the nsSNPs are found to be associated with cervical adenocarcinoma (CESC). Further, based on these findings we could suggest that tissue specificity of mutations in canSAR black data at two crucial positions 514 and 670 in TLR9 may result in moderate severity in UCEC.

Analyzing the gene expression profile for TLR9 gene in association with cancer

GEPIA2 database has been used to perform box plot analysis of TLR9 gene with respect to identifying the differential gene expression associated with cancer types. The box plot analysis revealed that CESC arises due to underexpression of TLR9, whereas UCEC results due to the overexpression of the TLR9 gene (Figure 6). Our results suggest that the majority of cancers occur due to the downregulation of TLR9; alternatively, few are resulted by the upregulation of TLR9. A similar finding has been already confirmed by some previous experimental studies where the downregulation of the TLR9 expression is associated with the establishment of oral cancer (Han et al. 2018). Further, our finding has been also strengthened by a study where a population of chronic HBV carriers showed suppressed expression and function of TLR9, assisting in immune escape strategy for HBV (Tout et al. 2018).

Figure 6. Gene expression analysis of TLR9. Figure showing detailed Box plot analysis diagram of the TLR9 gene expression in the case of cervical squamous cell carcinoma (CESC) and uterine corpus endometrial carcinoma (UCEC) for both tumor (red color box) and normal (dark grey color box) samples. The analysis showed that UCEC shows the overexpression of the TLR9 gene, whereas underexpression of the gene is observed in CESC.

Cumulative TLR9 gene expression and survival analysis in CESC and UCEC patients

The online available GEPIA database has been used to perform detailed survival analyses of patients affected by CESC and UCEC. Depending on the median TLR9 gene expression level, patients are categorized into two groups namely high or low expression level groups. Our results suggest that the abnormal expression of the TLR9 gene is involved in the poor prognosis of patients in association with CESC and UCEC. Additionally, the patients associated with CESC showed long but similar survival time (approximately 200 months) in both under and overexpression of the TLR9 gene (Figure 7(A)). However, patients associated with UCEC showed a shorter survival time (approximately 125 months) under expressed conditions of the TLR9 gene, whereas a comparatively longer survival time (approximately 140 months) has been observed in the overexpression of the TLR9 gene (Figure 7(B)). The TLR9 gene expression and survival analysis in CESC and UCEC patients suggest the clinical role of TLR9 in the association of various types of cancers.

Figure 7. Patient’s time-dependent percentile survival rate analysis. The figure shows the overall percent survival rate of the patient for both (A) cervical squamous cell carcinoma (CESC) and (B) uterine corpus endometrial carcinoma (UCEC). The blue color represents low expression and the red color represents the high expression of TLR9. The patients associated with CESC showed longer but similar survival time (approximately 200 months) in both under and overexpression of the TLR9 gene. The patients associated with UCEC showed a shorter survival time (approximately 125 months) during underexpressed conditions, whereas comparatively longer survival time (approximately 140 months) during the overexpression of the TLR9 gene.

Identifying the effect of most deleterious nsSNPs and expression of TLR9 gene with their neighboring genes

The muTarget tool identifies the effect of mutation on gene expression and gives an expression plot. Here, in the current hypothesis we have used two independent analyses i.e., ‘Genotype’ and ‘Target’ for the TLR9 gene to predict the detailed role of the mutant gene on its expression. In this respect, we performed a ‘Genotype’ run for the TLR9 mutant gene. For uterine cancer, the outcomes of the ‘Genotype’ hypothesis revealed the alteration in the expression of MKNK2, MED18, AEN, C19orf25, OSGIN1, and other genes (Figure 8). However, cervical cancer mutation in TLR9 could not alter the expression of neighboring genes. Next, we performed a ‘Target’ run by selecting TLR9 as the target gene for both cervical and uterine cancers. Our results showed gene alteration for cervical cancer (SOS1, NMT2, EPB41L3, and PZP) where a higher expression is observed for the MUC4 gene (Figure 9(A)). Similarly, for uterine cancer the expression of TLR9 can be influenced by a mutation in various genes i.e., TP53, AGO1, GRM8, SCAF4, and MROH8 (Figure 9(B)). Interestingly, we observed that gene TP53 showed a maximum change in gene expression for uterine cancer as compared to other counterparts’ genes. The comparative outcomes of our muTarget analysis suggest that the expression of various genes has been affected by point mutations in the TLR9 gene. In this respect, with regard to the case of uterine cancer, the expression of MKNK2, MED18, AEN, C19orf25, and OSGIN1 genes have been altered along with mutations in the TLR9 gene (detailed about potential genes are incorporated in Supplementary Data S2).

Figure 8. Genotype analysis of TLR9 gene. The figure represents a genotype run analysis to identify the changes in the TLR9 gene expression for uterine cancer. The figure shows that the expression of MKNK2, MED18, AEN, C19orf25, OSGIN1, etc. neighboring genes in uterine cancer is downregulated in presence of the mutated TLR9 gene (red) and normal expression in presence of wild-type (green) TLR9 gene expression.

Figure 9. Target run analysis of TLR9. Target run analysis is performed to identify mutations in crucial genes that can alter the expression of the TLR9 gene. (A) Highlighted mutations in significant genes belong to cervical cancer. The figure shows the downregulation of the TLR9 gene in presence of mutated SOS1, NMT2, EPB41L3, and PZP and MUC4 genes in cervical cancer. (B) Highlighted mutations in significant genes belong to uterine cancer. The figure shows the expression of the TLR9 gene is downregulated by a mutation in different genes i.e., AGO1, GRM8, SCAF4, and MROH8 except for mutated TP53 where the TLR9 expression is upregulated.

In the current study, we used TLR9 as the target gene where we found that its expression is changed by numerous mutated genes. Accordingly, the tool predicts the mutations in SOS1, NMT2, EPB41L3, and PZP genes, which might change the gene expression of TLR9 in cervical cancer patients. Further, in a uterine cancer patient, mutations in TP53, AGO1, GRM8, SCAF4, and MROH8 genes change the gene expression of TLR9 (details about potential genes are incorporated in Supplementary Data S2). Since all identified genes are known to play a crucial role in various types of metabolic pathways and biological processes, mutations in these key genes could enhance the possibility to cause various types of cancers. Accordingly, based on the comparative outcomes of muTarget analysis, it is suggested that mutations in any interacting gene partners in the gene-gene network also affect their neighboring genes, which could lead to cancer formation.

TLR9 are significant molecules that recognize and interact with DNA viruses. The presence of missense SNP in a gene responsible for amino acid substitutions causes functional changes or deleterious effects in the protein. The improper function of the TLR9 protein could further lead to inappropriate activation of the downstream immune pathway.

Here, we present the first systematic and exhaustive study using various in silico approaches leading to the identification of the most deleterious nsSNPs in the TLR9 gene. With the help of various computational tools and algorithms, we reported the eight most deleterious nsSNPs (P139H, R257C, C265Y, L283P, G514D, L544Q, H566Y, and W670R) in the TLR9 gene. The deleterious nsSNPs at highly conservation sites act as potentially damaging SNPs and could influence protein stability and its functional features. Further, our study identifies two crucial nsSNPs G514D and W670R, specifically associated with the severity of UCEC. Further, the key findings of the current study have been strengthened by some similar experimental studies (as discussed below) where the clinical roles of a few SNPs related to the TLR9 gene are highlighted.

Several previous studies suggest that the presence of SNPs led to hampering of the DNA recognition ability of the TLR9 receptor. In this respect, few of the reports highlighted the medical importance of TLR9 where its disease susceptibility and role are affected by the presence of SNPs. In a Swedish study, the clinical role of TLR9 SNPs (A1635G and G1174A) has been discussed where these SNPs are found to be associated with HIV-1 infection (Bochud et al. 2007). Another study showed that the TC genotype at locus -1237 T/C (rs5743836) of TLR9 in mononuclear cells gets upregulated expression and promiscuous response towards the CpG along with profound IL6 and B-cell proliferation (Carvalho et al. 2011). Similar SNP has been studied in Brazilian women to find the placental infection, where the presence of the C allele decreased the infection of human cytomegalovirus and another study showed a higher frequency of allele with allergic bronchopulmonary aspergillosis (Carvalho et al. 2008; Sánchez et al. 2021). A recent study on cervical cancer and cervicitis in India showed a decreased risk of HPV16 infection with three SNPs – 1486 T/C(rs187084), 2848 G/A (rs352140), and 1174 A/G (rs352139) of the TLR9 gene (Pandey et al. 2021). Similarly, a study on Han Chinese women with HPV16 infection carrying GA/AA genotype of G2848A (rs352140) showed an increased risk of cervical cancer as compared to women carrying GG genotype (Jin et al. 2017). A simultaneous contradictory meta-genomic study on rs352140 SNP of TLR9 gene with vulnerability to cervical carcinoma showed no association with the disease (Parsaeian et al. 2021).

Based on our outcomes and in association with the previous experimental reports highlighting the clinical role of SNPs in TLR9, we could conclude that this is the first systematic study to identify the possible deleterious role of nsSNPs in the TLR9 gene related to cancer. Therefore, the present study might be able to delineate the guidelines for identifying the most damaging SNPs and predicting the risk factors for cancer and disease susceptibilities. The current study has been focused on the most deleterious nsSNPs associated with endometrial and cervical cancer; however, the role of these SNPs could also be checked for other severe diseases. The 3D coordinates for human TLR9 proteins have not been solved by any experimental methods; therefore, the modeled structures have been used for detailed structural assessment analyses that might also slightly fluctuate the accuracy of our concluding remarks. Although, the findings of this study have been calculated by using various popular tools and algorithms, however, an extensive experimental and large-population-based study along with clinical trials are needed before implementing the key findings of this study for clinical purposes.

Materials and methods

The protein sequence of human TLR9 (Q9NR96) has been retrieved from the UniProtKB database and its related TLR9 gene has been retrieved from the NCBI database. The TLR9 gene-related 263 somatic mutations reported SNPs retrieved from the GDC portal of the TCGA database (https://www.cancer.gov/tcga). In the current study, we have used various bioinformatics tools and algorithms (SIFT, SNAP2.0, PROVEAN, and PhD-SNP) to identify the most deleterious nsSNPs. The results of these tools analyze the phenotypic effect of amino acid changes on protein functions as neutral/tolerated or deleterious/effect based on the score (Bromberg and Rost 2007; Sim et al. 2012) (Supplementary Data S1). Align GVGD works based on multiple sequence alignment and biophysical properties of protein residues where score C15 is less likely to affect and C65 is more likely to affect the protein function (Tavtigian et al. 2006). PolyPhen-2 identifies the effect of amino acid substitution on structural and functional features of protein as neutral, possibly damaging, or probably damaging (Ramensky et al. 2002). PMut tool is a neural network-based classifier that annotates pathological mutations on protein from 0–0.5 (neutral) to 0.5–1 (pathological) (Costa et al. 2005). Additionally, PROVEAN is also used to identify the effect of respective nsSNPs as deleterious or neutral. If the retrieved score is ≤ 2.5, then it is treated as a deleterious effect or otherwise as a neutral effect (Choi et al. 2002). Subsequently, PhD-SNP has been used to identify the role of respective nsSNPs as neutral or deleterious in homo sapiens related to disease formation (Capriotti et al. 2006). Next, DynaMut2, SDM, and I-Mutant Suite are employed to identify the impact of nsSNPs on protein stability based on a free energy change value (DDG) (Capriotti et al. 2005; Dehury et al. 2019, 2017, 2020; Pandurangan et al. 2017; Rodrigues et al. 2021). The Consurf server identifies the conservation profile of mutated residues at the nsSNP positions in protein molecules, which are represented in a color scheme (Berezin et al. 2004). Further, NetSurfP-2.0 predicts the secondary structure, solvent, and surface accessibility (Klausen et al. 2019). Further, HOPE is used to find respective functions such as changes in hydrophobicity, size, charge, etc. of mutation on the protein (Venselaar et al. 2010). Additionally, the STRING database has been used to predict related partners of the TLR9 gene (Jensen et al. 2009).

A suitable template of Homo Sapiens TLR9 for homology modeling has been searched by using PSI-blast against the PDB database. Further, the popular biomolecular modeling software Modeller10.1 has been used for homology modeling of wild and mutant proteins (Webb and Sali 2016). Qualitative analysis, refinement, and structural validation of each modeled structure have been done by using the SAVES server and ProSAweb-based tool (Wiederstein and Sippl 2007). Further, PyMOL molecular visualizations tool (pymol.org) has been used for molecular visualizations and structural superimposition (Schrödinger and DeLano 2020). Additionally, the RING (Residue Interaction Network Generator) 3.0 server has been used to identify noncovalent interactions at the atomic level in protein structures based on geometrical parameters (Clementel et al. 2022).

The present study is designed to identify and delineate the critical role of specific nsSNPs, which could be associated with endometrial and cervical cancer. In this respect, cBioPortal and canSAR Black tools have been used to access cancer genomics data and identify specific nsSNPs associated with cancer (Cerami et al. 2012). Further, GEPIA has been used to correlate gene expression and the overall survival rate (Tang et al. 2019). Next, we have used the muTarget tool for analyzing the correlation of gene expression (Genotype and Target run) and mutation based on the TCGA data (Nagy and Győrffy 2021) (detailed methodology has been given as in Supplementary Data S1).

Ethical approval

This article does not contain any studies with human participants or animals performed by any of the authors.

Authors’ contributions

Experiments, data curation, data analysis, literature review, and writing original draft: H.G.; data curation, data analysis, literature review, and writing original draft: V.V.V.; supervision, review, and editing: S.A.H.; conceptualization, supervision, investigation, project administration, funding acquisition, data interpretation, validation, review, and editing: M.B.

Abbreviations
CpG	=	cytosine-phosphate-guanine
TLRs	=	Toll-like receptors
nsSNPs	=	Nonsynonymous single nucleotide polymorphisms
SIFT	=	Sorting Intolerant from Tolerant
SNAP2.0	=	Screening for non-acceptable Polymorphism v2
PROVEAN	=	Protein variation effect analyzer
PhD-SNP	=	Predictor of human deleterious single nucleotide polymorphism
HOPE	=	Have (y) Our Protein Explained
STRING database	=	Search Tool for the Retrieval of Interacting Genes/Proteins
PDB	=	Protein Data Bank
TCGA	=	The Cancer Genome Atlas
GDC	=	Genomic Data Commons
RING	=	Residue Interaction Network Generator

Acknowledgments

Heena Gautam acknowledges the Council of Scientific and Industrial Research, India 09/1173(001)/2017-EMR/I for financial assistance as a Senior Research Fellowship. Ved Vrat Verma acknowledges the Indian Council of Medical Research, India, BMI/11(27)/2020 for financial assistance as a Research Associate Fellowship.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The details of third-party software and databases used in this study are provided in Supplementary Table S1.