1 Overview

This project aims to improve the activity of Arabidopsis thaliana nicotinamide synthase (AtNAS1) through protein engineering strategies to overcome its C-terminal autoinhibition. The ultimate goal is to construct an engineered strain of Bacillus subtilis that can promote plant iron and zinc accumulation, thereby providing a new solution to the global problem of "hidden hunger." The experimental phase aims to generate high-performance AtNAS1 mutants through site-directed mutagenesis and functional validation. However, this process faces challenges such as a large number of mutation sites and unclear mechanisms, making traditional experimental screening labor-intensive and inefficient.

To overcome these bottlenecks, we introduced a multiscale computational biology model designed to provide accurate predictions for experimental studies. This computational framework comprises three logically connected submodules: First, molecular docking simulations are used to elucidate the binding modes of AtNAS1 mutants with their substrates and identify key interacting residues. This provides a precise theoretical basis for target selection for site-directed mutagenesis, guiding rational design and reducing experimental blindness. Secondly, through wet enzymatic reactions, the purified protein experimental data were quantitatively fitted, and the catalytic efficiency of the wild-type and mutants was accurately compared, such as the level of enzyme activity. This allowed for further comparison, screening the most promising mutants from a functional perspective and verifying the accuracy of the computational predictions. Finally, molecular dynamics simulations were used to dynamically demonstrate on a nanosecond to microsecond timescale how the mutations induce protein conformational changes, disrupt autoinhibitory interfaces, and affect substrate access. This approach elucidated the underlying mechanisms of enhanced or reduced mutant activity, providing a mechanistic supplement to the experimental results.

2 Introduction

Iron (Fe) and zinc (Zn), essential trace elements for the human body, play key roles in numerous life processes, including oxygen transport, DNA synthesis, immune function, and neural development[1]. However, globally, approximately two billion people are at risk of iron and zinc deficiency due to inadequate dietary intake, a phenomenon known as "hidden hunger" and particularly severe in developing regions where staple foods such as maize are the primary source[2]. The low bioavailability of iron and zinc in plant-based foods is partly due to the presence of anti-nutritional factors such as phytic acid, but the fundamental reason lies in the low absolute concentrations of these trace elements in the edible parts of crops[3]. Traditional nutritional interventions, such as supplementation and food fortification, are plagued by high costs, poor compliance, and low sustainability. Therefore, increasing the trace element content of staple crops through biofortification breeding is recognized as the most promising and sustainable solution[4].

Nicotinamide synthase (NAS) is a core enzyme in the plant iron and zinc metabolic pathway, catalyzing the synthesis of the non-protein amino acid nicotinamide (NA). NA is a key chelator of metal ions, such as trivalent iron and zinc. In plants, it is primarily responsible for long-distance transport of metal ions (e.g., translocation through the phloem to the grain) and regulation of intracellular homeostasis, making it crucial for the absorption, transport, and accumulation of metal elements[5]. The Arabidopsis thaliana AtNAS1 gene is a key member of this family[6]. However, recent studies have revealed an autoinhibitory domain at the C-terminus of the AtNAS1 protein, which inhibits its own enzymatic activity through intramolecular interactions. This explains why, in previous studies, simply overexpressing wild-type AtNAS1 had limited effect on increasing iron and zinc content in rice endosperm. Therefore, precise regulation of AtNAS1 enzymatic activity through protein engineering, rather than simply increasing its expression, has emerged as a new strategy to overcome this bottleneck and achieve the goal of iron and zinc biofortification.

This study aimed to rationally design and validate AtNAS1 mutants with enhanced enzymatic activity by integrating computational biology with experimental biology. To this end, we employed a multi-scale computational modeling approach, selected for the following reasons: First, molecular docking was used to predict the interactions between the AtNAS1 protein, its substrate (S-adenosylmethionine), and its autoinhibitory domain. This approach allows us to intuitively identify key interface residues where autoinhibition occurs and provides precise targets for site-directed mutagenesis aimed at deregulating autoinhibition, thus avoiding blind mutation screening. Second, enzymatic reaction bridges the gap between protein structure and function. We can quantify the effects of mutations on substrate affinity and catalytic efficiency, objectively selecting optimal mutants from a functional perspective. Finally, molecular dynamics (MD) simulations are crucial for revealing dynamic mechanisms. Enzyme function stems not from static structure but from the temporal motion of its atoms. MD simulations can reveal how autoinhibition is achieved through conformational changes, how mutations perturb the protein's dynamic network, and how the accessibility of the substrate binding channel is affected, thereby providing a mechanistic explanation for the docking and enzymatic experimental results.

3 Materials and Methods

3.1 Molecular Docking

Molecular Docking

3.2 Enzyme Reaction

Enzyme Reaction

3.3 Molecular Dynamics Simulation

Molecular Dynamics Simulation

For details, see 4.1, 4.2, and 4.3 on the Experiments page.

4. Results

4.1 Docking Results

Fig. 1 Two-dimensional, three-dimensional, and surface model images of the wild-type and mutant AtNAS1 proteins docked with SAM.

Two-dimensional and three-dimensional results show that wild-type AtNAS1 primarily forms multiple hydrogen bonds with SAM through residues such as GLU77, SER133, and LEU138, resulting in a relatively stable binding. The binding site of the truncated AtNAS1 shifts to residues such as GLN196 and ASP198, resulting in a more exposed binding pocket and reduced stability. Among the point mutants, T287A interacts with ASN285 and VAL284, but the binding is weak; T288A and G289A maintain moderate levels of interaction; and TRG→AAA forms multiple hydrogen bonds around residues such as GLY277, GLY278, and GLY286, resulting in a deeper and tighter binding pocket (Fig. 1).

The binding energy results were consistent with the structural observations: the wild-type binding energy was -6.0 kcal/mol, while the truncated forms decreased to -5.3 kcal/mol. T287A exhibited the weakest binding energy (-4.5 kcal/mol), while T288A and G289A exhibited -5.7 kcal/mol, and TRG→AAA exhibited the best binding energy (-6.5 kcal/mol). Overall, this suggests that the TRG→AAA mutation enhances SAM binding stability and may improve AtNAS1 enzymatic activity (Fig. 2).

Fig. 2 Binding energies between wild-type and mutant AtNAS1 proteins and SAM ligands.

4.2 Enzyme-catalyzed reaction

We assessed the enzymatic activity of nicotinamide synthase (NAS1) mutants in Functional Tests Part 1. The enzymatic activities of the wild-type and various mutants (Truncated, T287A, R288A, G289A, TRG_to_AAA) of NAS1 were evaluated by measuring the initial reaction rates (ΔA/Δt) and calculating the enzyme activities (nkat/mg) based on the parameters of the reaction system (Fig. 3B). Compared to the slope and correlation coefficient of AtNAS1-Wildtype, the enzyme activities of the different mutants varied significantly, with some mutants, such as AtNAS1-T287A and AtNAS1-R288A, exhibiting markedly reduced activities. In contrast, AtNAS1-Truncated and AtNAS1-TRG_to_AAA demonstrated steep slopes and high enzyme catalytic efficiency (Fig. 3A). The AtNAS1-TRG_to_AAA mutant protein, unlike AtNAS1-Truncated, retains the complete C-terminal region of the NAS1 protein, thereby avoiding potential functional redundancies associated with the truncated form. Notably, the TRG_to_AAA mutant also exhibited a relatively higher level of nicotinamide synthesis activity compared to the truncated protein (Fig. 3B). Therefore, the TRG_to_AAA mutant may possess greater potential for application than the truncated variant.

Fig. 3 Kinetic and Thermodynamic Analysis of Engineered Variants of AtNAS1 protein. A. Enzyme Kinetics of AtNAS1 mutated protein; B. Enzyme Activity Analysis.

4.3 MD Simulation Results

Wildtype:

Fig. 4A:Radius of Gyration; B:RMSF (Root Mean Square Fluctuation); C:Rg (Radius of Gyration); D:SASA (Solvent Accessible Surface Area) ; E:Hbond (Hydrogen Bonds)

After 20 ns, the RMSD values stabilize, indicating that the protein has reached a stable state after ligand binding. The ligand does not move significantly within the binding site, indicating relatively stable binding. The RMSF values for most residues are low, indicating a relatively stable protein structure. The RMSF at the C-terminus is high, possibly representing a flexible region of the protein. After 10 ns, the Rg value decreases and stabilizes, indicating that the protein structure becomes more compact after ligand binding. After 10 ns, the SASA decreases, indicating that the protein wraps more tightly around the ligand, further confirming stable binding. After 10 ns, the number of hydrogen bonds stabilizes, indicating that the protein and ligand form stable hydrogen bonding interactions (Fig 4A-E).

Fig. 5 A&B:Gibbs (Gibbs Energy Landscaping) Gibbs (Gibbs Energy Landscaping); C:Average binding free energy D:Residue decomposition energy

Energy troughs (dark blue areas) correspond to high-probability conformational clusters and represent stable states (e.g., protein-ligand binding states).Energy plateaus or saddle points (yellow/red areas) represent transition states or energy barriers, separating different conformational clusters and reflecting the kinetic bottlenecks of conformational transitions (Fig5A-B).

The total binding free energy is -36.83 kcal/mol, indicating that the ligand and protein can spontaneously bind and that the complex is thermodynamically stable. This provides an energetic theoretical basis for further studies of the biological functions of this ligand-protein complex (Fig 5C). The energy values corresponding to protein residues such as THR-139, SER-133, and ARG-233 are relatively negative, indicating that these residues contribute relatively strongly to the stabilization of the protein-ligand interaction(Fig 5D).

AtNAS1_TRG_to_AAA:

After 10 ns, the RMSD values stabilized. Although there were some fluctuations, they were generally small. This indicates that the protein reached a stable state after ligand binding, and the ligand did not move significantly within the binding site, indicating relatively stable binding.

The RMSF values for most residues were low, indicating a relatively stable protein structure. The RMSF at the C-terminus was high, possibly representing a flexible region of the protein. After 10 ns, the Rg value decreased and stabilized, indicating that the protein structure became more compact after ligand binding. After 10 ns, the SASA decreased, indicating that the protein wrapped more tightly around the ligand, further confirming stable binding. After 10 ns, the number of hydrogen bonds stabilized, indicating that the protein and ligand formed stable hydrogen bonding interactions (Fig. 6A-E).

Fig. 6A:Radius of Gyration; B:RMSF (Root Mean Square Fluctuation); C:Rg (Radius of Gyration); D:SASA (Solvent Accessible Surface Area) ; E:Hbond (Hydrogen Bonds)

Energy troughs (dark blue areas) correspond to high-probability conformational clusters and represent stable states (e.g., protein-ligand binding states). Energy plateaus or saddle points (yellow/red areas) represent transition states or energy barriers, separating different conformational clusters and reflecting the kinetic bottlenecks of conformational transitions (Fig 7A-B).

The total binding free energy is -17.06 kcal/mol, indicating that the ligand and protein can spontaneously bind and that the complex is thermodynamically stable. This provides an energetic basis for further studies of the biological functions of this ligand-protein complex (Fig 7C).

The energy values corresponding to protein residues such as LEU-276, PHE-292, and GLY-277 are relatively negative, indicating that these residues contribute relatively strongly to the stabilization of the protein-ligand interaction (Fig 7D).

Fig. 7 A&B:Gibbs (Gibbs Energy Landscaping) Gibbs (Gibbs Energy Landscaping); C:Average binding free energy D:Residue decomposition energy

AtNAS1_Truncated:

Fig. 8A:Radius of Gyration; B:RMSF (Root Mean Square Fluctuation); C:Rg (Radius of Gyration); D:SASA (Solvent Accessible Surface Area) ; E:Hbond (Hydrogen Bonds)

After 20 ns, the RMSD value stabilizes. After approximately 40 ns, fluctuations begin, but the overall fluctuations are relatively small, indicating stability. The RMSF values for most residues are low, indicating a relatively stable protein structure. The RMSF at the C-terminus is high, likely due to a flexible region of the protein. After 10 ns, the Rg value decreases and stabilizes, indicating a more compact protein structure after ligand binding. After 10 ns, the SASA decreases, indicating a tighter wrapping of the protein around the ligand, further confirming stable binding (Fig 8A-E).

Fig. 9 A&B:Gibbs (Gibbs Energy Landscaping) Gibbs (Gibbs Energy Landscaping); C:Average binding free energy D:Residue decomposition energy

Energy troughs (dark blue areas) correspond to high-probability conformational clusters and represent stable states (e.g., protein-ligand binding states). Energy plateaus or saddle points (yellow/red areas) represent transition states or energy barriers, separating different conformational clusters and reflecting the kinetic bottlenecks of conformational transitions (Fig 9A-B)..

The total binding free energy is -13.63 kcal/mol, indicating that the ligand and protein can spontaneously bind and that the complex is thermodynamically stable. This provides an energetic theoretical basis for further studies of the biological functions of this ligand-protein complex (Fig 9C).

The energy values corresponding to protein residues such as ASP-198, HIS-120, and GLY-227 are relatively negative, indicating that these residues contribute relatively strongly to the stabilization of the protein-ligand interaction (Fig 9D).

5. Discussion

Through a multi-level integration of molecular docking, molecular dynamics simulations, and enzyme kinetics experiments, this study successfully revealed the interaction mechanism between the AtNAS1 protein and its mutants and the SAM ligand, and verified the potential application value of the TRG→AAA mutant. Molecular docking results showed that wild-type AtNAS1 forms multiple hydrogen bonds with SAM through residues such as GLU77, SER133, and LEU138, resulting in a relatively stable binding (binding energy -6.0 kcal/mol), while the binding site of the truncated mutant shifts to residues such as GLN196 and ASP198, resulting in a more exposed binding pocket and reduced stability (binding energy -5.3 kcal/mol). The point mutant T287A binds the weakest (-4.5 kcal/mol), T288A and G289A bind moderately (-5.7 kcal/mol), and the TRG→AAA mutant forms a deeper and tighter binding pocket through residues such as GLY277, GLY278 and GLY286, exhibiting optimal binding energy (-6.5 kcal/mol). This prediction is consistent with the results of enzyme kinetic experiments: both the wild type and TRG→AAA show high catalytic efficiency (steep slope and high nkat/mg value), while the activity of T287A and R288A is significantly reduced. Molecular dynamics simulations further confirmed binding stability: the wild-type showed stable RMSD, Rg, SASA, and hydrogen bond counts after 10 ns, indicating a compact and stable structure after ligand binding (binding free energy -17.06 kcal/mol), with key residues LEU-276, PHE-292, and GLY-277 contributing significantly. The TRG→AAA mutant showed stable RMSD after 20 ns, and although there were fluctuations after 40 ns, the overall structure was compact (binding free energy -13.63 kcal/mol), with residues ASP-198, HIS-120, and GLY-227 playing a major role. The truncated variant showed high C-terminal flexibility and lower binding stability.

These results highlight the complementary advantages of multiple approaches: molecular docking provides initial binding conformations and energy predictions, molecular dynamics simulations reveal the dynamic stability of the binding process and the energetic contributions of key residues, and enzyme kinetics experiments directly verify the functional output. This study not only elucidates the catalytic mechanism of AtNAS1 but also provides clear guidance for subsequent experimental design: first, the TRG→AAA mutant should be selected as a priority candidate for in vivo functional validation (e.g., plant validation experiments) due to its high binding stability and enzymatic activity; second, key residues identified by molecular dynamics (such as SER-133 and ARG-233) can serve as targets for further rational design; finally, the structural instability of truncation mutants suggests that the C-terminal domain is crucial for the overall protein conformation, and subsequent studies should avoid blind truncation. This work demonstrates the potential of combining computational and experimental approaches in enzyme engineering and provides the iGEM team with a standardized research paradigm from prediction to validation.

6. Limitations

While this study provides comprehensive insights by integrating computational and experimental approaches, several limitations must be acknowledged. First, due to the lack of experimental crystal structures of AtNAS1 and its mutants, structure predictions relied heavily on homology modeling and molecular docking simulations. While these methods are widely used, they can introduce errors in side-chain packing, loop modeling, and ligand pose estimation, potentially affecting the accuracy of binding site characterization. Second, key kinetic parameters (such as SAM binding constants and catalytic rates) were partially derived from literature values ​​or computational predictions rather than direct experimental measurements, which can lead to biased quantitative energy assessments. Third, the timescale of molecular dynamics simulations is limited to 100 nanoseconds, which may not capture rare but critical conformational transitions occurring outside this window. Furthermore, simulations were performed under ideal conditions that may not fully reflect the physiological environment.

7. Future Work

Based on these findings, we propose several directions for further research. First, computational parameters should be optimized by extending molecular dynamics simulations to the microsecond timescale to capture rare conformational changes and improve sampling accuracy. This should be accompanied by incorporating well-defined solvent models and physiological ion concentrations to better simulate cellular states. Second, further integration with wet experiments is crucial: determining the crystal structures of key mutants (e.g., TRG→AAA) by X-ray crystallography to validate the predicted binding sites; measuring binding affinities using isothermal titration calorimetry (ITC) or surface plasmon resonance (SPR) experiments to validate calculated binding energies; and conducting in vivo experiments in plant models to assess the metal accumulation phenotype and its physiological relevance.

8. References

1. Prashanth L, Kattapagari K K, Chitturi R T, et al. A review on role of essential trace elements in health and disease[J]. Journal of Dr. YSR University of Health Sciences, 2015, 4(2): 75-85.

2. Al Mamun M A, Ghani R B A. The role of iron and zinc in cognitive development of children[J]. Asian Journal of Medical and Biological Research, 2017, 3(2): 145-151.

3. Durrani A M, Parveen H. Zinc deficiency and its consequences during pregnancy[J]. Microbial biofertilizers and micronutrient availability: The role of zinc in agriculture and human health, 2021: 69-82.

4. Galani Y J H, Orfila C, Gong Y Y. A review of micronutrient deficiencies and analysis of maize contribution to nutrient requirements of women and children in Eastern and Southern Africa[J]. Critical Reviews in Food Science and Nutrition, 2022, 62(6): 1568-1591.

5. Nguyen P, Grajeda R, Melgar P, et al. Effect of zinc on efficacy of iron supplementation in improving iron and zinc status in women[J]. Journal of Nutrition and Metabolism, 2012, 2012(1): 216179.

6. Zhu C, Sanahuja G, Yuan D, et al. Biofortification of plants with altered antioxidant content and composition: genetic engineering strategies[J]. Plant biotechnology journal, 2013, 11(2): 129-141.