Using yeast two-hybrid system and molecular dynamics simulation to detect venom protein-protein interactions.

Proteins and peptides are major components of snake venom. Venom protein transcriptomes and proteomes of many snake species have been reported; however, snake venom complexity (i.e., the venom protein-protein interactions, PPIs) remains largely unknown. To detect the venom protein interactions, we used the most common snake venom component, phospholipase A2s (PLA2s) as a "bait" to identify the interactions between PLA2s and 14 of the most common proteins in Western diamondback rattlesnake (Crotalus atrox) venom by using yeast two-hybrid (Y2H) analysis, a technique used to detect PPIs. As a result, we identified PLA2s interacting with themselves, and lysing-49 PLA2 (Lys49 PLA2) interacting with venom cysteine-rich secretory protein (CRISP). To reveal the complex structure of Lys49 PLA2-CRISP interaction at the structural level, we first built the three-dimensional (3D) structures of Lys49 PLA2 and CRISP by a widely used computational program-MODELLER. The binding modes of Lys49 PLA2-CRISP interaction were then predicted through three different docking programs including ClusPro, ZDOCK and HADDOCK. Furthermore, the most likely complex structure of Lys49 PLA2-CRISP was inferred by molecular dynamic (MD) simulations with GROMACS software. The techniques used and results obtained from this study strengthen the understanding of snake venom protein interactions and pave the way for the study of animal venom complexity.

Introduction

Snakebite remains a major public health problem (Kasturiratne et al., 2008), particularly in impoverished rural communities (Gutiérrez et al., 2017). The World Health Organization therefore recognized snakebite as a high priority neglected tropical disease (WHO, 2018). The venom proteome and transcriptome of many snake species have been reported, such as European viper (Leonardi et al., 2019) and Sidewinder Rattlesnakes (Hofmann et al., 2018). Recently, Suryamohan et al. (2020) completed transcriptome and genome sequencing of Indian cobra. Moreover, the pathophysiological effects of many individual snake venom components have been characterized. However, snake venom is a complex mixture of proteins and peptides that are stored in the gland lumen to exert a wide range of toxic actions during envenomation; therefore, we hypothesize that venom proteins interact with each other to make a cocktail of proteins and peptides and synergistically exert their toxic effects. To our knowledge, in addition to the extensively reported lysing-49 PLA2s (Lys49 PLA2) homodimers (Almeida et al., 2016, Salvador et al., 2019), aspartic acid-49 PLA2(Asp49) homodimers (Corrêa et al., 2008) and Lys49-Asp49 heterodimer (Mora-Obando et al., 2014), little is known about the venom protein-protein interactions (PPIs), particularly the non-homologous venom protein interactions. Understanding the snake venom PPIs is critically important for deciphering the synergistic actions which exert various pathophysiological effects of snakebites.

Western diamondback rattlesnake (Crotalus atrox) is likely responsible for most snakebite fatalities in northern Mexico, and the second greatest number in the U.S.A. after Eastern diamondback rattlesnake (Crotalus adamanteus) (Campbell and Lamar, 2004). The proteome (Calvete et al., 2009) and incomplete transcriptome (Rokyta et al., 2011, Jia et al., 2020) of C. atrox venom protein have been reported; however, its venom PPIs remain poorly understood. The lack of a deep understanding of how C. atrox venom proteins interact is an important roadblock in advancing efforts to uncover the corresponding mechanisms of the synergistic effects for developing antivenom to treat snakebites. In our previous work (Jia et al., 2020), we identified 14 protein transcripts from the venom of C. atrox by using reverse transcription polymerase chain reaction (RT-PCR) and modeled their three-dimensional (3D) structures using predicted amino acid sequences and computational approach (MODELLER). These 14 transcripts are particularly important for us to further understand the complexity of C. atrox venom as well as to uncover the toxicity of each venom component by producing recombinant venom proteins. Therefore, to advance our knowledge about the venom complexity of C. atrox, in this study we attempted to detect the venom PPIs by using these 14 venom protein transcripts including PLA2s and cysteine-rich secretory protein (CRISP) transcripts. Snake venom PLA2, the most extensively studied venom components, is a major protein in most snake venoms. Snake venom PLA2s were classified into Group I (Elapidae) or Group II (Viperidae) of PLA2 superfamily that is composed of 16 Groups (Dennis et al., 2011). Group II snake venom PLA2s were further subdivided into at least two subgroups: the catalytically active Asp49-PLA2s which have an aspartic acid residue at the position 49, and catalytically inactive lys49-PLA2 homologue (or PLA2-like myotoxins) that the aspartic acid residue at the position 49 is replaced by Lysine (Renetseder et al., 1985, Maraganore and Heinrikson, 1986). In addition to the primary catalytic function, snake venom PLA2s display an array of pathophysiological effects including neurotoxicity, myotoxicity, cardiotoxicity, hemolytic activity, hypotensive activity, anticoagulant activity, etc. (Reviewed in Kini (2003)). CRISPs, single chain polypeptides with molecular weights of ~20–30 kDa, were found in mammalian epididymis and the immune system, and then in many snake venoms (Yamazaki and Morita, 2004). The overall structure of CRISP family proteins is well-conserved, consisting of an N-terminal domain (PR-1), and a conserved divalent metal-ion binding site that connected PR-1 to C-terminal domain (CR) (Shiol et al., 2019). The physiological functions of mammalian CRISPs, such as CRISP-1, -2 and -3, are associated with reproduction, cancer, and immune responses, while snake venom CRISPs inhibit ion channels (reviewed in Tadokoro et al. (2020)).

There are many molecular biology approaches that can be used to detect venom PPIs (Ho et al., 2002, Jares-Erijman and Jovin, 2006, Michnick et al., 2010). Among them, yeast two-hybrid (Y2H) analysis, since its inception (Fields and Song, 1989), has been extensively used to detect the physical interactions between any two proteins in different species. More importantly, Y2H analysis is amenable to being used in a high throughput setting, allowing a protein of interest to be tested for interactions with many proteins (Galletta and Rusan, 2015). For example, PPIs in yeast (Yu et al., 2008), nematode (Simonis et al., 2009), plant (Trigg et al., 2017, Marshall et al., 2019), human (Rolland et al., 2014) and bacterium (Rajagopala et al., 2014) were identified by high throughput Y2H analysis. By using Y2H analysis coupled with various computational programs including template-based modeling – MODELLER, protein-protein docking and molecular dynamics (MD) simulations, here we report the most common snake venom protein – phospholipase A2s (PLA2s) interacting with themselves and Lys49 PLA2 interacting with venom CRISP, as well as the most likely conformation of Lys49 PLA2-CRISP.

Materials and methods

Venom protein transcripts

Fourteen of the most common venom protein transcripts including Asp49 PLA2, Lys49 PLA2, snake venom metalloproteinase I (SVMPI), snake venom metalloproteinase II (SVMPII), snake venom metalloproteinase III (SVMPIII), snake venom serine proteinase (SVSP), L-amino acid oxidase (LAAO), C-type lectin, crotamine, three-finger toxin (3FTx), vespryn, CRISP, epidermal growth factor-like domain protein (EGF) and phospholipase B (PLB) used in this study, were obtained from our previous results (Jia et al., 2020). All translated amino acid sequences of these protein transcripts are available in our previous publication (Jia et al., 2020). All chemicals were purchased from Sigma-Aldrich (St. Louis, MO, USA).

Yeast two-hybrid analysis

Y2H activation and binding domain vectors (pAD and pBD) were purchased from Takara Bio USA, Inc (Mountain View, CA, USA). Polymerase Chain reactions (PCR) were carried out using cloned transcripts (Jia et al., 2020) as templates and Phusion High-Fidelity DNA Polymerase purchased from ThermoFisher Scientific (Waltham, MA, USA) by standard molecular biology methods (Table 1. Primer sequences). The PCR products of 14 mature transcripts were individually purified from agarose gel. Each transcript was ligated into EcoRI and BamHI sites of pAD vector as “prey” by using either restriction enzyme based-method or homologous recombination-based approach if enzymatic sites were present in the transcript sequence (such as EGF and SVSP). Using the same enzymatic method, transcripts for Lys49 PLA2, Asp49 PLA2 and CRISP were cloned into EcoRI and BamHI sites of pBD vector as “bait”. Empty vectors (pAD, pBD) served as negative controls. The conjunctions of recombinant DNAs were confirmed by Sanger sequencing, and the confirmed “bait” and “prey” constructs were co-transformed into yeast (Saccharomyces cerevisiae) strain (Y2HGold) cells (Takara Bio USA, Inc) based on the Takara instruction manual. The co-transformed cells were subsequently spread on double dropout (DDO, without tryptophan and leucine) medium and incubated at 30 °C until the colonies grew (~2 days) for confirming the successful pairwise transformation. Plasmid DNAs were extracted from co-transformed cells to verify the presence of both “bait” and “prey”. Verified single co-transformed colonies were diluted to optical density OD600 = 1 in distilled water. Equal amounts (5 μl) of each colony were spotted on DDO and quadruple dropout (QDO) medium (without tryptophan, leucine, histidine, and adenine). QDO plates were incubated at 30 °C for 55 h to detect PPIs and photographed. X-alpha-Gal (40 μg/ml) was added in QDO medium for ‘reciprocal’ Y2H (rY2H) (swapping “bait” and “prey” venom proteins between pAD and pBD vectors). At least 3 different colonies from each co-transformation were spotted on both DDO and QDO mediums to validate the reproducibility of results.

Table 1

PCR primer sequences.

Transcript	Primer pair	Primer sequences (5′–3′)
Lys49 PLA2	Lys49-EcoRI FLys49-BamH I R	GGTGAATTCAGCCTGGTCGAATTGGGGAAGGTGGATCCTCATTAGCATGTATCTGGCTTCTT
Asp49 PLA2	Asp49-EcoRI FAsp49-BamH I R	GGTGAATTCAACCTGCTGCAATTCAACAAGGTGGATCCTCATTAGCATTTCTCTGAAGGGT
SVMP I	MPI-EcoRI FMPI-BamHI R	GGTGAATTCGTGAATGATTATGAAGTAGGTGGATCCTCATTAAGCCTCCAAAAGTTCATT
SVMP II	MPII-EcoRI FMPII-BamHI R	GGTGAATTCATAATCCTGGAATCTGGGAGGTGGATCCTCATTAGCCATAGAGGCCATTT
SVMP III	MPIII-EcoRI FMPIII-BamHI R	GGTGAATTCATAATCCTGGAATCTGGGAGGTGGATCCTCACTAGTAGGCTGTAGCCA
C-type lectin	Lect-EcoRI FLect-BamHI R	GGTGAATTCGATTGTCCCTCTGGTTGGTGGATCCTCACTATGCCTTGCAGACGA
3FTx	3FTx-EcoRI F3FTx-BamHI R	GGTGAATTCCTGGAATGTGAAGCATGCAAGGTGGATCCTCATTAAGCGTTGCACAGGTT
CRISP	Cri-EcoRI FCri-BamHI R	GGTGAATTCAGTGTTGATTTTGATTGGTGGATCCTCACTATATTATTTTATTT
Vespryn	Ves-EcoRI FVes-BamHI R	GGTGAATTCGATGTGACGTTTGACTCAAGGTGGATCCTCATTAAAGAGTTGTGAGT
Crotamine	Cro-EcoRI FCro-BamHI R	GGTGAATTCCAATCACAGTGTGAACAGGTGGATCCTCATTATTTTCCAATTTTGCT
LAAO	LAAO-EcoRI FLAAO-BamHI R	GGTGAATTCATGTCTTCTGTGACAGTTGGTGGATCCTCATTAAAATTCATTGTCAT
PLB	PLB-EcoRI FPLB-BamHI R	GGTGAATTCGATATCCACTATGCTAGGTGGATCCTCATCACAGCACTGGTTTCA
EGF*	EGF-pBD FEGF-pBD R	CATGGAGGCCGAATTCCTCTGGGCACCTCCGAGCAGGTCGACGGATCCAAAGTCTTCATGGGTA
EGF*	EGF-pAD FEGF-pAD R	GGAGGCCAGTGAATTCCTCTGGGCACCTCCGACGAGCTCGATGGATCCAAAGTCTTCATGGGTA
SVSP*	Sp-pBD FSP-pBD R	CATGGAGGCCGAATTCGTCGTTGGAGGTGATGAAGCAGGTCGACGGATCCATGGGGGGCAGGTCGCAT
SVSP*	SP-pAD FSP-pAD R	GGAGGCCAGTGAATTCGTCGTTGGAGGTGATGAACGAGCTCGATGGATCCATGGGGGGCAGGTCGCAT

*Homologous recombination-technique was used for both EGF and SVSP due to either EcoRI or BamHI site presenting in the open-reading frames.

Molecular dynamics simulation

Amino acid sequences of Lys49 PLA2 and CRISP were blasted against Protein Data Bank (PDB) (Berman et al., 2007). Crystal structures with PDB codes for 6CE2 and 1RC9 corresponding to Lys49 PLA2 and CRISP were selected as the templates for modeling the 3D structures of Lys49 PLA2 and CRISP by the most widely used template-based protein structure modeling software, MODELLER version 9.24 (Ŝali and Blundell, 1993, Webb and Ŝali, 2014), using the procedures detailed in previous work (Jia et al., 2020). But this time, we simulated 100 candidate models for each venom protein and selected the one with the lowest discrete optimized protein energy (DOPE) score for subsequent protein-protein docking programs. The modeled structures of Lys49PLA2 and CRISP were evaluated through the following criteria: 1) DOPE score (Webb and Ŝali, 2014), the lowest DOPE score model represents the more accurate protein model at its native conformation; 2) overall quality factor (ERRAT) (Colovos and Yeates, 1993), checking if the score is above the expected accuracy of 70% of residues for medium resolution structure; 3) Verify3D (Eisenberg et al., 1997), testing if more than 70% residues having compatibility between the 3D model and the amino acid sequence (1D); 4) Z-score (Wiederstein and Sippl, 2007), testing if the protein model predicted falls within range of high quality experimental structure with a similar size and shape; and 5) Ramachandran (φ/ψ) plot (Ramachandran et al., 1963), checking if interrogated phi and psi dihedral angles of more than 90% of C-alpha residues were within the protein model. All 3D structures of Lys49 PLA2, CRISP and Lys49 PLA2-CRISP were revealed by UCSF Chimera program (Pettersen et al., 2004). Further, three top listed docking programs including ZDOCK (Pierce et al., 2014), HADDOCK (van Zundert et al., 2016) and ClusPro (Kozakov et al., 2017) were used to predict the binding modes of Lys49 PLA2-CRISP interaction using default settings except HADDOCK by activation of random patches. The top docked modes from each docking program were selected for MD simulations by using GROMACS 2020 package (Hess et al., 2008). The CHARMM36 force field (Best et al., 2012) and the standard TIP3P water model (Jorgensen et al., 1983) were chosen for all MD simulations. 1) Each complex was embedded in the center of a dodecahedron water box with a minimum distance from the complex to the box boundary of 10 Å. 2) The complex was solved with water modeled by the TIP3P force field. 3) The system was further neutralized by proper NaCl solution. 4) Energy minimization was carried out for each complex using a steepest-descent integrator to reach negative potential energy and the maximum force for less than 1,000 kJ/(mol.nm) on any atom. 5) To relax the protein complex, the equilibration was performed with positional restraints on all heavy atoms of complex under two ensembles, NVT (constant number of particles, volume and temperature) under constant temperature at 300 K for 100 picoseconds (ps), and NPT (constant number of particles, pressure, temperature) under constant pressure at 1 bar for 100 ps. 6) The production MD simulation was conducted for 1 nanosecond (ns) for initial screening, then 10 ns for further confirmation of structure stability by using the same settings as the previous equilibration but without restraints. The complex structure stability of Lys49 PLA2 interacting with CRISP was assessed by measuring room-mean-square deviation (RMSD) for backbone atoms of each complex. The RMSDs were illustrated with R-program (Team, 2013).

Results

Venom Lys49 PLA2 interacts with CRISP

We screened C. atrox venom proteins encoded by 14 of the most common transcripts for detecting the venom PPIs using the most common venom component, PLA2 including Lys49 PLA2 and Asp49 PLA2 as “bait” by Y2H analysis. The screening results showed that PLA2s interact with themselves, which can be served as perfect internal positive controls in Y2H analysis because native venom PLA2 dimers have been extensively reported, including C. atrox Lys49 PLA2 homodimer (Brunie et al., 1985) and Asp49 PLA2 homodimer (Keith et al., 1981). Moreover, Lys49 PLA2 strongly interacts with venom CRISP, while Asp49 PLA2 weakly interacts with 3FTx (Fig. 1A). To eliminate the false positives, we deployed “rY2H” analysis by swapping “bait” and “prey” proteins between pAD and pBD vectors (i.e., pBD-Lys49 PLA2 × pAD-CRISP, and pBD-CRISP × pAD-Lys49 PLA2). Briefly, 1) we repeated above screening results but individually and with more stringent detection by adding X-alpha-Gal. The results are reproducible in that PLA2S form dimers and Lys49 PLA2 interacts with CRISP (Fig. 1B). 2) To further verify the genuine interaction between Lys49 PLA2s and CRISP, we switched Lys49 PLA2 to pAD and CRISP to pBD vectors (Fig. 1C), yielding the same results as in Fig. 1B except the ‘unilateral’ interaction between Asp49 PLA2 and CRISP (Fig. 1C). Empty vectors (pAD and pBD) served as negative controls showing no interactions with any venom proteins (Fig. 1 A, B and C).

Fig. 1

Yeast two-hybrid analysis. Colony growth on DDO medium (without leucine and tryptophan) indicates the successful binary co-transformations, while colony growth on QDO medium (without histidine, leucine, tryptophan, and adenine) shows the protein-protein interactions (PPIs). A, using pBD-Asp49 PLA2 and pBD-Lys49 PLA2 as “bait” to detect the venom PPIs between PLA2s and 14 venom proteins in pAD vector (empty pBD and pAD served as negative control). B, C, the results of ‘reciprocal’ Y2H, swapping venom proteins from pBD (B) to pAD vector (C), show that PLA2s form dimers and Lys49 PLA2 interacts with CRISP on QDO/X-alpha-Gal medium.

The stable and possible binding mode of Lys49 PLA2-CRISP

Y2H technique can be used to identify PPIs, but it is unable to reveal the complex 3D structure of PPIs. To display the 3D structure of Lys49 PLA2-CRISP interaction, we first built the models for Lys49 PLA2 and CRISP based on the crystallographic structures of MjTX-I (PDB ID 6CE2) from Bothrops moojeni venom (Salvador et al., 2018) and Stecrisp (PDB ID 1RC9) from Trimeresurus Stejnegeri venom (Guo et al., 2004), respectively. The amino acid sequences of the models shared higher sequence identity (more than 70%) and overall sequence coverage of 100% to that of the templates (Table 2). The RMSD values after superimposition and aligning all atoms of models (Lys49 PLA2 and CRISP) to the templates (6CE2 and 1RC9) were 0.154 and 0.200 Å, respectively. This suggests very little deviation between carbon main chain atoms of models and templates, demonstrating that modeled structures are near to the native structures. To further evaluate the accuracy and reliability of modeled structures, we used various stereochemical parameters including ERRAT, vertify3D, Z-Score, and backbone φ/ψ angles. The results shown in Table 3 support that predicted models are of good quality and satisfied the quality checks. After choosing the validated models, we conducted the protein-protein docking, followed by MD simulations. The top complex structures of 11, 9 and 8 were generated by ClusPro, HADDOCK and ZDOCK docking programs, respectively. Subsequently, the great challenge is to identify the correct binding modes of Lys49 PLA2 and CRISP interaction. Fortunately, MD simulation is a well-established technique and can be used for identifying the possible correct binding mode by assessing the time-resolved motions and stability of protein-protein complex structures at atomic resolution (Radom et al., 2018, Sakano et al., 2016, Dror et al., 2012, Perilla et al., 2015). Therefore, we determined which complex structure forms the stable conformation by using MD simulations. All above 28 docked structures were run for 1000 ps (Fig. 2A) by using GROMACS. After MD simulations, the three docking poses, model3, cluster21 and complex6 with average complex backbone RMSDs of 0.229 ± 0.034 nm2, 0.245 ± 0.054 nm2 and 0.363 ± 0.058 nm2, respectively, show the most stable binding modes in each group (Fig. 2D, 2E, 2F). We further conducted MD simulations for these three top complex structures by extending time to 10 ns, resulting that model3 is the most stable complex structure and therefore is the most likely docking mode for Lys49-CRISP interaction, whereas complex6 and particularly cluster21 quickly shifted from initial mode at the time of 2 ns (Fig. 2G). During 10 ns simulations, the average RMSD of backbone atoms of the three systems (Lys49 PLA2, CRISP and model3) indicated that CRISP (0.21 ± 0.02) had very similar RMSD values compared to 0.23 ± 0.02 for model3, whereas Lys49 PLA2 had even lower RMSD (0.14 ± 0.03), suggesting that three systems show similar behavior (Fig. 2G). Additionally, there are 3 hydrogen bonds formatted between Lys49 PLA2 and CRISP within 3 Å in model3, and no inter-residue-hydrogen bonds formatted in cluster21 and complex6 (Fig. 2D, 2E, 2F).

Table 2

BLASTP analysis of Lys49 PLA2 and CRISP and the selected templates.

Models	Templates	Identity (%)	Similarity (%)	Cover (%)	Mas score	E-values
Lys49 PLA2	6CE2	72	80	100	176	4e−58
CRISP	1RC9	81	91	100	391	4e−140

Table 3

Evaluation of Lys49 PLA2 and CRISP models.

Models	DOPE	ERRAT (%)	Verify3D (%)	Z-score	φ/ψ-plot (%)
Lys49 PLA2	−11141.5	82.3	87.6	−5.0	91.4
CRISP	−22896.4	86.6	88.2	−5.7	92.8

DOPE (Discrete Optimized Protein Energy) scores were generated via many iterations by MODELLER (Webb et al., 2014) scripts; more negative DOPE score values tend to correlate with more native-like models. ERRAT scores over 80% indicate that only a few residues have an elevated error function when compared with similar experimental structures. Verify3D scores over 80% signify that the amino acids have compatibility between the 3D model and the amino acid sequence. Z-score is used to assess if the knowledge-based potential could recognize a native fold from the other alternatives (acceptable range from −12 to 12). Ramachandran plot is used to test if protein amino acid residues are in most favored phi (φ) and psi (ψ) dihedral angle region. ERRAT, Verify3D, Z-score and Ramachandran plot were calculated by PROCHECK (Laskowski et al., 1993).

Fig. 2

Prediction of Lys49 PLA2-CRISP complex structure by molecular dynamics (MD) simulations. Room-mean-square deviation (RMSD) of backbone atoms from docked modes of (Lys49-CRISP) as a function of simulation time. A, B and C: During the simulation (1 ns), some complexes in each docking method (ClusPro, ZDOCK and HADDOCK) move away from the initial pose (the trajectory drift away from the initial structure), while the others remain close. The former could be identified as wrong docking models, and the latter as possible correct models (e.g., model3, cluster21, complex6). D, E and F: The 3D complex structures of model3, cluster21 and complex6 (CRISP in magenta, Lys49 PLA2 in blue, and hydrogen bonds denoted by black line in model3), visualized by UCFC Chimera. G: The comparison of structure stability of Lys49 PLA2, CRISP and three complexes in extended time (10 ns). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

Venom protein interaction

Proteins and peptides constitute the complex mixture of snake venom and understanding the venom complexity is critical for developing antivenom for the treatment of snakebites. However, to the best of our knowledge, there is no systematic report of even a single venom complexity, partially due to the lack of high-throughput technique for detecting venom PPIs. Fortunately, since its inception, Y2H analysis has been extensively used to detect any PPIs. Certainly, like all the other approaches, there are some limitations of Y2H analysis. The major drawback of Y2H analysis is that it generates false positives. Many researchers utilized different approaches such as different Y2H systems (Caufield et al., 2012), pull-down analysis (Rimbault et al., 2019), size exclusion (Kirkwood et al., 2013, Busch et al., 2018), etc. as complementary approaches to validate the Y2H results. Among them, a widely used method is to increase the stringency of the interactions. In the present study, apart from using stringent interactions by adding X-alpha-Gal, we deployed “rY2H” method to validate the true interactions. We only consider the PPIs from “rY2H” results as the guanine interactions, e.g., interactions between Lys49 PLA2 and CRISP, as well as Lys49 PLA2 homodimer and PLA2 heterodimer (Fig. 1A, 1B, and 1C). Braun et al., 2009, Caufield et al., 2012 also successfully validated the true PPIs by swapping bait-prey proteins in Y2H systems. PPIs from ‘unilateral’ Y2H need to be further confirmed such as the interactions of Asp49 PLA2 and CRISP in Fig. 1C but not in Fig. 1B, as well as Asp49 PLA2-3FTx interaction in Fig. 1A.

PLA2 enzymes are one of the most common proteins in most snake venoms and are responsible for diverse toxicities including neurotoxicity, myotoxicity, blood coagulation, etc. (Kini, 2003, Lomonte et al., 2003). We predicted that PLA2 interacting with other venom proteins plays a key role in pathophysiological effects. Surprisingly, our findings demonstrated that venom PLA2 only interacts with few venom proteins such as PLA2 itself, CRISP, and possibly with 3FTx. This result implies that venom toxins probably either execute toxic effects individually, or transiently interact with other molecules to exhibit synergistic activity during envenomation. The most common function of snake venom CRISP is to inhibit ion channels, but no snake venom CRISP so far have proved lethal to mammals (reviewed in Tadokoro et al. (2020)), whereas Lys49 PLA2s were extensively reported that they possess myotoxic activity; thus, we speculate that the interaction of Lys49 PLA2 and CRIDP probably produce more potent myotoxicity or strong binding capacity to various target molecules such as ion channels to affect cellular signaling. Our findings are in good agreement with prior conclusion that most snake venom toxins exhibit their pharmacological activities on their own; however, some proteins form covalent/non-covalent complexes with other proteins to exhibit more potent pharmacological activities (reviewed in Doley and Kini (2009)). It appears that snake venom proteins preferentially form homodimers or heterodimers with different members of the same venom protein family, and they scarcely interact with non-homologous venom proteins. Understanding this principle is critical for the development of universal antivenom or complement methods to treat snakebite, the neglected tropical disease that results in high mortality (138,000 deaths per annum) and morbidity (~400, 000 cases per annum) (Gutiérrez et al., 2017).

Complex structure of venom protein interactions

While the successful detection of venom PPIs is important, the study of venom protein interactions at the structural level plays a crucial role in structure-based antivenom development because the pathophysiological effects of venom proteins or protein complexes are known to be closely related to their 3D structures. To predict the complex structure of Lys49 PLA2-CRISP interaction, the accuracy and reliability of modeled individual structures of Lys49 PLA2 and CRISP were critically important for following docking and MD simulations. Both modeled structures of Lys49 PLA2 and CRISP are close to the native structures with less than 0.2 Å RMSD between templates and modeled structures. Furthermore, the various structure parameters indicate that the modeled structures of Lys49 and CRISP are close to the crystal experimental structures (Table 3). Subsequently, based on Critical Assessment of protein Structure Prediction (CASP)-Critical Assessment of PRediction of Interactions (CAPRI) (Lensink et al., 2017, Vangone et al., 2017, Kurkcuoglu and Bonvin, 2019, Agrawal et al., 2019), we employed the three top ranked protein-protein docking methods (ClusPro, ZDOCK and HADDOCK) for predicting the complex structure of Lys49 PLA2-CRISP interaction. To further identify the most likely docking mode of Lys49 PLA2-CRISP interactions, we performed the MD simulations and the results revealed that the model3 (Fig. 2D and G) shows stable complex structure in extended time (10 ns) and exhibits similar stable behavior with Lys49 PLA2 and CRISP (Fig. 2G) when comparing with cluster21 and complex6. MD simulation can evaluate the thermodynamic stability of protein complex structures in time-resolved motions, and the highly stable complex structures should represent the most likely biological form of functional structure (Krissinel and Henrick, 2007). In addition, MD simulation has been tested for several systems to identify the correctly docked modes (Dror et al., 2011, Shan et al., 2011, Buch et al., 2011), and it is commonly accepted that the correctly docked conformations appear to be thermodynamically stable (Quezada et al., 2017, Radom et al., 2018, Bhakat et al., 2018).

The study of venom protein interactions at the structural level plays a crucial role in the investigation of toxic systems and the development of antivenom. However, the current situation in animal venom research is that there are numerous venom protein transcripts and even protein sequences available in various databases. For example, there are 10 transcript sequences of Crotalus CRISP deposited in NCBI GenBank, but only one crystal structure of Crotalus CRISP is available in PDB (accessed in January 2021). Therefore, there is a large gap between the number of available venom protein sequences and their experimentally solved crystal structures. Developing crystal structures of venom proteins is a time consuming and tedious laboratory process. To fill out this gap, computational modeling of the 3D structure of complexes from their protein sequences play a central role in structural venom research. From our current study, we believe that the combination of experimental (Y2H) and computational (modeling) approaches appears promising for not only detections of venom PPIs, but also for dissection of venom complexity at the structural level.

Conclusions

Understanding the venom protein interactions is essential for deciphering venom complexity and designing specific chemical modifications to develop new reagents and therapeutics for the treatment of snakebites. We detected the venom protein-protein interactions of C. atrox and the results showed that venom PLA2s interact with themselves and Lys49 PLA2 interacts with CRISP. The thermodynamically stable conformations of Lys49-CRISP were further predicted by MD simulations. Based on the above and our other preliminary results (data not shown), we propose that, although snake venom is a cocktail of proteins and other molecules, most venom proteins individually exert their pathophysiological effects or are coupled with small molecules such as peptides and chemical compounds. Future studies will entail testing the synergistic activities of Lys49-CRISP complex.

CRediT authorship contribution statement

Paulina Kowalski: Methodology, Formal analysis, Investigation, review & editing. Ivan Lopez: Data curation, software, review & editing. Ying Jia: Supervision, Validation, Writing - original draft, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.