Overview

To effectively regulate the TvLac genes, we sought a transcriptional regulator that could facilitate laccase production in response to the presence of EE2 in the environment. Our research led us to identify a Tet-R-like protein, as described by Juan Ibero et al. (2021) [1]. This protein is known for its ability to regulate gene transcription in response to a compound similar to EE2, specifically estradiol (E2).

Tet-R proteins are characterized as dimeric entities composed of three subunits: one for DNA binding, another for dimerization, and a third for ligand reception. Upon ligand binding, a conformational change occurs that decreases the affinity of the DNA-binding subunit for DNA, thereby enabling RNA polymerase to transcribe downstream genes.

However, we encountered a limitation: this Tet-R-like protein was ineffective in regulating gene transcription in the presence of EE2, and its specific DNA-binding sequence was unknown.

To address this challenge, we employed protein engineering software to design a novel protein that could regulate gene expression in the presence of EE2, complete with a defined DNA-binding site.

Protein modelation

To kick off the protein engineering process, we first generated a 3D model of the protein using AlphaFold 3 [2], which is known for accurately predicting the PDB structures of dimeric proteins. To enhance the reliability of our model, we also created comparative models through Rosetta Comparative Modeling [3], leveraging five analogous structures obtained via NCBI Blastp: 2fq4 [4], 2id3 [5], 2zb9 [6], 3qbm [7], and 4jkz [8]. Following the generation of these models, we utilized CaviDB [9] to predict the allosteric site, successfully identifying a cavity on the protein's surface. Our final decision to adopt the AlphaFold 3 model was influenced by its higher confidence score of 0.9, along with its exceptional capacity to maintain the integrity of the cavity at the allosteric site, surpassing that of the Rosetta models. This preservation was crucial in our selection process. Additionally, we assessed the quality of the PDB structure using MolProbity, which indicated an extremely high quality for our protein model.

We compared different softwares for protein modeling by the cavity preservation and functional comparation of models

Protein cavity identification using the CaviDB server

Ligand docking

For the molecular dockings, we used AutoDock Vina [10] in conjunction with Chimera [11] to position the ligand within the cavity. After completing the initial docking, we refined the results using Rosetta Ligand-Docking [12]. The XML script we used for this refinement was inspired by Lemmon et al. (2011) [13] for flexible ligand docking with proteins. Below is the final XML script used for Rosetta Ligand-Docking:


<ROSETTASCRIPTS>
  <SCOREFXNS>
      <ScoreFunction name="ligand_soft_rep" weights="ligand_soft_rep"/>
      <ScoreFunction name="hard_rep" weights="ligand"/>
  </SCOREFXNS>

  <LIGAND_AREAS>
      <LigandArea name="inhibitor_dock_sc" chain="X" cutoff="6.0" add_nbr_radius="true" all_atom_mode="true" minimize_ligand="10"/>
      <LigandArea name="inhibitor_final_sc" chain="X" cutoff="6.0" add_nbr_radius="true" all_atom_mode="true"/>
      <LigandArea name="inhibitor_final_bb" chain="X" cutoff="7.0" add_nbr_radius="false" all_atom_mode="true" Calpha_restraints="0.3"/>
  </LIGAND_AREAS>

  <INTERFACE_BUILDERS>
      <InterfaceBuilder name="side_chain_for_docking" ligand_areas="inhibitor_dock_sc"/>
      <InterfaceBuilder name="side_chain_for_final" ligand_areas="inhibitor_final_sc"/>
      <InterfaceBuilder name="backbone" ligand_areas="inhibitor_final_bb" extension_window="3"/>
  </INTERFACE_BUILDERS>

  <MOVEMAP_BUILDERS>
      <MoveMapBuilder name="docking" sc_interface="side_chain_for_docking" minimize_water="true"/>
      <MoveMapBuilder name="final" sc_interface="side_chain_for_final" bb_interface="backbone" minimize_water="true"/>
  </MOVEMAP_BUILDERS>

  <SCORINGGRIDS>
      <ClassicGrid grid_name="classic" weight="1.0"/>
  </SCORINGGRIDS>

  <MOVERS>
      <Transform name="transform" chain="X" box_size="10.0" move_distance="0.2" angle="20" cycles="500" repeats="1" temperature="5"/>
      <SlideTogether name="slide_together" chains="X"/>
      <HighResDocker name="high_res_docker" cycles="6" repack_every_Nth="3" scorefxn="ligand_soft_rep" movemap_builder="docking"/>
      <FinalMinimizer name="final" scorefxn="hard_rep" movemap_builder="final"/>
      <InterfaceScoreCalculator name="add_scores" chains="X" scorefxn="hard_rep"/>
  </MOVERS>

  <PROTOCOLS>
      <Add mover_name="transform"/>
      <Add mover_name="slide_together"/>
      <Add mover_name="high_res_docker"/>
      <Add mover_name="final"/>
      <Add mover_name="add_scores"/>
  </PROTOCOLS>
</ROSETTASCRIPTS>

This Rosetta protocol begins by first transforming the ligand to explore different possible conformations. It then brings the ligand and protein closer together through a low-resolution docking phase, during which only the backbone of both the protein and ligand is moved. Once the ligand is in proximity to the protein, a high-resolution docking step is performed, every atom involved in the interaction to achieve a more precise fit. The process concludes with energy minimization of the entire protein-ligand complex, generating a score to assess the stability of the docking models. A lower score indicates a more physically and chemically stable complex.

After the docking procedure, we used Discovery Studio [14] to generate 2D interaction maps, providing a detailed view of how the protein and ligand interact. Our analysis showed that E2 successfully entered the cavity of the dimeric protein, obtaining a docking score of -655, while EE2 was unable to enter the cavity and received a slightly higher score of -635. The interaction between the ligand and the protein was predominantly mediated by hydrophobic interactions and hydrogen bonds.

Protein engineering

To initiate the mutations, we first focused on altering the DNA-binding subunit to enable our protein to bind to a known DNA sequence. Specifically, we replaced the FASTA sequence of our protein's DNA-binding domain with that of the TetR protein base on the PDB 1QPI [15], based on the engineering described by Dimas et al. (2019) [16]. After modifying the DNA-binding domain, we generated a PDB file using Alphafold 3, which predicted the DNA interaction with a confidence level of 0.9. This prediction utilized the DNA sequence to model the protein-DNA interaction.

For the mutation process, we followed the steps outlined in the accompanying workflow. Initially, we generated a 3D structural model of the protein with the desired mutation using AlphaFold 3, which provided us with a PDB file. This model served as the basis for subsequent docking studies. Using AutoDock Vina, integrated with Chimera, we introduced the ligand into the cavity of the protein structure altered by the mutation. This docking aimed to explore how the mutation affected the binding site and ligand interaction.

After initial docking, we refined the ligand-protein interaction using the Rosetta Ligand Docking protocol. This step allowed us to perform both low-resolution docking to position the ligand and high-resolution docking to adjust atomic interactions and optimize the fit. Energy minimization was applied, and a docking score was calculated for each model, where lower scores indicate greater physical and chemical stability of the complex.

We then analyzed the molecular interactions between the mutated protein and the ligand using Discovery Studio. This step provided detailed 2D interaction maps, allowing us to visualize the specific interactions, such as hydrogen bonds, hydrophobic interactions, and any changes induced by the mutation.

Initially, the native protein produced a docking score of -1240 with E2, indicating a stable interaction. In contrast, EE2 did not enter the binding cavity, resulting in no meaningful docking score. However, after completing the mutation cycle and re-evaluating the docking, we observed improved binding affinity. The mutated protein obtained a docking score of -1260 with E2, showing enhanced stability compared to the native protein. Additionally, the EE2 ligand, which previously did not enter the cavity, now docked successfully with a score of -1265.

This improvement in both ligand interactions was reflected in the subsequent interaction maps, highlighting a stronger and more favorable binding profile between the mutated protein and the ligands. The cycle was then repeated by introducing further mutations into the FASTA sequence of the protein, continuing the process of optimization.

In silico validation

For the molecular dynamics simulation, GROMACS [17] was used in conjunction with the CHARMM36m force field [18], specifically designed for protein simulations. The TIP3P water model, optimized for this force field, was employed, and an octahedral solvation box was utilized. Na and Cl ions were added to neutralize the system’s charge using the Monte Carlo method.

Energy minimization (EM), canonical ensemble (NVT), and isothermal-isobaric ensemble (NPT) were performed to equilibrate the system. The production phase was then conducted under NPT conditions for 150 nanoseconds, with the initial nanoseconds considered as a stabilization phase. The system was set up in the CHARMM-GUI server [19], where the necessary files for ligand interaction with the force field were generated by introducing the ligand to the cavities of the protein. The entire simulation ran for 150 nanoseconds performing the analysis via the build-in tools from GROMACS and the gmx_MMPBSA software for the calculation of free energy binding for the last 20 ns of the simulation [20].

We conducted four simulations under specific conditions, first we simulated the protein-DNA interactions without the ligand and then the protein-EE2-DNA interaction to simulate the conformational change induced by the ligand and its effect on the protein.

On the graph we can see a clear difference in the interaction between the protein and the DNA when the ligand is present in the complex, as seen in the graphs D and E the number of hydrogen bonds is reduced when the ligand is found in the cavity as compared to when it’s not, this is represented by the loss of affinity seen in the binding free energy calculated for the protein and the DNA, this indicates a conformational change as seen in the radius of gyration and the RMSF of the protein being significantly different when the EE2 is in the complex. This suggests a loss of affinity for DNA, indicating that the presence of EE2 in the allosteric site negatively affects the DNA binding site.

References

[1] Ibero, J., Galán, B., & García, J. L. (2021). Identification of the EdcR Estrogen-Dependent Repressor in Caenibius tardaugens NBRC 16725: Construction of a Cellular Estradiol Biosensor. Genes, 12(12), 1846. https://doi.org/10.3390/genes12121846

[2] Abramson, J., Adler, J., Dunger, J., Evans, R., Green, T., Pritzel, A., Ronneberger, O., Willmore, L., Ballard, A. J., Bambrick, J., Bodenstein, S. W., Evans, D. A., Hung, C., O’Neill, M., Reiman, D., Tunyasuvunakool, K., Wu, Z., Žemgulytė, A., Arvaniti, E., . . . Jumper, J. M. (2024). Accurate structure prediction of biomolecular interactions with AlphaFold 3. Nature, 630(8016), 493-500. https://doi.org/10.1038/s41586-024-07487-w

[3] Song, Y., DiMaio, F., Wang, R. Y. R., Kim, D., Miles, C., Brunette, T. J., ... & Baker, D. (2013). High-resolution comparative modeling with RosettaCM. Structure, 21(10), 1735-1742.

[4] Zhang, R., Wu, R., Moy, S., Cymborowski, M., Minor, W., Joachimiak, A., & Midwest Center for Structural Genomics (MCSG). (2006). The crystal structure of the transcriptional regulator (TetR family) from Bacillus cereus (PDB ID: 2FQ4). Protein Data Bank. https://doi.org/10.2210/pdb2fq4/pdb

[5] Grabowski, M., Chruszcz, M., Koclega, K.D., Cymborowski, M., Gu, J., Xu, X., Savchenko, A., Edwards, A., Joachimiak, A., Minor, W., & Midwest Center for Structural Genomics (MCSG). (2006). Crystal structure of transcriptional regulator SCO5951 from Streptomyces coelicolor A3(2) (PDB ID: 2ID3). Protein Data Bank. https://doi.org/10.2210/pdb2id3/pdb

[6] Okada, U., Kondo, K., Watanabe, N., Yao, M., & Tanaka, I. (2008). Crystal structure of TetR family transcription regulator SCO0332 (PDB ID: 2ZB9). Protein Data Bank. https://doi.org/10.2210/pdb2zb9/pdb

[7] Joint Center for Structural Genomics (JCSG). (2011). Crystal structure of a TetR transcriptional regulator (Caur_2221) from Chloroflexus aurantiacus J-10-FL at 1.80 Å resolution (PDB ID: 3QBM). Protein Data Bank. https://doi.org/10.2210/pdb3qbm/pdb

[8] Yang, S.F., Gao, Z.Q., He, Z.G., & Dong, Y.H. (2013). Crystal structure of ms6564 from Mycobacterium smegmatis (PDB ID: 4JKZ). Protein Data Bank. https://doi.org/10.2210/pdb4jkz/pdb

[9] Rueda, A. J. V., Bulgarelli, F. L., Palopoli, N., & Parisi, G. (2023). CaviDB: a database of cavities and their features in the structural and conformational space of proteins. Database, 2023. https://doi.org/10.1093/database/baad010

[10] Trott, O., & Olson, A. J. (2009). AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. Journal Of Computational Chemistry, 31(2), 455-461. https://doi.org/10.1002/jcc.21334

[11] Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M., Meng, E. C., & Ferrin, T. E. (2004). UCSF Chimera—A visualization system for exploratory research and analysis. Journal Of Computational Chemistry, 25(13), 1605-1612. https://doi.org/10.1002/jcc.20084

[12] Meiler, J., & Baker, D. (2006). ROSETTALIGAND: Protein–small molecule docking with full side‐chain flexibility. Proteins Structure Function And Bioinformatics, 65(3), 538-548. https://doi.org/10.1002/prot.21086

[13] Lemmon, G., & Meiler, J. (2011). Rosetta Ligand Docking with Flexible XML Protocols. Methods In Molecular Biology, 143-155. https://doi.org/10.1007/978-1-61779-465-0_10

[14] Jejurikar, B. L., & Rohane, S. H. (2021). Drug Designing in Discovery Studio. Asian Journal Of Research In Chemistry, 14(2), 135-138. https://doi.org/10.5958/0974-4150.2021.00025.0

[15] Saenger, W., Hinrichs, W., Orth, P., Schnappinger, D., & Hillen, W. (2000). Structural basis of gene regulation by the tetracycline inducible Tet repressor-operator system. Nature Structural & Molecular Biology, 7(3), 215-219. https://doi.org/10.1038/73324

[16] Dimas, R. P., Jordan, B. R., Jiang, X., Martini, C., Glavy, J. S., Patterson, D. P., Morcos, F., & Chan, C. T. Y. (2019b). Engineering DNA recognition and allosteric response properties of TetR family proteins by using a module-swapping strategy. Nucleic Acids Research, 47(16), 8913-8925. https://doi.org/10.1093/nar/gkz666

[17] Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A. E., & Berendsen, H. J. C. (2005b). GROMACS: Fast, flexible, and free. Journal Of Computational Chemistry, 26(16), 1701-1718. https://doi.org/10.1002/jcc.20291

[18] Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., De Groot, B. L., Grubmüller, H., & MacKerell, A. D. (2016). CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nature Methods, 14(1), 71-73. https://doi.org/10.1038/nmeth.4067

[19] Jo, S., Kim, T., Iyer, V. G., & Im, W. (2008). CHARMM‐GUI: A web‐based graphical user interface for CHARMM. Journal Of Computational Chemistry, 29(11), 1859-1865. https://doi.org/10.1002/jcc.20945

[20] Valdés-Tresanco, M.S., Valdés-Tresanco, M.E., Valiente, P.A. and Moreno E. gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS. Journal of Chemical Theory and Computation, 2021 17 (10), 6281-6291. https://pubs.acs.org/doi/10.1021/acs.jctc.1c00645

To address this challenge, we employed protein engineering software to design a novel protein that could regulate gene expression in the presence of EE2, complete with a defined DNA-binding site.

We compared different softwares for protein modeling by the cavity preservation and functional comparation of models

Rosseta script