Pathogenicity Prediction of Potential Variants in TULP1 Gene causing Hereditary RP: An In-silico Approach

Introduction

Eye is the main organ associated with vision and is a major part of the sensory nervous system. Poor eye vision or blindness can directly affect the quality of life. Several retinal damaging conditions including glaucoma, retinal vascular disease, and macular degeneration associated with aging, are known to affect approximately 60 million people worldwide. Some eye disorders occur due to certain environmental factors, for instance aging and poor health conditions. While, some are heredity traits where RP is among the most prominent. Approximately, one hundred (100) genes have been observed to cause hereditary RP. RP is a vast terminology used to group eye disorders with retinal damages leading to poor vision and in rare cases complete night blindness [1]. The prevalence ratio of RP is 1/3500 around the world [2].

There are two types of RP stated as non-syndromic RP and syndromic RP. Non-syndromic RP is related to poor vision and has been found to affect 70-80% of the people and is a more common type of RP. The incidence rate for non-syndromic RP is 1/5000, while carrier ratio observed is 1/1000. The other type of RP is syndromic RP where the retinal damage combines with some other systemic dysfunction. Usher syndrome is the typical example of syndromic RP in which vision loss conjugates with hearing disabilities [3].

The risk factor to transfer the hereditary disorders including hereditary RP as well as spastic paraplegia is high in South Asian region. As reported by statistical data, the autosomal recessive RP with the number of 60 genes is known to be involved [4].

TULP1 gene is an important member of TULP (Tubby-like protein) family and is found in retinal epithelial cells. TULP1, as a transmembrane protein, mainly functions in protein trafficking and has the tendency to bind with lipids forming phosphatidylserine, phosphatidylinositol, 3,4,5-bisphosphate, phosphatidic acid, 4,5-bisphosphate, 3,4-bisphosphate, phosphatidylinositol 5-phosphate, phosphatidylinositol 4-phosphate, and phosphatidylinositol 3-phosphate. TULP1 also functions in the cleaning of dead cells inside the retina and is involved in cell apoptosis/phagocytosis [5].

TULP1 gene with genomic size 14997 base pairs (15 exome count) is located on human chromosome 6, on the antisense strand of DNA (Deoxyribonucleic Acid) double helix. TULP1 gene encodes protein with 542 amino acids (Figure 1).

Figure 1. Pathogenic variants observed from cDNA (SNV) to protein. (An illustration of TULP1 gene (containing 15 exons, locates at 6p21.32) to protein translation (encoding 542 amino acid residues, region ranges between 298-542 covers a conserved C-Terminal domain and region between 1-289 linked with intrinsic disorder))

In the analysis of the current study, results were crossmatched using different tools and databases with different algorithms. The importance of algorithm designing was reflected upon for any biologist or healthcare expert. Not all the tools or software are 100% error free, however, still can be a helping hand for the scientists to make authentic predictions. Presently, scientists prefer dry lab practices before proceeding to wet lab experiments. This is because it is a more helpful approach in a limited time span for the scientists to predict the structure of any biological molecule even if it has not been discovered yet via using these tools or databases [6].

2. MATERIALS AND METHODS

TULP1 genomic variants were retrieved through gnomAD. Missense variants, splice site variants, frameshift, and stop-gained variants were analyzed for their pathogenicity grading using multiple tools. Protein functioning was observed through visualization tool. Streamline followed for the insilico approach is shown in Figure 2.

Figure 2. Flowchart to carry out TULP1 Gene Variants Analysis

2.1. Variants Retrieval

To study the effect of SNPs reported for TULP1 gene, these were retrieved from gnomAD gnomAD - Broad Institute, containing millions of SNVs (Single Nucleotide Variants) data available for human genome in both genome assemblies (GRCh37 and GRCh38) [7]. The number of variants obtained for TULP1 is shown in Figure 3.

Figure 3. Distribution of SNPs in TULP1 Gene

2.2. Analyzing Missense Variants

TULP1 missense variants were scored with a very low allelic frequency. Combined Annotation Dependent Depletion (CADD) tool was designed on machine learning (ML) algorithm combining more than 63 annotations. Variants with CADD C-Score range between 10-20 are considered highly pathogenic [8]. Variants with CADD C-score ≥ 10 were further analyzed using multiple tools, for instance CAPICE [9], REVEL, PMut, UMD-Predictor, and Predict SNP2.0 Mutation Assessor. Quality control of missense variants was verified using different parameters which included checking the protein stability using stability tools, (CUPSAT, iStable [10], DynaMut [11], MAESTRO, YASARA [12]) and Clashes finding inside the protein structure as well as PTM sites analysis through UCSF Chimera [13]. To observe the candidate gene functioning caused by variations, conservation ratio and secondary structures were also observed via ConSurf and NetSurfP-2.0 [14].

2.3. Splice Site, Frameshift, and Nonsense Variants Analysis

Scoring of TULP1 splice site variants was performed to screen out the highly pathogenic variants for disrupting the splicing mechanism. SPiCEv2.15 [15], Spliceman, Mutation T@ster [16] tools were used for variant analysis.

Frameshift and nonsense variants scoring was performed using CADD.

3. RESULTS

3.1. Analyzing Missense Variations

Through gnomAD, total number of 324 missense variants for candidate gene (TULP1) (allelic frequency <0.001 and allelic count < 50) were obtained. CADD scored 274 variants as pathogenic variants with C score ≥ 10. Most of the tools set a threshold value of 0.5, while variants score above 0.5 are considered pathogenic. To filter out the highly pathogenic variants, multiple tools were used, that is, CAPICE [9], REVEL, PMut, UMD-Predictor [17], Predict-SNP2.0, and Mutation Assessor. Lastly, only 20 missense variants with CAPICE score ≥ 0.5 were filtered out as "highly pathogenic" (Table 1). These twenty missense variants were further analyzed for candidate protein functioning. Graphical representation of CADD and CAPICE scoring of missense variants can be seen in Figure 4 (a).

Figure 4. Graphical Representation of pathogenic TULP1 genomic variants scoring (a) highly pathogenic missense variants using (CADD & CAPICE scoring). (b) protein Stability tools scoring vs missense variants (c) Splice site variants through SPiCE (Considering the Red line as threshold and dots above the threshold are variants affecting splicing while dots below threshold are less deleterious variants for splicing disturbance) (d)  Pathogenic frameshift variants scoring by CADD  (e) Stopgained variants by CADD

Table 1. TULP1 Pathogenic Missense Pathogenic Variants Scoring Obtained using Different Tools

(UMD predictor score predictions were divided by 100 to homogenize the data.)

Mutations with Asterisk (*) sign have already been reported in ClinVar and the rest are novel.

3.3.1. Effect of Missense Variations on Protein Stability. Highly pathogenic missense variants were analyzed by CUPSAT, i-Stable, DynaMut, MAESTRO, and YASARA scoring methods for stability check against each variant. Negative scores indicate the destabilizing effect of mutations on protein structure and vice versa for positive scores. In YASARA, variants with positive scores show decreasing protein energy due to the variant (Table 2).

Graphical representation of stability scores by tools used can be seen in Figure 4(b).

3D structure was modeled for TULP1 protein with missense variations through UCSF Chimera (Figure 5).

Table 2. Stability Score using Stability Tools to Analyze the Effect of Missense Mutations on Protein

Figure 5. TULP1 3D Structure Visualization by Angle Rotation of 180 Degrees (red area marked the mutated sites)

3.3.2. Clashes Finding due to Missense Variations. Clashes can be simply termed as unnecessary/unwanted contacts due to mutations,thereby changing the confirmation. Clashes were observed using UCSF Chimera.

Clashes were determined in mutations reported; p.Gly368Arg in which no clashes/contacts were found in wild type structure, while for mutated type 46, contacts were found.

p.Glu511Gly: no clashes were found for both wild type and mutant type, p.Tyr501Phe: no clashes were found for wild type. Whereas, 9 contacts including interaction with water molecule were observed for mutant type (Yellow color indicating the wild type of residue and purple color showing the mutant residue whereas red lines showing the contacts found, blue colored balls showing the water molecules) (Figure 6).

Figure 6. Clashes Findings due to Missense Mutations in Candidate Protein Structure using UCSF Chimera

3.3.3. PTM Site Disruption by Missense Variants. PTM sites inside candidate protein were confirmed using PROSITE database.

Phosphorylation sites with different protein kinase were observed and visualized in protein structure using UCSF Chimera. Some of the missense variations were found to disturb some phosphorylation sites, hence hindering the cellular processes (Figure 7).

Figure 7. PTM sites Disruptions Visualization through UCSF Chimera

3.3.4. Protein Conservation and Secondary Structure Prediction. Candidate protein (TULP1) has a single C-Terminal domain. To observe the conservation, ConSurf was used for analysis. ConSurf prediction scoring is done in terms of color grading (1-9) in which 1-4 grades are variables, grade 5 as average, and 6-9 color grade indicates highly conserved amino acid residues inside protein sequence.

NetSurfP-2.0 predicts surface accessibility as RSA (Relative Surface Accessibility), ASA (Absolute Surface Accessibility) and also provides information for secondary structures (coil, helix, strand) present in a protein. For secondary structures (loops or coils) inside TULP1 protein and to observe that amino acids are buried/exposed to the surface in sequence, NetSurfP-2 was used. Table 3, visuals seen in Figure 8.

Table 3. TULP1 Protein Conservation Prediction due to Missense Mutations using ConSurf and Solvent Accessibility, Secondary Structure Prediction by NetSurfP-2.0

Figure 8. TULP1 Protein Secondary Structure Prediction by NetSurfP-2.0. Surface Accessibility=Red high ground indicating exposed amino acids while sky blue downward high areas showing the buried residue in protein, considering threshold as 25%. Secondary Structure= straight pink line depicting coil, orange spirals are helix, colored arrow (indigo) indicating strand. Disorder= just under the secondary structure prediction line is thick black line showing the chances of disordered residues, the thicker the line, the higher the risk of disordered residues there. (as shown in figure 1)

3.2. Splice Site Variants Analysis

By gnomAD, 57 splice site variants were obtained for candidate gene (TULP1). Seven (7) variants were canonical splice site variants (splice site donor=4, splice site acceptor=3). The remaining 50 variants were scored using CADD. Thirteen (13) variants with CADD score ≥ 15 were scored by splicing tools (SPiCE v2.1.5, Spliceman, Mutation T@ster) (Table 4 & Figure 4 (c)).

Table 4. TULP1 Pathogenic Splice Site Variants Score

Asterisk (*) marked protein substitutions have already been reported in ClinVar, while the rest of mutations are novel.

3.3. Frameshift and Nonsense Variants Scoring

Nine (9) frameshift variants and seven (7) nonsense variants with the Allelic frequency filter (≤0.002) were retrieved for TULP1 and scored using CADD. Nine (9) frameshift variations were filtered as highly pathogenic (Table 5). While, stop-gained variations were scored high as compared to frameshift variants. (Table 6) Figure 4(d) and Figure 4(e) show the CADD scoring for frameshift and nonsense variants, respectively.

Table 5. Pathogenic Frameshift Variants in TULP1 Gene

Table 6. Pathogenic Stop-gained Variants in TULP1 Gene

Asterisk (*) marked protein substitutions have already been reported in ClinVar, while the rest of mutations are novel.

4. DISCUSSION

TULP1 gene variants are very rare and could be responsible for causing autosomal recessive RP. Protein native structure is important for proper functioning. If any sudden change or mutation occurs inside the amino acid sequence of protein, it directly affects its function. The wild structure type is more stable and mutation inside a protein structure can destabilize it. Stability is important for proteins" proper functioning [18].

TULP1 gene variants (missense, frameshift, nonsense, splice site) were scored according to their degree of pathogenicity using CADD. CADD scores the variants on three different ranks, that is, highly pathogenic, moderate, and lower pathogenic, referring to the main score PHRED (C Score). Variants with C score range between 10-20 are considered highly pathogenic [8]. CAPICE (Consequence Agnostic Pathogenicity Interpretation of Clinical Exome) based on ML algorithm outperformed compared to other tools, due to a low false positive rate [9].

After selecting highly pathogenic missense TULP1 variants, the amino acid substitutions were tested for stability. DynaMut predicts scores in two-step algorithms. The first step is utilizing "Holdout Random Sampler" which calculates the energy change represented as Cumulative distribution function. While in the second step, there comes the authentic stability score using cumulative distribution function which trains the final algorithm-"Neural network model" [11]. istable works on "Support Vector machine algorithm". Due to its uniquely designed algorithm, this tool proved a better option as compared to other tools" output for stability check. Both structural or sequential information of protein can be used as input methods [10] (Figure 4 (b)).

Protein with its native folds is crucial for its proper functioning. Mutations change its structure and, hence function. TULP1 encodes 542 amino acids contained protein structure which was visualized using Chimera. 3D structure with pathogenic missense mutations was modeled [13] (figure 5). Variants predicted as pathogenic for TULP1 were also visualized, disturbing PTM sites by 3D modeling (Figure 7).

TULP1 protein is a transmembrane protein involved in cellular trafficking. Therefore, its nearby interactions, cellular contents, and surface accessibility were observed using NetSurfP-2.0 [14]. In ConSurf, MSA is performed followed by constructing a phylogenetic tree and it also shows the conserved regions inside protein structure with colors [19] (Figure 8).

To find out about variants which can directly affect the splicing mechanisms, splicing variants retrieved through gnomAD were analyzed by Spliceman, SPiCE, and Mutation T@ster.

Spliceman based on "distance matrix calculation methods" provides the output as L1 distance. If position vectors are close enough, it denotes high ranks (high L1 distance). The higher the rank, the higher the chances of disturbing splicing sites [20]. SPiCE (Splicing Prediction in the Consensus Elements) tool designed on ML classifier "Logistic Regression" predicting two scores MaxEntScan(MES) and Splice-Site Finder-Like(SSF-Like) also provides output in a graphical form. In graphical representation (Figure 4c), variants above threshold line are considered as highly pathogenic. While, dots denoting the variants below the threshold line are considered neutral. Geneticists can easily apply this in diagnostics for the analysis of NGS results [20].

Probably, more than 100 genes have been found to be involved in causing RP. In the sub-continent, the custom of consanguineous marriages raises the risk of hereditary autosomal recessive RP. ABCA4, RHO, SPATA7, CRB1, and TULP1 gene families are important in Pakistan and Indian populations. This is because they can be the carrier to cause autosomal recessive RP [21].

Some novel genes, for instance ENSA, DDR1, CCDC188, and DACT2 were found to have an important role in vision cycle and ion channeling in photoreceptor cells. In 2020, a study was conducted to highlight the correlation of these four genes (ENSA, DDR1, CCDC188, DACT2) and the study concluded that the two genes (ENSA, DDR1) co-express and correlate in causing RP [22].

Recently, research was conducted to determine the pathogenic variants in TULP1. The results concluded two heterozygous mutations p. Arg482Trp and p. Leu504fs* as highly pathogenic and can cause autosomal recessive RP [23].

Finding the cure for hereditary disorders has always been a challenging task for the scientists. However, the options which have captured the attention of scientists include, cell therapy, retinal prostheses, and gene therapy. Gene therapy involves removing the genetic errors using viral or non-viral vectors followed by gene silencing pathways. Gene therapy implemented using phase3 voretigene neparvovec vector showed a better efficacy rate against inherited RPE65 mediated retinal dystrophy [24].

In 2020, a finding based on statistical data revealed that probably, 2.7 billion people around the globe (about 36% population) may be the healthy carriers for having at least one of the gene variants which may cause RP in future generations [25]. If scientists consider the cases of visual impairment/inherited blindness for detailed experiments, then this could be the road map in discovering new loci and genes causing RP [26].

4.1. Conclusion

RP is generically referred to as a group of heterogeneous eye disorders, where its complicated genetic patterns make it a complex class to study. TULP1 gene is expressed in rod and cone cells of the retina and contributes to 1% autosomal recessive RP. TULP1 gene variants were analyzed through an insilico approach to study the molecular basis of RP and highly pathogenic TULP1 variants were predicted. Predicted pathogenic variants for TULP1 can provide new insights for the scientists to design a treatment line for RP if caused by TULP1 gene variants. Candidate gene variants were cross-checked with different pathogenicity, predicting algorithms using different tools or databases to make predicted pathogenic scoring more authentic and reliable. Predicted pathogenic variants can be confirmed following wet lab experiments from blood sampling to sequencing.

Conflict of Interest

The author of the manuscript has no financial or non-financial conflict of interest in the subject matter or materials discussed in this manuscript.

Data Availability Statement

The data associated with this study will be provided by the corresponding author upon request.

Funding details

This research did not receive grant from any funding source or agency.