Thyroid peroxidase (TPO) is essential for the synthesis of thyroid hormones. However, specific mutations render TPO antigenic and prone to autoimmune attacks leading to thyroid cancer, TPO deficiency, and congenital hypothyroidism (CH). Despite technological advancement, most experimental procedures cannot quickly identify the genetic causes of CH nor detect thyroid cancer in the early stages.
We performed saturated computational mutagenesis to calculate the folding energy changes (∆∆G) caused by missense mutations and analyzed the mutations involved in post-translational modifications (PTMs).
Our results showed that the functional important missense mutations occurred in the heme peroxidase domain. Through computational saturation mutagenesis, we identified the TPO mutations in G393 and G348 affecting protein stability and PTMs. Our folding energy calculations revealed that seven of nine somatic thyroid cancer mutations destabilized TPO.
These findings highlight the impact of these specific mutations on TPO stability, linking them to thyroid cancer and other genetic thyroid-related disorders. Our results show that computational mutagenesis of proteins provides a quick insight into rare mutations causing Mendelian disorders and cancers in humans.
Thyroid peroxidase (TPO) is an essential member of human peroxidase that acts as an oxidoreductase. TPO requires oxidization by hydrogen peroxide to carry out its function. The oxidized TPO plays a role in the production of thyroid hormones by catalyzing the oxidation of iodide into iodine. TPO further couples the iodine atoms with tyrosine residues (Habza-Kowalska et al. 2019). TPO forms an integral part of the plasma membrane through its intracellular, transmembrane, and extracellular domains. The extracellular domain further embeds the MPO-like, EGF-like and CCP-like domains. TPO also plays a role in embryonic hemopoiesis and response to oxidative stress (Mancini et al. 2016). Abnormal TPO activity or complete absence of TPO activity causes Hashimoto's disease, Grave's disease, total iodide organification defect (TIOD), congenital goiter, and congenital hypothyroidism 1 . Further, inhibition of TPO activity perturbs its antioxidative property, which drives inflammation. Several single nucleotide polymorphisms of TPO have been identified to be responsible for congenital defects.
In unrelated patients with TPO defects, the missense mutations N307T, V433M, P499L, and C808R, were presented as clinical evidence through Gel Electrophoresis of amplified TPO genes 2 . Another study concluded that heterozygous mutations would partially affect TPO function, while the homozygous inheritance of the mutation completely denatures the TPO protein 3 .
Kotani and others linked D240N to TIOD, causing a lack of thyroxine formation and TPO activity in the thyroid gland 4 . The localization of TPO offers it to serve as an antigen, which makes it a target for autoantibodies resulting in autoimmune thyroid diseases. Two missense mutations, R665W and G771R, caused a localization defect of TPO which led to congenital hypothyroidism 5 .
Congenital hypothyroidism (CH) is a universal neonatal disorder that can cause intellectual disability in newborns. CH is related to the absence of thyroid hormones (THs), categorized into permanent or transient cases according to clinical practice, and classified into primary, secondary, and tertiary cases according to the etiology 6 . Primary cases are caused by thyroid dysgenesis (TD), representing 85%, or dyshormonogenesis (DH), representing 10 to 15%.
Secondary cases are caused by pituitary injury, and tertiary cases are caused by hypothalamus malfunction 7 . The global CH incidence ranges from 1: 3000 to 1: 4000 in newborns 8 . However, the report from several newborn screening programs in the United States and worldwide shows an increasing rate of CH incidence of 1: 2000 [9][10][11] .
Somatic mutations of TPO have also been linked to various types of cancer, predominantly papillary carcinoma, and follicular carcinoma. These somatic mutations are often point mutations that cause single amino acid substitutions, which drive tumorigenesis in the affected parathyroid or thyroid gland 12 . The incidence of thyroid cancer in the US has increased alarmingly in recent years. In 2020, 53990 new cases of thyroid cancer were reported, attracting the focus of health professionals. However, the American Cancer Society estimated 43720 and 2120 thyroid cancer new cases and deaths for 2023, respectively 13 . To a certain extent, the increased cases of thyroid cancer can be attributed to increased detection, environmental factors, such as air pollution, radiation exposure and lifestyle choices. The Catalogue of Somatic Mutations in Cancer (COSMIC) confirmed nine non-redundant somatic mutations found in the tissues of the thyroid and parathyroid glands of patients with tumor 14 .
Post-translational modification (PTM) of TPO is necessary for its function. Interestingly, TPO is a highly glycosylated protein with four participating asparagine residues N129, N307, N342, and N569. Heme binding is also required for TPO's function at the cell surface. TPO has residues that form disulfide bonds required for its stability. However, with less confidence, PhosphoSitePlus predicted six phosphorylation sites, T52, Y55, T57, S97, T561, and T891, and one ubiquitination site, K879, on the TPO protein 15 . According to BioMuta, a database for single-nucleotide variations in cancer, mutations of specific PTM sites on TPO have been linked to lung (Y55*), stomach (Y55H), pancreatic (T57M), uterine (S97P), melanoma (T891P) and colorectal cancers (T561M)) 16 .
The advancement of sequencing technology has made a lot of gene and protein variants available. However, many TPO variants annotated in multiple databases are classified as having unknown significance due to small samples and the slow rate of experimental analysis. Due to the role variants play in the progression of chronic and complex diseases, annotating specific variants involved in a particular condition has become increasingly essential. The recent effort of computational biologists has sped up the annotation of functional variants, and the development of robust algorithms has accelerated the accurate classification and prediction of protein variants.
Computational mutagenesis can generate thousands of variants and predict their effects on protein stability, interactions, and functions. Understanding these effects gives scientists insights into the mechanisms behind diseases or allows them to develop new drugs that target specific proteins or pathways.
This study primarily uses folding energy calculation to analyze the effect of computed missense mutations on TPO stability. Then, it combines different computational tools to annotate functional variants to present helpful information to clinicians and scientists. The availability of these results would assist physicians to make informed decisions regarding the biological function of these variants and prevent the onset of chronic and congenital diseases associated with TPO.
We used an AlphaFold's monomeric structure (Identifier: AF-P07202-F1) predicted using deep learning algorithms for TPO structure analysis. AlphaFold is a highly accurate and reliable 3-D protein prediction tool, and it predicts protein structures with unknown similar structures.
Additionally, AlphaFold's structures have a high global distance test (GDT) with experimental structures 17 . We used the canonical sequence isoform (P07202-1) of TPO available on the UniProt database 18 .
We used in-house PERL programming scripts to create a full non-redundant missense mutation list from the TPO canonical sequence. We utilized Foldx 5.0 suite to calculate the folding energy change (∆∆G) caused by each missense mutation 19 . The Van der Waals of all atoms (vdw), solvation energy of apolar (solvH) and polar groups (solvP), intra-molecular hydrogen bond (hbond), electrostatic energy (el), backbone entropy cost (Mc), side chain entropy cost (Sc), electrostatic interactions (Kon), extra stabilizing free energy from water molecules (wb), and loss of translational and rotational entropy (Str) contribute to the folding energy (∆G) of the protein.
Protein engineering techniques were used to determine the individual energy terms, which are the default terms adopted by Foldx. There are two steps involved; First step is to repair the wildtype structure command using the 'RepairPDB' command. The repaired structure generated as output in the first step is now used as input in the second step. In the second step, we used the 'BuildModel' command to generate ∆∆G for each missense mutation.
The mathematical equation for calculating folding energy change (∆∆G) caused by missense mutation is shown below:
A positive ∆∆G value indicates an increase in folding energy leading to a less stable protein. A negative ∆∆G value will result in a more stable protein as less folding energy is required. When compared to experimental values, Foldx calculated ∆∆G values deviated by 0.46 kcal/mol. We, therefore, used ~0.5 kcal/mol as threshold to place the ∆∆G values into five categories: -highly stabilizing (ΔΔG ≤ 2.5 kcal/mol), stabilizing (-2.5 < ΔΔG < -0.5 kcal/mol), neutral (0.5 < ΔΔG < +0.5 kcal/mol), destabilizing (+0.5 < ΔΔG < 2.5 kcal/mol), and highly destabilizing (ΔΔG ≥ 2.5 kcal/mol).
We obtained 130 Pathogenic TPO missense mutations from the Human Gene Mutation Database (HGMD) 20 and cross-referenced from the Clinical Variation database (ClinVar) (https://www.ncbi.nlm.nih.gov/clinvar/) and the genome aggregation database (gnomAD) 21 . We filtered data from gnomAD and ClinVar based on pathogenic missense mutations from multiple submitters or submissions having at least one star.
We used the Meta-SNP and Polyphen2 tools to predict pathogenicity among the missense mutations. Meta-SNP uses the RandomForest algorithm and combines the prediction power of PhD-SNP, SNAP, SIFT, and PANTHER. Meta-SNP provides a reliability index (RI) within the range of 0 -10 for each prediction. The higher the RI the more accurate the pathogenic predictions 22 . Polyphen2 uses a machine learning algorithm to classify missense mutations as damaging to protein function or not 23 .
We generated 17727 possible TPO missense mutations and calculated their effects on protein stability using the Foldx software. Figure 1A shows the five categories of the folding energy change (∆∆G) caused by the missense mutations. Of the 17727 missense mutations, ~65% would increase the wild-type ∆G by at least 0.5 kcal/mol and ~9% would decrease the wild-type ∆G by a minimum of -0.5 kcal/mol. These mutations have the potential to disrupt the folding process of the TPO protein. As shown in Figure 1B, the median value of the effects of these missense mutations is ~1.22 kcal/mol. This is an indication that a novel missense mutation would most likely destabilize the TPO protein. The line chart below shows saturated and alanine scanning mutagenesis at each residue position.
On average, mutations of the residues located within the N-and C-Terminal would decrease the ∆∆G of TPO and hence increase its stability. However, missense mutations in the Haemperoxidase, sushi (740 -795), and EGF-like (796 -839) domains would most likely increase the ∆∆G of TPO. The residues, G393 and G348, have the highest destabilization potential, while S402 and G493 show the highest stabilization potential. The multiple sequence alignment (MSA) of TPO in humans, pigs, mice, and rats shows the conservation of the G393 residue.
However, there is a variation of the S402 residue in the TPO sequence of pig, mouse, and rat. The heatmap below shows the missense mutations, within the top ten residues, with the highest and lowest average ∆∆G values. G393W and G393Y significantly increased ∆G by 89.88 kcal/mol and 80.66 kcal/mol, respectively. Notably, S402M and S402P decreased ∆G by 4.65 kcal/mol and 4.29 kcal/mol, respectively. Any substitution at G393, G348, G587, G637, and G167 resulted in a significant increase in ∆G (>2.5 kcal/mol). For S402, G493, G516, G861, and G395, most changes in residue will result in a reasonable decrease in the ∆G (< -0.5 kcal/mol).
The local structure of the S402 and G393 shows their sidechains and their mutant amino acids, M402 and W393, respectively. Residues S402 and G393 are located within the Haem peroxidase domain and are 12 Angstrom apart on the TPO structure.
Several tools are now available to predict whether a missense mutation is disease-causing or neutral. In this study, we queried Meta-SNP and Polyphen2 tools to predict the pathogenicity of the saturated missense mutations on the TPO protein. Meta-SNP is a highly accurate metapredictor of single-point variations, and it harnesses the strengths of other predictors, such as PANTHER, SIFT, SNAP, and PHD-SNP. PolyPhen2 predicts the functional effect of nonsynonymous single nucleotide polymorphisms in humans. Table 1 shows the prediction of pathogenicity of Top 5 destabilizing and stabilizing mutations from the Foldx outputs. The top 5 destabilizing mutations from our mutagenesis analysis are all predicted by PolyPhen2 and Meta-SNP as disease-causing or pathogenic. However, only S402I and G493L out of the Top 5 stabilizing mutations are consensually predicted to be pathogenic by Polyphen2 and Meta-SNP.
Moreover, the Meta-SNP reliability index (RI) is higher with destabilizing mutations compared to the stabilizing mutations.
We generated boxplots to show the similarities between the results of saturated mutagenesis on TPO stability and the pathogenic predictions from Meta-SNP and PolyPhen2 tools. We also computed the statistical significance indicating the difference in the categories of the ∆∆G effects. Analysis Of Variance (ANOVA) of all five ∆∆G effect categories indicates a significant difference, with p-values less than 0.05 (P < 2e-16). The figure below shows that the neutral (-0.5<∆∆G<=0.5 kcal/mol) ∆∆G effect is less pathogenic than the other ∆∆G effect categories.
However, the highly destabilizing (∆∆G>2.5 kcal/mol) effect shows the highest correlation with pathogenicity predicted by Polyphen2 and Meta-SNP tools.
The 130 TPO missense/nonsense mutations from HGMD were truncated into 128 missense mutations, excluding one nonsense mutation, P204X, and one out-of-coverage missense mutation, N916S.
Our prediction shows that 80% of the HGMD missense mutations destabilize the TPO protein.
However, we also found that almost 16% of the HGMD missense mutations stabilize the TPO protein. As shown in the lollipop Plot below, S37P has been detected to cause Hypothyroidism in humans and we predicted that S37P stabilizes TPO by -1.35 kcal/mol. Also, G348R significantly increased the folding energy of TPO by 43.5 kcal/mol and has been found in patients screened with Hypothyroidism.
Whole-Exome studies showed the presence of a few somatic mutations in cancer of the parathyroid and thyroid tissues. Table 3 shows the effect caused by the nine thyroid and parathyroid cancer mutations on TPO stability. 7 out of the 9 reported thyroid and parathyroid cancer mutations destabilize TPO by increasing the folding energy by at least 0.5 kcal/mol.
Interestingly, T725P stabilized TPO and H676Q had a neutral effect on TPO stability. The parathyroid Adenoma has one mutation, W671R (2.91 kcal/mol), which is highly destabilizing.
Of the eight cancer mutations in the thyroid, P651S (3.49 kcal/mol) has the highest folding energy change while T725P (-0.83 kcal/mol) has the lowest folding energy change. Figure 6 shows the sidechains of the mutations, W671R and T725P. The residue W671 is located inside the TPO structure where it interacts with other neighboring atoms to stabilize TPO. However, T725 is situated on the periphery of the TPO structure where it binds to water molecules.
The EVE tool is a robust tool that trained its deep generative model using evolutionary sequences of related proteins 24 . EVE scores range from zero to one, scores closer to zero are benign and scores closer to one are pathogenic. However, EVE scores are generally classified as Benign, uncertain, or pathogenic. Classifying a certain percentage of variants as uncertain increases the accuracy of predicting pathogenic of benign variants. Overall, EVE is a better predicting tool for interpreting the clinical significance of variants.
We analyzed curated Post-Translational Modification (PTM) sites of TPO from the GlyGen and PhosphoSitePlus (PSP) databases. The GlyGen database integrated multiple datasets to create a model used to study the functional properties of glycans, especially in N-linked glycosylation 25 .
The PTM sites from PSP were identified using Proteomics experiments 26 . In addition, the PTM annotations for TPO were also available in the UniProtKB/Swiss-Prot database. The PSP phosphorylation and Ubiquitination sites were from at least two High Throughput Papers (HTP).
Figure 6 shows that mutations in the N-linked glycosylation sites stabilize TPO, with the exceptions of N129P and N342P which remarkably destabilize TPO by 3.62 kcal/mol and 8.39 kcal/mol, respectively. Interestingly, most mutations at the phosphorylation sites will likely destabilize TPO except for residue T891 with a mean ∆∆G value of -0.27 kcal/mol. K879 is the only ubiquitination site on the TPO protein. Mutations on K879 also decrease the ∆G value except for K879P which increases the ∆G value by 3.57 kcal/mol.
The local environment of a modification site contributes to its morphology and biological function. Missense mutations of flanking residues in the local environment of PTM sites can disrupt protein functions, thereby driving disease mechanisms 27 . The flanking residues typically constitute the 7 residues downstream (towards the C-terminal) and 7 residues upstream (towards the N-terminal) of the modification site. We mapped disease-associated missense mutations to PTM sites on the TPO protein sequence and found that the PTM-related deleterious mutations are associated with different phenotypes such as Hypothyroidism, Thyroid Peroxidase Deficiency and Partial Iodide Organification Effect (Table 4).
Computational analysis of missense mutations provides a quick and reliable preliminary assessment of disease phenotypes in humans 28 . In our previous study, we highlighted the reliability of our approach to identifying missense mutations in the Myeloperoxidase (MPO) protein 29 . In the case of TPO, chemiluminescence immunoassay and other experimental assays were used to identify some thyroid-related diseases. However, the experimental assay requires many samples, and it is tedious, slow, and labor-intensive. Comparatively, the computational approach is fast, accurate, and as reliable for the experimental design and data interpretation 30 .
Computational mutagenesis can identify rare variants in individuals with underlining thyroid diseases and intervene using biomarkers to prevent thyroid disease or avert worse prognosis 31 .
In this study, we highlighted the top ten missense mutations predicted to alter the folding energy of the TPO protein. G393W is the most destabilizing missense mutation, and it is in the hydrophobic core of the TPO protein structure. The G393 residue is located within the Haem peroxidase domain and was shown to be evolutionarily conserved. The conservation of residues within families of protein indicates their vital roles 32 .
Existing pathogenicity prediction tools such as PolyPhen2 and Meta-SNP provided more functional insight. We showed that missense mutations that destabilize TPO by >2.5 kcal/mol are most likely disease-causing. However, some missense mutations stabilize TPO and are pathogenic. S37P stabilized TPO structure and was detected in patients with congenital hypothyroidism. The study discovered that a biallelic TPO mutation of S37P resulted in reduced TPO activity 33 . Further, the S37 is a conserved residue among humans, pigs, mice, and rats. In some missense mutations: A443V (∆∆G= 2.94 kcal/mol), R769W (∆∆G= 1.52 kcal/mol), N592S
(∆∆G = 3.36 kcal/mol), and N798K (∆∆G = -0.80 kcal/mol), TPO had a reduced peroxidase activity, but its expression level was unaffected 34 . Analysis of the somatic mutations found in the thyroid and parathyroid tumor samples shows that highly destabilizing mutations can drive tumorigenesis. Studies by the Sanger Institute Cancer Genome Project highlighted nine tumorigenic mutations on the thyroid gland and surrounding tissue. In our results, A373S and T725P involved in thyroid neoplasm are also implicated in thyroid dyshormonogenesis and hypothyroidism, respectively.
Missense mutations of evolutionarily conserved residues drive pathogenicity. Bioinformatic analysis of top target mutations using the EVE tool showed that mutated conserved residues are pathogenic. The EVE tool predicted T244E as pathogenic, despite stabilizing the TPO protein.
PTM sites of TPO play a crucial role in its overall function. For instance, Deglycosylation inhibits TPO enzymatic activity by disrupting its tertiary structure 35 . We analyzed the mutagenesis of the N-linked glycosylation, phosphorylation, and ubiquitination sites on TPO.
Our results showed that S97, T561, N342, and N129 are the most critical residues with the least tolerance to missense mutations. Further, we showed that missense mutations on residues within the local environment of PTM sites can result in mendelian human diseases.
In summary, we performed a saturated computational mutagenesis to understand the biological and functional impact of missense mutations on the TPO protein. We found that most of the disease-causing missense mutations can destabilize TPO. Further, the effect of missense mutations also relies on the evolutionary conservation of the residue and its closeness to PTM sites. Thyroid cancers can be detected early through blood tests or fine needle aspiration biopsy (FNAB) of suspicious regions. Thyroid cancer can also be detected through molecular screening of related mutations. Understanding genetic mutations driving abnormal activity of TPO can be a reliable genetic screening tool against thyroid cancer.
Table 1. Pathogenic prediction of top destabilizing and stabilizing missense mutations.