Structural changes in thyroid hormone receptor-beta by T3 binding and L330S mutational interactions

1.   Introduction

The point mutations on thyroid hormone receptor-beta (THR-β) isoform of the nuclear receptor super family cause generalized resistance to thyroid hormones (RTH). The genetic disorder, RTH has the clinical features of elevating thyroid hormone levels without suppressing thyroid stimulating hormone (TSH) in the blood serum. The RTH patients suffer from inadequate secretions by hypothalamus-pituitary-thyroid (HPT) axis. The RTH reduces triiodothyronine (T3) binding affinity to THR-β ligand binding domain (LBD) and it further resists the binding affinity to DNA binding domain resulting in gene repressions. The wild type THR-β is able to select the coactivators for the basal transcription whereas mutated THR-β recruits corepressors inhibiting the transcription factors. Thus, being linked with genetic defects, the THR-β mutations cause neoplasia and several endocrine disorders. Goitre formation, dizziness, speech disorder and tachycardia are common causes of RTH. All of these clinical facts are associated with the previously identified THR-β mutations [1]–[5]. L330S was identified in an 11-year-old boy with RTH showing the symptoms of both hypo- and hyperthyroidism, goitre and speech disorder [6]. The L330S mutation was also detected in a 19-year Thai woman with RTH and goitre [7].

To explore the mutational impacts on THR-β gene of the patients with RTH, it is necessary to compare the structural and physical properties of the mutated THR-β with that of wild type liganded and unliganded THR-β associated with euthyroid and hypothyroid states. For this purpose, we can perform MD simulations of each system and analyze root mean square deviation (RMSD), radius of gyration (RG), root mean square fluctuation (RMSF), classical as well as quantum IR spectral density (IRSD), radial distribution function (RDF), solvent accessible surface area (SASA), Ramachandran plots for the dihedral-angle distributions and interaction or internal energies. These types of analysis were also performed by Khan, et al., 2016 to explore the mutational effects in Ras-related protein [8] and by Rajendran, et al., 2012 to investigate the molecular mechanism of laminopathy due to a point mutation R482W in lamin A/C protein [9]. The SASA calculations of the residues in protein folding and unfolding states ranging from 0–300 Ų were reviewed by Ausaf Ali, et al., 2014 [10]. The protein stability influenced by buried and surface point mutations is due to physicochemical, energetic, and conformational properties of amino acid residues as reflected by the mutant positions on Ramachandran plots [11]–[13].

2.   Methodology

The wild type (WT) native structure of T3-liganded thyroid hormone receptor-β (THRT3-WT) was taken from the crystal structure code: 3GWS.pdb [14]. The mutated (MT) nuclear receptor (THRT3-MT) was prepared after the point mutation on 330-codon replacing LEU by SER with the help of PyMOL-Mutagenesis. Topology and parameters required to prepare the complete simulated systems were used from the CHARMM force fields [14]–[17]. The T3 hormone being the hetero-molecule, it was suitably parameterized to define numerical constants, atomic masses and charges required to evaluate forces and energies [18]. Starting from the generation of protein structure files (psf), all the simulations for energy minimization, equilibration and production runs were performed with nanoscale molecular dynamics (NAMD-2.12) [19]. The simulated systems were analyzed with visual molecular dynamics (VMD-1.9.3) [20].

Each of the molecules, i.e. unliganded THR-WT, THRT3-WT and L330S-mutant THRT3-MT was solvated with 17220 water molecules (TIP3P) under the spherical boundary conditions and neutralized with Na⁺ and Cl⁻ ions in the concentration of 0.15 mol/L. Each system had the boundary conditions: radius 37.48 Å, center coordinates (4.14, 25.05, 20.60 Å), force constant 10 kcal mol⁻¹Å⁻², and exponent 2 for the first potential to be applied. The TCL script file called configuration file was prepared to provide the simulation details. The Langevin piston temperature is 310 K and damping coefficient is 1 ps⁻¹. The Langevin dynamics [21],[22] was implemented with the integrator parameter of 2 fs/step keeping rigid bonds of all H-atoms. Multilevel summation method (MSM) was active with the grid spacing of 2.5 Å for the electrostatic force evaluation [23]. In order to evaluate the electrostatic and van der Waals interactions, the force field related terms such as 1–4 scaling, cut-off, switching and pair-list distances were set as 1.0, 12 Å, 10 Å and 14 Å, respectively. Then, the energy of each system was minimized up to 3000 CG steps and its equilibration run was performed up to 20 ns at the body scale temperature 310 K. The new trajectories of the constituent atoms were calculated by using the velocity Verlet algorithm [24].

The conformational changes in the unliganded, T3-liganded and mutated THR-β were analyzed from the nature of RMSD, RG, RMSF and SASA plots related to their protein-backbone, T3-hormone and residues of interest, i.e. L330 and L330S. Ramachandran plots were used to observe the distributions of dihedral angles and the possible steric hindrances among α-helices, β-sheets, side chains and residues during the MD simulations. With the help of NAMD Energy Plugin (Version 1.4), the interaction energies (coulombic and van der Waals) between the native or mutated residues on the 330 codon and T3-hormone were observed in the related system. Furthermore, H-bond distribution, RDF, IRSD and internal energy of T3 and L330 codons were analyzed in both wild type and mutated systems.

As per analysis of the structural trajectory, RMSD vs time graph has been plotted to explain the conformational stability of the molecular system during the MD simulations. The RMSD is a numerical measure of the difference between the protein structures defined by Equation 1 [25],[26].

where N is the number of atoms whose positions are to be compared, r_i is the position of atom i and r_{i, ref} is the position of atom i in the reference configuration defined by X-ray structure of the protein system.

In addition to RMSD, the equilibrations of the system have been analyzed with RMSF and RG throughout the simulations. The RMSD from the average over time can be referred to as the RMSF. The RMSF of residue i is given by Equation 2 [26].

where N_t is the number of time steps giving new configurations in the simulations, r_i(t) is the position of residue i, ⟨ri⟩=1Nt∑tj=1Ntri(tj) is the average position of residue i at time t_j.

The RG is a measure of compactness describing the extent of collapse or extension of the protein backbone [26],[27]. It represents the deviation of atoms in a molecule from its centre of mass defined by the Equation 3.

where N is the total number of protein atoms, r_i is the position of atom i, rcm=1M∑t=1Nmiri is the centre of mass, m_i is the mass of atom i and M is the total mass of the protein atoms.

Solvent accessible surface area (SASA) is the area described by the locus of the center of the water molecule when it rolls along the protein complex by making the maximum permitted van der Waals contacts without penetrating any other atom. The SASA is calculated by using the Equation 4 [10], [28].

where D = ΔZ/2 + Δ′Z, R = r + r_s = sum of van der Weals radius of the atom and chosen radius of the solvent molecule, L_i = length of arc drawn on a given section i, Z_i = the perpendicular distance from the center of the sphere to the section i, ΔZ = the spacing between the sections, and Δ′Z = ΔZ/2 or R−Z_i, whichever is smaller. The summation ranges over all of the arcs drawn for the given atom. Here, the solvent accessibility is defined by SA = 100.SASA/4πR².

In this way, NAMD calculates the spherical atomic radial distribution function (RDF/g(r)) between the atom coordinates in two selections over a given trajectory using GUI Plugin (Version 1.3) in VMD as explained by Levine, et al., 2011 [29]. The normalized distance-dependent RDF for the atom pairs with types i and j is defined by the Equation 5 [30].

where scaling to 4πr² accounts for the volume of the spherical shell of radius r and thickness dr. The number N_i,j(r) of atom pairs i,j at a distance between r and r+dr is found by counting the occurrences Ni,j(r)=∑i∑rδ(|ri−rj|,r) where the delta function is δ = 1 if r ≤ |r_i − r_j| ≤ r+dr and otherwise it is zero. Here, the bin size dr is chosen with high resolution between r_min and r_max including the atom pairs.

NAMD computes spectral densities from trajectories using time series data created with the measure dipole command in VMD-GUI, i.e. IR Spectral Density Calculator Plugin (Version 1.3). The details of theory and methods of NAMD simulations to calculate vibrational spectra of proteins were explained by Malolepsza, et al., 2014 [31]. The classical IR spectral density, I_cl(ω) is given by Fourier transform of the auto-correlation function G_cl(t) of the dipole moment operator [32].

where G_cl(t) = tr{ρ(β)μ(0)μ(t)†}. Here, μ(0) is the dipole moment operator at the initial time and μ(t) is the same operator at the time t, and ρ(β) is Boltzmann density operator with the statistical parameter β = 1/K_BT.

3.   Results and discussion

The protein structure of wild type THR-β that we take for the simulation has 3895 atoms with total mass 27640 amu. The active ligand bound to THR-β LBD is T3 hormone which has 35 atoms with mass 651 amu. In the point mutation L330S, LEU (C₆H₁₃NO₂) of mass 131 amu is replaced by SER (C₃H₇NO₃) of mass 105 amu.

The constancy of RMSD (≤ 3 Å) and RG (≈ 19 Å) of the protein backbone over the time frames of equilibration runs show the conformational stability and validity of force field topology and parameters used for MD simulations of wild type and mutated receptor systems. Some residues with more fluctuating RMSD/RG show their dynamic properties which are responsible for T3 binding and dissociation from THR-β LBD. The values of RMSD and RG of the molecules with the standard deviations of the data are listed in the Table 1. During the 20 ns equilibration runs in neutral water-ion environment, RMSD of 330-residue in liganded THRT3-WT is smaller and less fluctuating than that in unliganded THR-WT and L330S-mutant THRT3-MT (Figure 1). T3 is more fluctuating due to smaller RG of L330S in mutated LBD than in wild type LBD. The RMSD/RG of the protein backbone is almost same with small fluctuations in these three cases (Table 1).

Root mean square fluctuation (RMSF) of the mutated residue L330S is smaller than the wild type L330 provided with larger fluctuations in the unliganded system as shown in the Figure 2. The RMSFs of all residues in unliganded, liganded and mutated THR-β systems are small lying below 1.0 Å during their equlibrations as shown in the Figure 3. Related to the 330-codon, the RMSFs are 0.46 ± 0.10 Å for LEU of THR-WT, 0.41 ± 0.07 Å for LEU of THRT3-WT and 0.15 ± 0.04 Å for SER of THRT3-MT (Table 2). The positions of L330S mutational sites and T3-hormone in THRT3-MT after their 20 ns equilibrations are shown in Figure 4. During the MD simulations of wild type and mutated THR-β systems, distance between L330 and T3 in THRT3-WT varies as 6.36 ± 0.33 Å and the distance between L330S and T3 in THRT3-MT varies as 7.11 ± 0.28 Å. Such distance is smaller and more fluctuating in THRT3-WT than in THRT3-MT (Figure 5). This is one of the indicators of T3 resistant features of mutated nuclear receptors causing RTH.

Radial distribution function (RDF) or pair auto-correlation function (g(r)) in a system of particles measures probability of finding a particle at r distance away from the reference particle [33]. The distribution function g(r) between the iodine atom I1 (T3) and the oxygen atom OG (L330S) in THRT3-MT is shifted left from that in its wild type THRT3-WT within the range of 3.0–5.5 Å. In this way, the RDF between I3 (T3) and CG2 (I431V) in THRT3-MT has smaller height within the flat range of 4.5–8.5 Å than that in its wild type THRT3-WT as shown in the Figure 6. In other words, the number density of the mutant residue is less than that of the native wild type residue in order to respond for T3-hormone and TREs.

Ramachandran plots of dihedral angles (φ, ψ) for THR-WT, THRT3-WT and L330S mutant THRT3-MT in the last frame of 20 ns equilibrations are shown in Figure 6. The residues of parallel, anti-parallel and right twisted β-sheets lie in the blue region of I; right handed α-helix lie in the blue region of II, and left handed α-helix lie in the green region of III quadrant [34],[35]. The amino acid GLY existing without side chain has no steric hindrance so that it is allowed in the white region of IV quadrant. All the amino acids except GLY are sterically disallowed in the white regions. The most allowed regions for α-helices and β-sheets are the blue regions. The green region is allowed for those amino acids which consist of atoms with short van der Waals radii. By comparing the plots of Figure 7, we observe that the (φ, ψ) values of L330S in THRT3-MT scatter more towards the disallowed regions than that of L330 in THRT3-WT. This analysis verifies that the larger steric hindrances occur in the unliganded and mutated THR-β resulting gene repressions and the related thyroid disorders. The clinical evaluations show that THRT3-WT, THR-WT and THRT3-MT are associated with normal thyroid functions, overt hypothyroidism and RTH, respectively.

The SASA of 330-codon in THRT3-WT is observed to be different in the extension radius of 1.4 Å than that in THRT3-WT and THRT3-MT (Figure 8). As depicted in the Table 2, The SASA of the 330-codons in THR-WT, THRT3-WT and THRT3-MT are 290.96 ± 3.51, 285.62 ± 3.74 and 223.89 ± 3.41 Ų, respectively which lie in the range of 0-300 Ų reported by Ausaf Ali, et al., 2014 [10]. The lower SASA indicates the higher thermodynamic stability of the protein. The SASA of T3-hormone is observed to be zero indicating the hydrophobicity of LBD of THR-β. The mutation on the solvent accessible residue causes the conformational change by partial unfolding of the protein [36]. Such effect is more intensive if the mutation occurs on the surface residue than on the interior of the protein. Since the SASA of L330 has been found to be decreased relative to its wild type, such mutation is confirmed as a buried mutation of the THR-β gene.

The IR spectral densities using total dipole moment of the atoms in the 330-residue are calculated in the frequency range of 0 to 2500 cm⁻¹ by adjusting the time interval between frame to frame calculations as 0.01 ps in the last 1 ns of the equilibrating simulations. Such spectral distributions are different for THRT3-WT and THRT3-MT as shown in Figure 9. The maximum IR-intensity of L330 in the wild type receptor is 0.0091 at 1228.56 cm⁻¹. Again, the maximum IR-intensity of L330S in the mutated receptor is 0.0081 at 1163.55 cm⁻¹. The distribution has larger number of high intensity peaks for L330 in THRT3-WT than for L330S in THRT3-MT.

Total non-bonding interaction energy is the sum of electrostatic and van der Waals interaction energies. The changing values of interaction energy between 330-codon and T3-hormone in wild type and L330S mutant THR-β systems are plotted w.r.t. the frame order in the Figures 10 and 11. These are averaged over the time frames from the last 1 ns of the equilibrating simulations. The mean ± SD of electrostatic (E_elect) and van der Waals (E_vdw) interaction energies are calculated by using the NAMD-Energy-GUI adjusted with switching and cut-off distances of 10 Å and 12 Å, respectively. As per the interactions between L330 and T3 in THRT3-WT, the values of E_elect and E_vdw are −1.58 ± 0.72 kcal/mol and −4.22 ± 0.84 kcal/mol whereas for the interactions between L330S and T3 in THRT3-MT, these are −2.37 ± 1.52 kcal/mol and −2.26 ± 0.85 kcal/mol, respectively (Table 2). The Figures 9 and 10 also demonstrate that E_elect is more fluctuating with larger negative value whereas E_vdw is less fluctuating with smaller negative value when LEU is mutated with SER in 330-codon of the THR-β gene. These factors are also significantly important to explore the insights into the RTH.

Internal energy of the biomolecular system is the sum of bonding and non-bonding potentials. Coulomb energy is the largest and its fluctuation during MD simulations makes the total internal energy fluctuate about its mean value. Using NAMD-Energy-GUI with switching and cut-off distances of 10 Å and 12 Å, the internal energy (U_i) of 330-residue, bounded T3-hormone and protein of THR-WT, THRT3-WT and THRT3-MT are calculated from the last 1ns of their equilibrations and the average values for the corresponding systems are listed in the Table 2. The variations of internal energy w.r.t. the simulation time are plotted in Figure 12. The plots for the individual potential terms show that the major contributor of fluctuations in internal energy is the electrostatic energy of the related system. The value of U_i is 30.03 ± 3.03 kcal/mol for L330 in THRT3-WT and 43.70 ± 2.25 kcal/mol for L330S in THRT3-MT, i.e. the internal energy increases when LEU is mutated with SER in the 330-codon of THR-β LBD. In addition to these mutational impacts, thermophysical properties such as temperature echo, heat capacity, thermal diffusivity and thermal conductivity are also important to evaluate hormone-receptor interactions in different THR isoforms [37],[38].

4.   Conclusions

Almost constant values of RMSD (≤ 3 Å) and RG (≈ 19 Å) of wild type as well as mutated THR-β show their conformational stability during MD simulations. The point mutation L330S with largely fluctuating RMSD shows its resistive nature to the T3 binding in and dissociation from THR-β LBD. The (φ, ψ) values in Ramachandran plots of L330 in unliganded THR-WT and L330S in mutated THRT3-MT scatter more towards the disallowed regions than that of L330 in liganded THRT3-WT. This analysis verifies that the larger steric hindrances occur in the unliganded and mutated THR-β as the basis of gene repressions and the related thyroid disorders. The smaller SASA of L330S mutation than its wild type residue indicates that it is stable buried-muation. The SASA of T3 is zero verifying its hydrophobicity in THR-β LBD. Comparing RDFs between the atoms of T3 and mutational residues, the probability of finding the atomic particles within the domain range of 3 to 10 Å is shifted closer in THRT3-MT than in THRT3-WT. The IR spectral distributions in the range of 0-2500 cm⁻¹ suggest that the L330-residue in THRT3-WT has larger number of high spectral peaks with greater change in dipole moments than the L330S-residue in THRT3-MT. Related to the 330-codon of THR-β LBD, electrostatic interaction energy is more fluctuating with larger negative value, van der Waals interaction energy is less fluctuating with smaller negative value and total internal energy is less fluctuating with larger positive value when LEU is mutated with SER. These are the strong parameters that explain the nature of point mutations in THR-β gene.