Mutation Spectrum in TPO Gene of Bangladeshi Patients with Thyroid Dyshormonogenesis and Analysis of the Effects of Different Mutations on the Structural Features and Functions of TPO Protein through In Silico Approach

Although thyroid dyshormonogenesis (TDH) accounts for 10-20% of congenital hypothyroidism (CH), the molecular etiology of TDH is unknown in Bangladesh. Thyroid peroxidase (TPO) is most frequently associated with TDH and the present study investigated the spectrum of TPO mutations in Bangladeshi patients and analyzed the effects of mutations on TPO protein structure through in silico approach. Sequencing-based analysis of TPO gene revealed four mutations in 36 diagnosed patients with TDH including three nonsynonymous mutations, namely, p.Ala373Ser, p.Ser398Thr, and p.Thr725Pro, and one synonymous mutation p.Pro715Pro. Homology modelling-based analysis of predicted structures of MPO-like domain (TPO142-738) and the full-length TPO protein (TPO1-933) revealed differences between mutant and wild type structures. Molecular docking studies were performed between predicted structures and heme. TPO1-933 predicted structure showed more reliable results in terms of interactions with the heme prosthetic group as the binding energies were -11.5 kcal/mol, -3.2 kcal/mol, -11.5 kcal/mol, and -7.9 kcal/mol for WT, p.Ala373Ser, p.Ser398Thr, and p.Thr725Pro, respectively, implying that p.Ala373Ser and p.Thr725Pro mutations were more damaging than p.Ser398Thr. However, for the TPO142-738 predicted structures, the binding energies were -11.9 kcal/mol, -10.8 kcal/mol, -2.5 kcal/mol, and -5.3 kcal/mol for the wild type protein, mutant proteins with p.Ala373Ser, p.Ser398Thr, and p.Thr725Pro substitutions, respectively. However, when the interactions between the crucial residues including residues His239, Arg396, Glu399, and His494 of TPO protein and heme were taken into consideration using both TPO1-933 and TPO142-738 predicted structures, it appeared that p.Ala373Ser and p.Thr725Pro could affect the interactions more severely than the p.Ser398Thr. Validation of the molecular docking results was performed by computer simulation in terms of quantum mechanics/molecular mechanics (QM/MM) and molecular dynamics (MD) simulation. In conclusion, the substitutions mutations, namely, p.Ala373Ser, p.Ser398Thr, and p.Thr725Pro, had been involved in Bangladeshi patients with TDH and molecular docking-based study revealed that these mutations had damaging effect on the TPO protein activity.


Introduction
Congenital hypothyroidism (CH) is defined as insufficient production of thyroid hormone at birth [1]. The global frequency of CH is 1 in 3000-4000 whereas a pilot study in Bangladesh reported an incidence rate of 1 in 1300 [2,3]. Thyroid dyshormonogenesis (TDH) which results due to defect in the pathway of thyroid hormone biosynthesis accounts for 10 to 20% cases of CH [4]. Till now, mutations in seven genes, namely, NIS (sodium iodine symporter, SLC5A5), PDS (Pendrin or SLC26A4), TPO (thyroid peroxidase), TG (thyroglobulin), IYD (iodotyrosine deiodinase, DEHAL1), DUOX2, and DUOXA2 (dual oxidase), have been reported to be involved in pathogenesis of TDH [5]. However, mutations in TPO gene which may result in total iodide organification defect (TIOD) or partial iodide organification defect (PIOD) have been described as the most common forms of TDH, and to date more than 60 different mutations have been described [6,7].
TPO is the major enzyme in thyroid hormone biosynthesis and it catalyzes both iodination and coupling of iodotyrosine residues in TG (thyroglobulin). It is a glycosylated hemeprotein of 110 kDa bound at the apical membrane of thyrocytes [8]. The single gene encoding TPO is located on chromosome 2p25 containing 17 exons and it spans at least 150 kb [9]. It is a homodimer protein (each monomer consists of 933 amino acid residues) and contains a peroxidase domain, three additional extracellular domains, a transmembrane helix, and a short C-terminal intracellular tail [10]. Although a low-resolution crystal structure of TPO has been reported, its high resolution structure remains to be determined [11,12]. The closest known homologues of TPO are lactoperoxidase (LPO), myeloperoxidase (MPO), and eosinophil peroxidase (EPO) with a sequence identity of 48%, 47%, and 47%, respectively. X-ray crystallographic structures have been determined for MPO and LPO [13,14], providing a platform for investigating the structural basis of TPO. Nevertheless, Sarah et al. investigated plausible modes of TPO structure and dimer organization through in silico approach [15]. However, how certain mutations affect the TPO protein structure and its enzyme activity is still unknown.
Investigation of mutational spectrum is paramount for genetic disorder study as it gives insight into the disease pathogenicity as well as severity. Screening and identification of mutational spectrum in the TPO gene of patients with TIOD and PIOD have been reported in different countries of the world like Argentina, Netherlands, Japan, Portugal, and China [6,[16][17][18][19]. In addition, changes in the TPO enzyme were assayed in vitro to compare mutant and wild type activities by several study groups and demonstrated mild to severe TPO enzyme inactivity for some mutations [20,21]. Though the prevalence of CH is more than twice the global incidence rate in Bangladesh, molecular basis of CH is still unknown in this country. Moreover, the effects of different mutations in the TPO gene on enzyme activity have not been investigated. In this study, we investigated the spectrum of mutations in TPO gene of patients with TDH and explored the possible effect of these mutations on the structure of TPO protein and TPO's MPO-like domain through in silico approach such as homology modelling, molecular docking followed by qunatum mechanics/molecular mechanics (QM/MM) and molecular dynamics (MD) simulation.

Study Participants.
A total of 36 confirmed cases of congenital hypothyroid patients with dyshormonogenesis were enrolled at the clinical care settings of National Institute of Nuclear Medicine and Allied Sciences (NINMAS) and department of Endocrinology, Bangabandhu Sheikh Mujib Medical University (BSMMU), Dhaka, for their follow-up examinations. Written informed consent was signed by the parents or guardians of the patients. 3 mL of blood was collected in an EDTA tube from each patient. Ethical clearance was obtained from the ethical committee of BSMMU and University of Dhaka.

DNA Isolation and Polymerase Chain Reaction (PCR)
Amplification. Genomic DNA was isolated from the EDTA blood by using a Qiagen DNAeasy mini kit. The isolated DNA was then amplified by PCR using TPO gene specific primers that together covered from Exon 8 to Exon 14, since global data showed that most of the common mutations in the TPO gene of the patients with congenital hypothyroidism were confined in this region. The primer sequences are listed in Table 1.
To amplify the desired target sequence of TPO gene, PCR amplification was conducted on a thermal cycler (Bio-Rad, USA). The final reaction volume was 10 l for each of  [25][26][27]. To validate the structure, the PDB format of the predicted 3D structure was submitted to both the servers. The result included the percentage of the amino acids having the average 3D-1D score of ≥ 0.2 for Verify3D server. The RAMPAGE webserver performs Ramachandran Plot Analysis by providing results including the percentages of the amino acid residues within the favored, allowed, and outlier regions. The results of the Verify3D and the RAMPAGE webservers were summarized in the results section.

Optimization of Heme Using Quantum Mechanical (QM) Calculations.
Since TPO is a heme containing protein, we investigated the effect of mutation on binding affinity of heme with TPO and its interactions with specific amino acid residues. The initial structure of heme was obtained from the Protein Data Bank (PDB) database, (PDB ID: HEM). Quantum mechanics (QM) calculation had been used to optimize the heme. The QM calculation was conducted using density functional theory (DFT) method employing Becke's (B) exchange functional combining Lee, Yang, and Parr's (LYP) correlation functional, widely known as B3LYP density functional in Gaussian 09 program package [28][29][30]. SDD basis set had been used for optimization of the heme [31,32].

Molecular Docking of Heme with the Predicted 3D
Structures of TPO. For molecular docking, Autodock Vina [33,34] protocol was employed. The molecular docking approach involved the prediction of the interaction between a small molecule and a protein at the atomic level, which provided us the opportunity to investigate the behavior of small molecules in the binding site of target proteins and to elucidate the fundamental biochemical processes. Before docking processes, knowing the location of ligand binding sites in the target protein could significantly increase the docking efficiency [35]. The binding sites of heme with TPO had been identified. It was found that Asp238, His239, Arg396, Glu399, and His494 were present in the active site of TPO 1-933 ; i.e. MPO-like domain and these amino acid residues were involved with heme interactions and thus the catalytic activity of TPO [36,37] in the x, y, and z directions, respectively. Grid box value center and grid box size was also optimized for full-length TPO   (Table 2).

Visualization and Analysis of Docking Results.
The molecular docking results of heme with TPO 142-738 WT, TPO 142-738 MT1, TPO 142-738 MT2, and TPO 142−738 MT3 and also for the whole protein TPO 1-933 were visualized and analyzed using the PyMol (version 2.0) and BIOVIA Discovery Studio 2017 software [38,39]. The binding affinities of heme with wild type TPO (TPO WT) compared to 3 mutant proteins (TPO MT1, TPO MT2, and TPO MT3) were observed. Moreover, the corresponding non-bond interactions of amino acids with heme were also studied for wild type and mutant 3D structures.

QM/MM and Molecular Dynamics (MD) Simulation.
To validate the molecular docking results, further analysis by quantum chemical method QM/MM and molecular dynamics (MD) simulations were employed for full-length TPO proteins (wild type, MT1, and MT3) since published data showed MT1 and MT3 had severe effect on enzyme activity than MT2. The QM/MM calculations of the selected proteinligand complexes were performed by a two-layer ONIOM method available in the Gaussian09 software package [40][41][42][43][44][45][46], QM and MM regions have been shown in supplementary Figure 3. The heme molecule was included in the QM region and semiempirical PM6 level of theory [47] was considered due to large structure of the heme molecules. The regions of full-length protein were computed in the MM region and the Universal Force Field (UFF) was used for the energy minimization. The total ONIOM energy of the entire system is as follows: (see [48]) The real system consists of all the atoms and is calculated only at the MM level. The model system consists of the part  of the system (such as heme) that is treated at the QM level [40]. For molecular dynamics simulation, YASARA dynamics program was employed and [48,49]. AMBER14 force field was considered for all calculations. The size of the cubic simulation box was 167.17Å * 167.17Å * 167.17Å and 151,343 water molecules were added to maintain a solvent density of 1.0 g/ml. The total number of atoms in the system was 469,240. Due to the large size of the protein, we have performed 5000 ps MD simulation [50]. For short-range van der Waals and Coulomb interactions, a cut-off radius of 8.0Å was considered. Long-range electrostatic interactions were taken into consideration using the Particle Mesh Ewald algorithm [50]. MD simulation was performed for 5000 ps at 310K having a time step of 1.25 fs and the simulation snapshots were saved at every 100 ps. To check whether the ligand stayed in the same binding pocket after 5000 ps of simulation, the heme ligands were retrieved from the simulated ligandprotein complexes and molecular docking was performed again using Autodock Vina protocol as mentioned earlier. The docked structures were analyzed using the same protocol as stated before.

Demographic Information of the Study Participants.
A total of 36 confirmed cases of congenital hypothyroidism with dyshormonogenesis were enrolled in this study. Among 36 patients, 15 (41.67%) and 21 (58.33%) were female and male, respectively, with an average age of 7.58 ± 4.56 years.

Analysis of TPO Gene for Identification of Molecular
Basis of Hypothyroidism. As mutations in the TPO gene are commonly associated with thyroid dyshormonogenesis and related complications [18,51], we opted to analyze the TPO gene to find whether there were mutations in this gene of the study participants. Upon Sanger sequencing of specimens from the patients with thyroid dyshormonogenesis targeting exon-8 to exon-14 which are commonly reported in TPO-associated thyroid dyshormonogenesis, mutations were detected in all 36 samples. A total of four mutations, namely, c.1117G>T (p.Ala373Ser), c.1193G>C (p.Ser398Thr), c.2145C>T (p.Pro715Pro), and c.2173A>C (p.Thr725Pro) were identified in the study participants ( Table 3). The first two of these four mutations were detected in exon-8, whereas the remaining two mutations were detected in exon-12. Even though c.2145C>T (or p.Pro715Pro) of the four mutations is a synonymous point mutation and is innocuous, the other three nonsynonymous mutations, namely, p.Ala373Ser, p.Ser398Thr, and p.Thr725Pro, had previously been reported in the patients with thyroid dyshormonogenesis and the reaction kinetics catalyzed by the mutant TPO enzyme proved to be similar with nonenzymatic reaction rates by several other studies [20,52,53].

Prediction of 3D Structures of Myeloperoxidase-(MPO-) Like Domain (TPO 142-738 ) and Full-Length TPO Protein
(TPO    Figure 1). To investigate the effects of the mutations on the full-length TPO protein structure and functions, we also predicted the 3D structures for  Table 4.  Table 5.

Optimization of Heme.
Density functional theory using B3LYP/SDD was employed for the optimization of the heme prosthetic group. Slight changes in bond distances and bond angles were observed between crystal and DFT structures which are presented in Table 6 and Supplementary Figure 1. Moreover, of the aforementioned 11 residues, 4 residues including His239, Arg396, Glu399, and His494 had been reported to correspond to amino acid residues in MPO protein which were crucial for enzyme activity [36,37]. For TPO 1-933 MT1, a total of 12 residues were found to interact with heme. Among those residues, 2 residues, namely, His239 and His494, interacted with heme through hydrogen bonds and 5 residues including His239, Phe243, Arg396, Phe524, and Leu567 interacted through hydrophobic interactions (Table 7, Figure 2 and Supplementary Table 2). Total number of interactions decreased significantly for TPO 1-933 MT1 compared to the TPO 1-933 WT and the Glu399 residue which is crucial for interaction was also absent in the TPO 1-933 MT1. On the other hand, when TPO 142-738 MT1 predicted structure was used for molecular docking, a total of 19 residues were found to interact with heme (Table 7, Figure 3 and Supplementary Table 1). Among those 19 residues, 4 residues including Met231, Gly234, Ser402, and Gly493 interacted with heme through hydrogen bonds, the residues Gly234, Gln235, Val400, Phe490, Arg491, Ile497, Phe523, Leu560, Leu564, and Leu575 interacted through hydrophobic interactions, and the residue Glu399 interacted through electrostatic interaction with the heme (Table 7). However, the crucial interactions of heme with His239, Arg396, and His494 were absent in the TPO 142-738 MT1. The structure-based docking of both TPO 1-933 and TPO 142-738 suggested that MT1 mutation was damaging for the TPO enzyme activity.

Molecular Docking, Visualization and Analysis of the
For TPO 1-933 MT2, a total of 20 amino acid residues were found to interact with the heme. Among those residues, 3   residues including His239, Arg491, and Arg582 interacted with heme through hydrogen bonds, 9 residues including His239, Val400, Arg491, His494, Ile497, Phe523, Leu564, Val566, and Leu575 interacted through hydrophobic interactions, and 2 residues including Arg396 and Glu399 interacted through electrostatic interactions (Table 7, Figure 2 and Supplementary Table 2). All four crucial interacting residues, namely, His239, Arg396, Glu399, and His494, which are important for MPO-like domain activity in TPO enzyme, were found to interact with heme in TPO 1-933 MT2. On the other hand, when TPO 142-738 MT2 predicted structure was used for docking, a total of 20 residues were found to interact with heme. Among those 20 residues, 3 residues including Arg491, Asn579, and Arg582 interacted with heme through hydrogen bonds, the residues Phe243, His494, Phe523, Phe524, Leu564, Phe565, Val566, Leu567, and Leu575 interacted through hydrophobic interactions and the residue Glu399 interacted through electrostatic interaction (Table 7, Figure 3 and Supplementary Table 1). However, the crucial interactions of heme with residues His239 and Arg396 were absent in the TPO 142-738 MT2. The TPO 1-933 structure-based docking showed that the mutations did not result in any major changes in interactions and all major interacting residues were present, whereas the TPO 142-738 structure-based docking suggested result which was similar to TPO 1-933 except that the major interacting residues His239 and Arg396 were absent. Together, the docking results are indicative of the fact that this mutation might have an association with PIOD. For TPO 1-933 MT3, a total of 21 residues were found to interact with the heme. Among those residues, 3 residues including Gln246, Arg491, and Ser568 interacted with heme through hydrogen bonds, 11 residues including Phe243, Val400, Arg491, His494, Ile497, Phe523, Phe524, Leu560, Leu564, Val566, and Leu575 interacted through hydrophobic interactions, and the residue Glu399 interacted through electrostatic interaction (Table 7, Figure 2 and Supplementary  Table 2). Two crucial interacting residues, namely, His239 and Arg396 which are important for MPO-like domain activity in the TPO enzyme, were absent in TPO 1-933 MT3. On the other hand, when TPO 142-738 MT3 predicted structure was used for docking, a total of 16 residues were found to interact with heme. Among those 16 residues, 4 residues including His239, Arg396, Arg491, and Arg582 interacted with heme through hydrogen bonds, the residues Phe243, Phe523, Phe524, Leu564, Phe565, Leu567, and Leu575 interacted through hydrophobic interactions and the residue Glu399 interacted through electrostatic interaction (  [20] and our molecular docking-based study showed that full-length TPO 1-933 MT2 (p.Ser398Thr) structure interacted to all the crucial amino acids in the catalytic site of TPO; thus we selected these TPO 1-933 MT1 (p.Ala373Ser) and TPO 1-933 MT3 (p.Thr725Pro) mutant cases for further analysis to validate the molecular docking results and to compare the amino acid interactions with the wild type structure. The interactions between the amino acid residues of TPO 1-933 WT, TPO 1-933 MT1, and TPO 1-933 MT3 proteins and the heme groups were further studied by QM/MM and the mode of interactions is depicted in Figure 4. It was observed that the interacting residues were the same as obtained from the initial docking results.
To investigate the structural changes during molecular dynamics simulation, the protein-ligand complex after 5000 ps of simulation was superimposed on the initial docked protein-ligand complex. The superimposed structures of  Figure 5. It was observed that the heme ligand was found within the catalytic sites.
As protein flexibility can give rise to difference in binding interactions of a ligand with its target protein, the retrieved protein structure after 5000 ps simulation was again docked with the heme ligand and the results are depicted in Figure 6. From the analysis of molecular docking after performing molecular dynamics of the protein structures, we observed that although heme interacted with TPO 1-933 WT through all the crucial amino acids (Asp238, His239, Arg396, Glu399, and His494), while it interacted with TPO 1-933 MT1 and TPO 1-933 MT3 through Glu399 and His494 residues only. Thus in the mutant cases, there were no interactions for the other 3 amino acid residues, namely, Asp238, His239, and Arg396. As these 5 amino acid residues are important for the catalytic activity of the protein, the absence of interactions with one of these crucial amino acids could affect the functional activity of the protein.

Discussion
Though congenital hypothyroidism (CH) is the most common preventable disorder, newborn screening is not a regular practice in Bangladesh. Due to late initiation of treatment, many late-diagnosed hypothyroid patients experience various typical signs and symptoms of hypothyroidism even though they receive regular levothyroxine treatment. Although several genes have been reported to be involved in thyroid dyshormonogenesis (TDH), mutation in the TPO gene is frequently described with mild to severe repercussions resulting in partial iodine organification defect (PIOD) to total iodine organification defect (TIOD) [21]. In this present study, we investigated (a) the mutational spectrum in the TPO gene of TDH patients and (b) the influence of specific mutation on TPO protein structure by means of in silico approach. To the best of our knowledge, this is the first molecular investigation on genetic etiology of TDH in Bangladesh.
In this study, we analyzed exons 8-14 of TPO gene of the TDH patients as the previous studies had reported that this region was crucial for the enzymatic activity and mutations in this region may result in absence or reduction in TPO activity [20,59]. Analysis of 36 specimens revealed three nonsynonymous mutations including p.Ala373Ser (TPO MT1), p.Ser398Thr (TPO MT2), and p.Thr725Pro (TPO MT3) and one synonymous mutation p.Pro715Pro. The identified nonsynonymous mutations had previously been reported to be pathogenic or disease-causing mutations [51,53]. Moreover, a cloning-based study involving aforementioned mutations by Guria et al. reported that these mutations could result in low expression of TPO mRNAs as well as a reduction in TPO enzyme activity [20]. They also demonstrated that mutation p.Ala373Ser (TPO MT1) was more damaging than mutation p.Ser398Thr (TPO MT2). This phenomenon could be explained by a change in aliphatic amino acid to hydroxylcontaining amino acid at 373 th position of TPO protein, whereas the other mutation (p.Ser398Thr) did not result in such a shift from one group of amino acids to another. However, the mutation in exon-12, namely, p.Thr725Pro (TPO MT3), could result in a failure of TPO protein to shift to its active state as it is well reported that threonine is the phosphorylation site for protein activation [60,61]. Moreover, p.Thr725Pro had been reported to be associated with autoimmune hypothyroidism [62]. As CH is genetically and phenotypically diverse, molecular studies may provide additional information for diagnosis, classification, and prognosis of the disease. Particularly, in patients with normal thyroid gland morphology, it could be difficult to distinguish between thyrotropin resistance and dyshormonogenesis; and molecular genetic studies may reveal true etiology of the disease in these cases. In this study, at least one of this three damaging mutations was found in a homozygous state or two of them were in a heterozygous state in each of the study participants. We previously reported a mutational hot-spot in the HBB gene and established a cost-effective molecular method (high resolution melting curve analysis) for screening of HBB gene mutations in Bangladesh [63]. Such a cost-effective approach can be adopted for TDH patients and carrier screening targeting the TPO genes in Bangladesh.
Substitution of an amino acid affects the shape, function, or binding properties of a given protein. With growing importance of genetics and genomics in health sector, considerable efforts have been devoted to linking human phenotypes to genotypic variations at the nucleotide level and associated changes in 3D protein structure [64,65]. Since the identified mutations were pathogenic, we wanted to investigate how these mutations were related to dyshormonogenesis by affecting the structural integrity and function of TPO protein. To date, there is no high resolution X-ray crystallographic structure for any of the TPO proteins (wild type or mutated) in the public databases. However, there were structures for closely related proteins, namely, myeloperoxidases (MPO) and lactoperoxidases (LPO) [11,12]. Thus we predicted and validated the three dimensional (3D) structures of TPO (WT and MT) protein using various bioinformatics tools. We studied the effects of mutations on the structural integrity, arrangement of various domains, and folding patterns of the TPO protein. To further study mutational effects, we performed molecular docking of heme in the catalytic site of TPO to investigate how these mutations could affect the activity of TPO enzyme. We obtained the binding energies of heme for the wild type and the mutant structures of TPO protein and found a reduction in binding affinities in the mutant structures compared to the wild type one. Further, we analyzed the detailed results of molecular docking by observing the non-bond interactions of heme with specific amino acid residues to understand the effects of mutations on the functions of TPO protein.
In our study, we applied our bioinformatics approaches on the MPO-like domain of TPO (TPO 142-738 ) and full-length TPO protein (TPO 1-933 ) targeting the non-bond interactions between the heme group and amino acid residues as the mutations identified in this study were found in this region and this MPO-like domain was the catalytic site of TPO enzyme [54]. We obtained the predicted structures for the wild type and the mutant proteins from the I-TASSER server for TPO 142-738 and all the structures had significant confidence score (C-score) suggesting that the structures were valid for further studies [22]. However, I-TASSER predicted structures for TPO 1-933 had lower C-score value and the reasons may be due to prediction for a very large protein, because I-TASSER predicts structures which are based on iterative threading and also homology modelling and as TPO has no crystallographic structure, this still remains a challenge [22,24]. We also verified the structures using the Verify3D and the RAMPAGE web servers and both of the servers gave satisfactory results [27,66]. As heme is crucial for the catalytic activity of TPO, the molecular docking of heme with the wild type and the mutant TPO structures could help us to understand the effects of mutations on the functions of TPO protein. We observed a decrease in heme binding energies in the cases of mutant TPO proteins suggesting that the catalytic activity of mutant TPO might have been hampered. Further investigation gave us information about the heme interactions with specific amino acid residues of TPO protein.
The heme prosthetic group of wild type structure interacted with all the important amino acid residues including His239, Arg396, Glu399 and His494 through non-bond interactions.
But in the cases of mutant structures, some of the important amino acid interactions were absent, suggesting that the mutations might have damaging effects on heme interactions and thereby affecting the catalytic activity of TPO enzyme. Guria et al. demonstrated that MT1 (p.Ala373Ser), MT2 (p.Ser398Thr), and MT3 (p.Thr725Pro) had damaging effect on TPO mRNA expression and enzyme activity [20]. However, MT1 (p.Ala373Ser) and MT3 (p.Thr725Pro) showed iodination reactions similar to the nonenzymatic reaction rate, suggesting that these two mutations were more damaging than the MT2 (p.Ser398Thr) which showed more efficient iodination reaction than MT1 and MT3 but less efficient than the wild type TPO protein [20]. This finding is consistent with our molecular docking-based study targeting both the MPO-like domain of TPO (TPO 142-738 ) and the full-length TPO protein (TPO 1-933 ) predicted structures. MT1 and MT3 had more enhanced influence on the interactions between several crucial residues and the heme group, whereas MT2 had less influence on the interaction between the TPO protein and the heme prosthetic group. Different studies showed that quantum chemical methods such as DFT or QM/MM can also be employed for investigating the intermolecular interactions [67][68][69]. Thus, to validate the results of molecular docking, we performed QM/MM calculations on the fulllength wild type TPO protein and two mutant proteins (MT1 and MT3) that had severe damaging effect on the catalytic activity of TPO.
The heme group was treated using a semiempirical approach of calculations and the protein was treated using an approach of molecular mechanical calculations. The structures obtained after performing QM/MM calculations were taken into consideration for further analysis of the presence of non-bond interactions of the active site residues with the heme group. The interactions observed were exactly similar to those obtained from molecular docking, suggesting that the docking results obtained were valid.
In cellular environment, protein structures are always in a dynamic condition. To mimic the cellular environment, we performed molecular dynamics simulation in a cubic simulation box containing water molecules and NaCl and a pH of 7.4 was maintained. After performing 5000 ps simulation under the influence of AMBER14 force field, the final protein structures were superimposed on the corresponding docked protein structure. The results infer that the proteinligand complexes studied were considerably stable over time as there was very small change in the coordinates of the heme group after the simulation. To observe the effect of protein dynamics, the protein structures were retrieved from the 5000 ps simulation snapshot and again docked with the heme ligand. The finding was rather interesting as the wild type protein was still found to be interacting with the crucial amino acid residues. However, these important interactions were found to be absent for the mutant proteins that might be a potential cause of a decrease in the enzymatic activity of TPO protein.
To understand the in-depth characteristics of TPO protein of both wild type and mutant forms, we performed in silico approach to mimic cellular environment. However, various computational chemistry-based approaches like IR and Raman spectroscopy can be used to identify the changes in molecular level with higher specificity [70,71]. Such compositional analysis in various cellular conditions has become very popular for characterizing the biochemical changes in various disease conditions and also for the study on characterizing the spectrum of various hormones such as corticosteroids [70][71][72]. The study by Claudio et al. analyzed the molecular vibrational spectrum of thyroid tissues from normal and disease conditions which could ultimately represent the characteristics of secondary structures of proteins [73]. Although use of IR and Raman spectroscopy could offer more insights into the physiological conditions of TDH patients carrying mutations in the TPO gene, the present study was not subjected to such approaches because we did not have such facilities.
In this present study, we had used computational approaches such as molecular docking, QM/MM and molecular dynamics simulation to investigate the effect of mutations on TPO protein interactions. The present study could be useful for future studies including study of the effect of mutations on TPO dimer organization and how mutations in the TPO gene could lead to a change in the TPO protein conformation.

Conclusion
Our study investigated the genetic etiology of Bangladeshi patients with TDH, which may further help us to screen and categorize the disease. Three pathogenic mutations were observed in the patients with TDH including p.Ala373Ser, p.Ser398Thr, and p.Thr725Pro. In silico based study revealed that p.Ala373Ser and p.Thr725Pro mutations resulted in a significant change in interactions between the amino acid residues of TPO protein and the heme group, whereas the p.Ser398Thr mutation could affect the TPO protein to a lesser extent. The results of this study may help better understand the correlation between specific mutation in the TPO gene and altered biological activities of the TPO protein as well as disease severity among the TDH patients. Furthermore, future study on dimer formation of TPO protein and functional activity study of TPO enzyme in physiological condition is expected to shed light on how mutation in the TPO gene can affect thyroid hormone biosynthesis pathway and find a mutation-based better treatment strategy for individual patients.

Data Availability
The data used to support the findings of this study are included within the article.

Ethical Approval
This study was approved by the Ethical Review Board for Human Studies of BSMMU and University of Dhaka.

Conflicts of Interest
The authors declare that they have no conflicts of interest.