The reported FTO SNPs and clinically important mutations were collected from dbSNP and ClinVar databases, respectively (24). The information about the functionally important mutations at the protein level was obtained from the UniProt database (25). Different mutations obtained from all databases were further verified by literature resources (8, 9, 26, 27). Tertiary structures of FTO protein are available at the protein data bank (PDB) with PDB codes: 3LFM and 5ZMD (21, 28). Since available structures at PDB are incomplete and have missing amino acids up to 35, therefore, we adopted homology and de novo modeling approaches to construct the full-length model using I-TASSER server (29).

Sequence-based pathogenicity prediction of various mutations was performed by the PredictSNP server which facilitates eight prediction tools: MAPP, PhD-SNP, PolyPhen-1, PolyPhen-2, SIFT, SNAP, PANTHER, and PredictSNP (30). It predicts the consequences of mutations on protein function based on different parameters inferred from evolutionary, physicochemical, or structural characteristics as previously described (22). Stabilities of different missense mutations were predicted through the Mutpred2, I-Mutant2.0, MUpro, mCSM, and iStable tools. MutPred2 predicts the pathogenicity of amino acid substitution and describes the molecular alterations likely altering the phenotype. It assesses the structural, functional, and phenotypic consequences of sequence variants. MutPred2 score > 0.5 indicated harmful mutation (31). I-Mutant2.0 is based on the support vector machine (SVM) that predicts protein stability upon substituting specific amino acid. I-Mutant2.0 predicts the direction of the free energy change (ΔΔG) values (kcal/mol). ΔΔG < 0 indicates the variant decreased the protein stability while a ΔΔG > 0 indicates the variant elevated protein stability (32). MUpro is also based on both the SVM and neural networks that predicts the effects of single-site mutation on protein stability. Free energy changes (ΔΔG) < 0 and ΔΔG > 0 indicate the variant with decreased and increased the protein stabilities, respectively (33). mCSM is a mutation cutoff scanning matrix tool which uses graph-based signatures to predict the impact of mutations on protein stability (34). Further, iStable server was also used to assess the stability of different mutants, and functionally important residues were predicted through an algorithm based on sequence conservation (35, 36).

The tertiary model of FTO obtained from I-TASSER was considered a wildtype (WT) and different mutant (MT) models were constructed through the PyMOL graphical software (The PyMOL Molecular Graphics System, Version 1.3 Schrodinger, LLC). Both WT and MT models were initially energy minimized by GROMACS (Groningen Machine for Chemical Simulation) suite (37) and later the structure quality of all models was examined through the PROCHECK, ProSA, ERRAT, and QMEAN servers (38–41).

Molecular dynamics (MD) simulations were performed using the GROMACS program (37). The Amber03 force field was applied to generate the topologies of WT and all MTs (42). All the protein systems (WT and MTs) were solvated in a triclinic cubic box using the TIP3P (transferable intermolecular potential 3P) water model, under the defined buffer system of 0.12 nm distance between the edge of the protein to the surface of the box. Each system maintained electroneutrality with the addition of sodium (Na⁺) and chloride (Cl^–) ions. Once the systems were neutralized, the energy minimizations were followed through the steepest descent (SD) method. Two common ensembles, NVT and NPT were employed using positional restraints for 100 and 500ps under a constant temperature of 300K and pressure of 1bar, respectively. The linear constraint solver (LINCS) algorithm was used to constrain all the bonds. The temperature and isotropic pressure of the system were maintained using modified Berendsen thermostat (v-scale) and Parrinello-Rahman methods, respectively. Further, the well-equilibrated complex molecules were subjected to a production run of 100ns time for each complex system in duplicate, by applying 2fs time steps. The analyses of trajectories were done through various utilities of GROMACS such as gmx energy, gmx rms, gmx rmsf, and gmx gyrate.

To obtain a precise view of the prominent motions in a trajectory generated during MD simulation, we employed essential dynamics (ED) also known as principal component analysis (PCA). ED is a computational approach to understanding the correlated motions of a given protein that relate to its biological functions (43). The covariance matrix was constructed after eliminating the rotational as well as translational motions and diagonalized to generate a set of eigenvectors with respective eigenvalues. ED analysis was performed by gmx covar and gmx anaeig utilities of the GROMACS suite which were used to calculate and analyze the eigenvectors and their corresponding eigenvalues. Eigenvectors also known as principal components (PCs) help us to gain insight into the direction of concerted motions of an atom while its corresponding eigenvalues represent the magnitude of displacement. ED analyses were confined to backbone atoms to avoid statistical noise and motions owing to biological events were measured by assessing the cosine content of respective PCs (44).

The functions of protein rely on its residues and interactions among them or with their surrounding environment. Each tertiary structure of protein exhibits its territory of network. To elucidate how non-synonymous mutations influenced the interatomic synergy, we studied residue interaction network (RIN) using the RINalyzer tool (45). Initially, MD optimized models were used as inputs to establish networks based on their physicochemical properties and their centralities. Network centrality was measured through weighted graph theory in which nodes represent residues and weight stand for the number of hydrogen bonds among the nodes. The obtained networks were analyzed and visualized in Cytoscape3.4 (46). We evaluated betweenness centrality (C_B) for both WT and MTs which helped us to identify the essential residues imparted during signaling.

Molecular docking is a key approach to studying the interaction between protein and small molecules such as ligand and the dynamics of that molecule within the binding pocket of the protein. Here, we used ensemble docking with 5 docking tools such as AutoDock Vina (47), QuickVina2 (48), SMINA (49), LeDock (50), and GNINA (51) and ligand binding affinity was measured in kilocalorie per mole. For all docking tools, except LeDock, the receptor and ligands were prepared in AutoDock tools as described previously (52, 53). Briefly, the receptor (WT and all MTs) and ligand were prepared by adding the polar hydrogen and Gasteiger charge, respectively. Dimension and coordinates for the ligand binding site were taken from the experimentally known binding site available at the PDB (5ZMD) structure (21). To cover the binding pocket, we created the grid box with the dimensions of 28 × 28 × 28 and centered at 80.533, 76.699, and 74.23 of x, y, and z coordinates, respectively. Receptor for LeDock was prepared in LePro module of LeDock tool and dimensions of grid box were set with the same size as mentioned above. Average docking scores and complexes with best binding poses for WT and all MTs were further analyzed. 3D and 2D plots of protein and protein-ligand systems were rendered in PyMOL (The PyMOL Molecular Graphics System, Version 1.3 Schrodinger, LLC), Schrodinger Maestro 12.8 (Maestro, Schrodinger, LLC, New York, NY, United States, 2021), and LigPlot + tools, respectively (54).

The variants or mutations of FTO gene investigated in this study were retrieved from the dbSNP database. The database exhibited approximately 94,000 SNPs out of which about 91,312, 2,133, 383, 169, 2, and 1 SNPs were intronic, non-coding transcript variants, non-synonymous (missense), synonymous, inframe deletion, and initiator codon variants, respectively (Figure 1A). For further analysis, we have carried missense mutations as these mutations led to changes in amino acid residues that directly affect the structure and function of a protein. Further, the percentages of transition and transversion within the missense mutations of the FTO gene were inspected and it was found that the G > A (20.6%) transition showed the highest prevalence followed by A > G (17.41%), C > T (15.83%), and T > C (8.17%) (Figure 1B) while transversion such as A > C (6.59%) and C > G (6.59%) contributed the highest percentage followed by C > A (5.8%), G > C (5.01%), T > A (4.22%), T > G (3.69%), G > T (3.43%), and A > T (2.63%) (Figure 1B).

After removing the redundancy, about 245 missense variants were further used for the PredictSNP server to examine the deleterious nature of mutations (Figure 1C and Supplementary Table 1). Out of 245 variants, 7 variants (R96M, Y108A, F114D, E234P, R316Q, R322Q, and C392D) were finally selected based on higher frequencies of deleterious predicted by all 8 SNP tools (PredictSNP, MAPP, PhD-SNP, PolyPhen-1, PolyPhen-2, SIFT, SNAP, and PANTHER) (Supplementary Table 2). Additionally, the pathogenic strength and stability of 7 variants were further inspected through multiple independent tools such as I-Mutant, MuPro, mCSM, and iStable. Mutpred2 showed out of 7 variants, 6 were scored higher than 0.5 suggesting the pathogenic nature of these mutants. The highest score was observed in the F114D (0.926) variant followed by Y108A (0.926) and R96M (0.921), respectively (Table 1). ΔΔG prediction by I-mutant showed that all the selected nsSNPs or variants showed decreased stability (ΔΔG < 0) (Table 2). However, the lowest stabilities were found in the F114D and Y198A variants with ΔΔG values −3.13 and −2.44, respectively. MU-Pro analysis suggested that all the selected 7 variants were found to decrease protein stabilities, where Y108A (ΔΔG: −2.20) and F114D (ΔΔG: −1.78) showed the highest protein disabilities. According to the MU-Pro and mCSM analyses, all the selected variants were found to decrease stabilities except for the E234P variant which has a stabilizing property as predicted by mCSM. In mCSM, destabilizing was observed in all variants. Among them, Y108A (ΔΔG: −2.696) and F114D (ΔΔG: −2.694) variants showed a highly destabilizing nature (Table 2). The iStable server showed that all 4 MTs, namely E234P, R316Q, R322Q, and C392D exhibited decreased stabilities. The 3 variants (R316Q, R322Q, and C392D) unanimously showed ΔΔG < 1 kcal/mol calculated by three tools (I-Mutant, MuPro, and iStable), which would predict to disturb the structure and function of the protein. Functional impacts of R96M, Y108A, and F114D MTs were previously done through a structural-based experiment while the functional impact of 4 MTs such as E234P, R316Q, R322Q, and C392D remains to be explored (21). Therefore, in this study, we have elucidated the structural and functional impacts of 4 MTs namely E234P, R316Q, R322Q, and C392D on the FTO protein. Further, sequence conservation analysis showed that the Glutamate, Arginine, and Cysteine residues positioned at 234, 316, 322, and 392, respectively, exhibited ≥ 0.8 scores and thus indicated that the above residues are functionally conserved (Supplementary Table 3).

FTO is ∼559 amino acid long protein (Uniprot: Q9C0B1) and its crystal structure is available at PDB which was resolved at ∼2.5Å resolution. However, the above reported structure has missing residues at flanking ends of N- (∼31) and C-terminals (∼62). Thus, we computationally modeled the full-length tertiary structure of FTO using the I-TASSER server service. The predicted tertiary model has −0.59 and 0.64 C- TM-scores, respectively, which indicate that the model is reliable and has a better global topology. Further reliability of the generated model was verified by superposing the template and model (Supplementary Figure 1). We omitted ∼31 and 38 residues long extreme N- and C-terminals as these residues formed highly disordered structures and also these regions do not affect the conformation and catalytic activity of the protein. The tertiary structure of FTO comprised 9 α-helices and 14 β-sheets. Double-stranded distorted β sheets reside solely in NTD along with 4 α helices (α1–α4) while CTD consists of 5 unique helices (α5–α9). NTD controls core catalytic activity with the involvement of CTD and encompasses β5–β14 (Figure 1D). CTD is primarily made up of α-helices in which α7, α8, and α10 formed a three-helix bundle. One end of the helix interacts with NTD thus involved in the stabilization of NTD conformation. Mutations in the CTD domain greatly affected the conformation and catalytic activity of the protein (28). Out of 4 selected mutations, 3 mutations (E234P, R316Q, and R322Q) were located at NTD while 1 mutation (C392D) was located at CTD (Figure 1E). The tertiary model obtained from I-TASSER was assigned as wildtype (WT) and structures of all 4 mutants (MTs) were prepared in the PyMOL tool using WT as templates. All WT and MT models were further subjected to quality assessment by evaluation of their 3D geometry and stereochemical properties. The model quality evaluation results suggested that WT and all MTs had a negligible number of amino acids at the disallowed region of Ramachandran plots and had high structure quality as indicated by > 85% of ERRAT scores for WT and all MT models (Supplementary Table 4).

The sequence level prediction of protein stability was examined through various tools as discussed above. The protein stabilities at the structure level of WT and all MTs were assessed through MD simulations which were performed in duplicate for 100 ns each (Supplementary Figure 2). During MD simulations various parameters such as RMSD, RMSF, and Rg were measured concerning equilibrium structures. RMSDs of WT and MTs were thoroughly stabilized within a time scale of 100 ns (Figure 2A). WT and all MTs showed steady behaviors of RMSDs after ∼60ns time. Hence, well-equilibrated trajectories were utilized for further analysis in all cases (Figure 2B). The average RMSD values of WT, E234P, R316Q, R322Q and C392D were approximately 0.50, 0.44, 0.39, 0.54, and 0.56 nm, respectively (Table 3). E234P and R316Q MTs showed higher RMSDs than WT with 13 and 22 nm of percent differences, respectively, while R322Q and C392D MTs showed lower RMSDs than WT with percent differences of −8 and −12 nm, respectively. The Above data showed that RMSDs of all MTs displayed stable behaviors with more or fewer RMSD values as compared to WT indicating that E234P, R316Q, R322Q, and C392D MTs largely affected the structure of the FTO protein. The effect of different mutations on the protein compactness and size was monitored by measuring the Rg at the function of time. WT and all MTs showed a similar pattern of Rg with average values being ∼2.64, 2.62, 2.62, 2.63, and 2.63 nm for WT, E234P, R316Q, R322Q, and C392D, respectively (Supplementary Figure 3). All MTs displayed similar Rg values with WT indicating that the compactness of the FTO protein remains unaffected upon these mutations.

Next, we wanted to check the stabilities and flexibilities of WT and MT proteins at residues and secondary structure levels. To accomplish this, we measured root mean square fluctuation (RMSF) and secondary structures (α-helix, β-sheet, bridge, turn, coil, and bend) of WT and MT proteins for individual residue and domain throughout the simulation period. The fluctuation intensities varied among different regions of WT and all MTs proteins (Figure 2C). The regions occupied around positions I (160-190), II (250-260), III (345-355), and (IV) (420-430 residues) showed > 0.5nm of RMSFs in WT and all MTs. WT showed maximum fluctuation in region I (∼1 nm) and E234P showed maximum fluctuation at region IV (∼0.6 nm) (Figure 2C). The R322Q MT showed maximum fluctuation in region I (∼0.9nm). The R316Q, R322Q, and C392D MTs showed fewer fluctuations in regions II and IV thus implying that these MTs decreased the flexibilities of the given segments of the FTO protein. Further to investigate the effect of different mutations on the secondary structures of the FTO protein, a dictionary of secondary structure of protein (DSSP) was performed during MD simulation. It was found that the coil contents in all MTs were increased as compared to WT (Supplementary Table 5). The increased coil content led to decreased protein flexibility which was also observed during the RMSF analyses.

Further, we wanted to study the structural and functional properties of various MTs along with WT by calculating the solvent-accessible surface area (SASA), protein–protein (intra), and protein–water (inter) hydrogen bonding. Overall, SASA was further divided into hydrophobic and hydrophilic SASA which showed steady behaviors throughout the simulation period in WT as well as in all MTs (Supplementary Figure 4). The total mean values of SASA for WT, E234P, R316Q, R322Q, and C392D were approximately 258.7, 250.6, 254.4, 253.4, and 252.7 nm², respectively (Figure 3A) All MTs displayed lesser values of SASA as compared to WT. This suggested that all MTs have a low surface area of exposure for the interaction with other molecules. WT, E234P, R316Q, R322Q, and C392D exhibited average hydrophilic SASA values of about 126.4, 124.6, 125.7, 126.6, and 124.6 nm², while hydrophobic SASA of all MTs remained similar to WT. Furthermore, SASA per residues basis was calculated, and it was found that GLN86, LYS121, PHE176, TYR220, GLU281, and ARG459 exhibited higher SASA values than the rest of the amino acid residues in WT and all MTs (Figure 3B).

Wild Type and MTs showed consistent behaviors of intra- and inter-H-bonds formation (Supplementary Figure 5). The average intra-H-bonds of WT, E234P, and C392D displayed a similar number of H-bonds (360) while R316Q (354) and R322Q (348) MTs showed a lesser number of intra-H-bonds as compared to WT (Figure 3C). Lesser number of inter-H-bonds were found in E234P (1090) MTs as compared to WT (1092), while R316Q (1101) and R322Q (1100) MTs exhibited higher mean values of inter-H-bonds (Figure 3D). The above results showed that R316Q and R322Q MTs exhibited lesser intra-H-bonds indicating the above MTs are more flexible and less stable. Moreover, these MTs exhibited higher inter-H-bonds indicating that the above MTs had higher interacting capacities as compared to WT.

The function of protein depends upon its motion which is elucidated by conformational sampling in phase space from MD simulated trajectories. To study the appropriate motions of the protein resulting from atomic trajectories, we performed essential dynamics and reduce the dimension of data using principal component analysis (PCA). R322Q MT exhibited relatively higher trace values of covariance matrix while E234P, R316Q, and C392 MTs consisted of lower matrix values as compared to WT suggesting that the R322Q MT accompanied higher motions as compared to WT (Supplementary Table 6). During PCA, we examined the first 30 eigenvectors with corresponding eigenvalues which showed cumulative percentage of more than 85% in WT and all MTs (Figure 4A and Supplementary Table 6). In conjunction with the first 30 eigenvectors, we observed that the first 3 eigenvectors or PCs shared more than 50% of the cumulative percentage, which plays a significant role in the motions of protein (Supplementary Table 6). Moreover, to further confirm the movement of protein is not due to random diffusion, we investigated cosine value, in which a lower value (<0.2) was observed in the first 2 PCs (Supplementary Table 7). Further, to understand the protein conformation, we projected eigenvector 1 versus eigenvector 2 and eigenvector 2 versus eigenvector 3 in phase spaces (Figures 4B,C). In all 3 projections (PC1, PC2, and PC3), E234P occupied a smaller subspace in comparison to WT suggesting that it restricted the motions of the FTO protein (Figures 4B,C). Furthermore, to understand the motions at domain or 3D structural levels, we utilized 30 frames of the first eigenvector by sequentially superimposing the corresponding frames (Figure 5). We found that the motions were found at both N- and C-terminal domains. MTs such as E234P and R316Q (Figures 5B,C) have fewer motions in both domains, while R322Q and C392D MTs (Figures 5D,E) showed more motions as compared to WT (Figure 5A). Interestingly, mutations in the N-terminal domain (NTD) affected the motions of both NTD and the C-terminal domain (CTD). Motions were limited to loop and turn regions of N- and C-terminal domains (Figures 5a–e). The results of PCA are quite interesting as it could differentiate the properties of different MTs instead of having similar patterns of RMSD (Figures 6A,B and Supplementary Table 8), RMSF (Figures 6C,D), and Rg (Figures 6E,F) of NTD and CTD.

Betweenness centrality (C_B) analysis was performed to understand the key residues involved in the function of the protein. During this analysis, the residues in the protein were represented as nodes while the connections among these residues were denoted as edges. Here, we calculated 3 main centralities: betweenness centrality (C_B), closeness (C_C), and degree (C_D) for WT and all MTs out of which C_B has the significance of calculating the importance of each residue in the signaling of proteins upon different mutations. We mapped the residues having C_B ≥ 0.05 in WT and all MTs and found a similar number of residues involved in signaling (Figure 7A). Similarly, we have calculated the difference between the C_B values of WT and MTs and the residues satisfying C_B ≥ 0.02 (Figure 7B and Supplementary Table 9). These residues were mainly confined to the catalytic core and NTD and CTD of the FTO protein (Supplementary Figure 6). In this analysis, 28 and 40 residues play a significant role in the signaling of E234P and C392D MTs, respectively (Figure 7B and Supplementary Table 9), while no residues satisfying the conditions ≥ 0.02 existed in R316Q and R322Q MTs. This suggested that substitution of Arginine at 316 and 322 with Glutamine largely influenced the residue signaling of the FTO protein. Moreover, all individual MTs affected the dynamics and stabilities of different MTs (Supplementary Figure 7)

Furthermore, the residues-residues interactions network of WT and all MTs were plotted to investigate different types of bonding. To accomplish this task, a Cytoscape plugin with RINalyzer was employed in which nodes were represented as residue and different types of bonding such as van der Waals (VDW), hydrogen bonds (H-bond), pi-pi, and ionic interactions were represented as edges (Supplementary Figure 8). Further, a 2D interaction map of mutated residues was analyzed to examine the different types of interactions with surrounding residues and compared with native (Figure 7C). The van der Waals and hydrogen bonds were mainly affected in the E234P MT. The VDW interactions were reduced in all MTs except R316Q while R316Q displayed a similar number of VDW interactions with respect to WT (Figure 7C and Supplementary Figure 8). However, the H-bonds were decreased in R316Q and R322Q MTs as compared to WT while the H-bonds in C392D MT remained unaffected. The Pi-pi stack bonding was also reduced in R322Q MT. The residues network analysis results indicated that all non-covalent interactions were affected and impaired in almost all MTs.

FTO and its substrate binding were monitored by calculating the protein-ligand binding affinity using the docking approach which was measured in kcal/mol. N6-Methyl-deoxy-adenosine-5′-monophosphate (6mA) is one of the major substrates of FTO and we used this as a ligand for measuring the binding affinity upon different mutations. A combination of 5 docking tools was employed for WT- and MT-6mA complexes and it was found that all docking tools showed a similar pattern of docking scores even when all docking tools have different scoring functions (Figure 8A). The average docking scores of WT, E234P, R316Q, R322Q, and C392D were −6.2, −5.02, −5.1, −5.6, and −5.5 kcal/mol, respectively (Supplementary Table 10). The docking results indicated that all MTs exhibited high docking scores compared to WT suggesting that all MTs showed low binding affinities to the substrate than WT. Further to understand the residues participating in both hydrophobic and hydrophilic interactions during ligand binding, we analyzed the 3D and 2D plots of high-scored protein-ligand complexes (Figure 8B). The results suggested that 18 bonds (hydrophobic:15; hydrophilic:3) existed between the WT- and C392D-6mA complexes (Figure 8B). R322Q-6mA formed the maximum number of bonds (hydrophobic: 17; hydrophilic: 3) (Figure 8B) while E234P-6mA formed the minimum number of bonds (hydrophobic: 9; hydrophilic: 0) (Figure 8B). Further, different MTs led to a change in the electrostatic and hydrophobic nature of the core and binding regions of the protein compared with the WT-substrate complex (Supplementary Figure 9). The protein-ligand interaction studies indicated that binding affinities of almost all MT-ligand complexes were decreased which suggested that binding of the FTO substrate was largely affected by the given point mutations.

Stabilities and rescoring of binding energies of docking complexes (WT and MTs) were assessed through 100ns MD simulation (Figure 9). Binding energies of all docking complexes were calculated by g_mmpbsa script using the last 20ns-stabled MD simulated trajectories (55). Initially, the stabilities of WT- and MT-6mA complexes were assessed by monitoring the protein and drug RMSD behaviors followed by inspecting the Rg and RMSF of the protein in complexed form. RMSD of WT, R316Q, R322Q, and C392D showed steady behavior and stabilized after 40ns time while RMSD of E234P was stabilized after 70 ns simulation time (Figure 9A). RMSD of all WT and MT proteins in the complexed form displayed consistent and equilibrated behaviors. The stability of 6mA substrate in protein complexes was examined by measuring the ligand RMSD which showed steady behaviors after 45 ns simulation time in all cases (Figure 9B). Further, the compactness of protein in complexed form was monitored through the inspection of Rg behaviors and demonstrated that Rg of WT and all MT proteins displayed consistent and steady behaviors (Figure 9C). RMSF of protein in complexed form showed similar patterns as observed in the apo form (Supplementary Figure 10). The stability analyses of WT- and MT-6mA complexes suggested that all complexes were well stabilized throughout the simulation period. Binding energies of WT- and MT-6mA complexes were studied by calculating all energy components such as van der Waals (VDW), electrostatic, polar solvation, and solvent accessible surface area (SASA) energies. Binding energies of WT-, E234P-, R316Q-, R322Q-, and C392D-6mA were −180.911 ± 2.1, −121.83 ± 11, −32.92 ± 18.2, −83.78 ± 2.114.7, and −23.92 ± 18.5 kcal/mol, respectively (Figure 9D). The lowest binding free energy of the WT-6mA complex as compared to all MT-6mA complexes demonstrated that all MT-6mA complexes have low binding affinities than WT which also agreed with the docking results and VDW energy was the major contributor.