Identification, description and structural analysis of beta phospholipase A2 inhibitors (sbβPLIs) from Latin American pit vipers indicate a binding site region for basic snake venom phospholipases A2.

Several snake species possess, in their circulating blood, endogenous PLA2 inhibitors (sbPLIs) with the primary function of natural protection against toxic enzymes from homologous and heterologous venoms. Among the three structural classes of sbPLIs - named α, β, and γ - the β class (sbβPLIs) is the least known with only four identified sequences, so far. The last class of inhibitors encompass molecules with leucine rich repeats (LRRs) motifs containing repeating amino acid segments. In the present study, we identified and characterized putative sbβPLIs from the liver and venom glands of six Latin American pit vipers belonging to Bothrops and Crotalus genera. The inhibitor from Crotalus durissus terrificus snakes (CdtsbβPLI) was chosen as a reference for the construction of the first in silico structural model for this class of inhibitors, using molecular modeling and molecular dynamics simulations. Detailed analyses of the electrostatic surface of the CdtsbβPLI model and protein-protein docking with crotoxin B from homologous venoms predict the interacting surface between these proteins.

Introduction

Endogenous phospholipase A2 inhibitors are oligomeric glycoproteins secreted by the liver into the circulating blood of several snake species. The primary purpose of these inhibitors, known as sbPLIs (an acronym for snake blood phospholipase A2 inhibitors), is self-protection against the eventual presence of toxic snake venom phospholipases A2 (svPLA2) in the bloodstream.

The identification of sbPLIs from different species led to their classification into three distinct classes − named α, β, or γ − according to the presence of structural domains previously described for mammal proteins (Ohkura et al., 1997). Up to now, the great majority of sbPLIs belongs to α or γ classes (Campos et al., 2016 and references therein). SbαPLIs hold a characteristic C-type lectin-like domain, while sbγPLIs have a particular cysteine pattern that forms three-finger motifs. On the other hand, sbβPLIs display tandem leucine-rich repeats (LRRs) that delineate a horseshoe shaped molecular structure. Such inhibitors have been solely purified from three Asian snake species, being the first of which from the viperid Gloydius brevicaudus (formerly Agkistrodon blomhoffii siniticus). Under its native form, Gloydius brevicaudus sbβPLI (GbsbβPLI) occurs as a trimer composed of heavily glycosylated subunits of 50 kDa each. Each monomer binds to one PLA2 molecule and GbsbβPLI exclusively inhibits basic PLA2 from viperid snake venoms (svPLA2s), showing no inhibitory effects over acidic or neutral svPLA2s (Okumura et al., 1998, Okumura et al., 2002). Two other sbβPLIs have been described in the nonvenomous snakes Elaphe climacophora and Elaphe quadrivirgata: EcsbβPLI and EqsbβPLI, respectively. Although these species lack the secretion of toxic svPLA2s, both inhibitors maintain the narrow specificity of GbsbβPLI against basic svPLA2 (Okumura et al., 2002, Shirai et al., 2009). The sole evidence of sbβPLIs in Latin American snakes so far was the identification of transcripts in the venom glands of Lachesis muta muta (Lima et al., 2011).

Considering the limited number of sbβPLIs described in literature, this study aims at identifying and describing novel sbβPLIs from liver and venom glands of Latin American pit vipers belonging to Bothrops and Crotalus genera. The sbβPLI from C. durissus terrificus (CdtsbβPLI) was used to construct the first in silico structural model of a sbβPLI. Through the molecular docking of CdtsbβPLI and crotoxin B, which is the major basic neurotoxin in the homologous venom, we predicted the amino acids that interact at the interface of the inhibitor-toxin complex.

Materials and methods

Snake tissue collection

Bothrops alternatus, B. jararaca, B. jararacussu, B. moojeni, B. neuwiedi and C. durissus terrificus were obtained from the Serpentarium of Fundação Ezequiel Dias (Belo Horizonte, MG, Brazil). The snakes were anesthetized according to the protocol approved by the Ethics Committee on Animal Use (CEUA/FUNED 022/2012). Liver and venom glands were collected in DEPC-treated tubes, quickly frozen in liquid nitrogen and stored at −80 °C until use.

RNA extraction and cDNA synthesis

Total RNA and cDNA synthesis was performed as previously described for α and γ sbPLIs from Latin American snake species (Estevão-Costa et al., 2008, Estevão-Costa et al., 2016). Briefly, total RNA was extracted from approximately 120 mg of snake tissue (liver, venom gland or both) with Trizol® (Invitrogen, USA). RNA integrity was tested by electrophoresis on 1% agarose gel. For cDNA synthesis we used the SuperScript® III First-Strand Synthesis System (ThermoFisher Scientific). In the first step, cDNA was synthesized with oligo (dT)12-18. In the second step, PCR was performed with sense signal peptide (3′ATGAAGTCTTCGGTTCCATCTC5′) and antisense carboxyl terminus (3′TTAGCAGGGACAAATTTGGGAT5′) oligonucleotides, based on the published nucleotide sequence of EqsbβPLI (Okumura et al., 2002). As an additional control of RNA integrity and PCR performance, the specific sbβPLI primers were replaced by sense and antisense beta-actin primers (Cat. G7540, Promega Co. USA), in the second step of the reaction. All the amplifications were performed in a TC-412 thermocycler (Techne, UK) under the following conditions: 5 min at 94 °C, 35 cycles of 1 min at 94 °C, 30 s at 55 °C and 30 s at 72 °C, followed by an extension period of 7 min at 72 °C. Fresh PCR products were cloned into pGEM-T vector (Promega, USA) according to the manufacturer's instructions. Insert-containing clones were isolated after confirmation by conventional PCR screening of transformed NM522 E. coli in the presence of T7/SP6 promoter primers. Aliquots of the amplification reactions were analyzed by electrophoresis on 1.0% agarose gel in TBE buffer, in the presence of ethidium bromide.

Primary and secondary structure analyses of sbβPLIs

DNA sequences were determined using the Big Dye Terminator Cycle Sequencing Kit on an automated ABI 3130 Genetic Analyzer (Thermo Fisher Scientific) and consensus sequences were obtained from a minimum of three complete reads in both directions. Primary sequence deductions, chemical protein properties calculations, secondary structure prediction and multiple sequence alignments using the ClustalW algorithm were performed with the MacVector 15.1.1 software (Mac Vector Inc., USA). N-glycosylation sites were predicted using NetNGlyc 1.0 software (http://www/cbs.dtu.dk) with default threshold (>0.50).

Molecular modeling and dynamics simulations

The in silico molecular model of C. durissus terrificus sbβPLI (CdtsbβPLI) was generated by threading and molecular dynamics simulation techniques (Saxena et al., 2009). Multiple alignments were performed with salign and align 2 scripts based on 4QXE_A (Score 130–139; E-value: 1e-18; Identity: 24%), 4BV4_R (Score: 128; E-value 1e-16; Identity: 20–23%), 1OZN_A (Score: 123–131; E-value: 1e-17; Identity: 24%) and 1P9A_G (Score: 120; E-value: 1e-16; Identity: 28%) templates. The alignments were manually adjusted, after incorporating the restriction of putative disulfide bonds. Fifty initial models of the selected sbβPLI primary sequence were generated by the Modeller v.9.14 software (Martí-Renom et al., 2000). The best initial model was selected based on QMEAN and Ramachandran plot values at QMEAN and RAMPAGE servers, respectively (Benkert et al., 2009, Lovell et al., 2003). Next, the initial model was submitted to molecular dynamic (MD) simulations using the GROMACS (Groningen Machine for Chemical Simulation) v.4.5.3 software (Pronk et al., 2013) in the presence of explicit water molecules. The protonation states of charged groups were set according to pH 7.0. Counter ions were added to neutralize the system and GROMOS 96 53a6 force field (Oostenbrink et al., 2005) was chosen to perform MD simulations. The minimum distance between any atom of the protein and the box wall was 1.0 nm. Energy minimization (EM) using a steepest descent algorithm was performed to generate the starting configuration of the system. After this step, 200 ps of MD simulations with position restraints applied to the protein (PRMD) were performed in order to gently relax the system. Then, 100 ns of unrestrained MD simulations were calculated to evaluate the stability of the structure. All MD simulations were carried out in a periodic truncated dodecahedron box under constant temperature (298 K) and pressure (1.0 bar), which were maintained by coupling to an isotropic pressure and external heat bath. Search for crystallographic data of sbβPLIs-homologous proteins was performed with the HHPred server (Söding, 2005). The final model was analyzed using Ramachandran plot (Lovell et al., 2003) and its overall quality was evaluated by Z- and QMEAN4 scores (Benkert et al., 2009, Wiederstein and Sippl, 2007). The solvent-accessible surface area (SASA) of the in silico model was calculated by AreaImol software (Saff and Kuijlaars, 1997). Electrostatic potential surfaces were generated by APBS (Adaptative Poisson-Boltzmann Solver) electrostatic calculations (Baker et al., 2001), available at Chimera v.1.10 (Pettersen et al., 2004), after the transformation from PDB to PQR file using the online server PDB2PQR (Dolinsky et al., 2004).

Protein-protein docking predictions

The interface interaction between the CdtsbβPLI in silico model and the crystal structure of crotoxin B (CB) isoform CBc (PDB ID 2QOG) was computationally predicted using docking algorithms available at HADDOCK 2.2 (Van Zundert et al., 2016). The docking protocol consisted in three stages (Dominguez et al., 2003): (i) randomization of orientations around its mass center and rigid body energy minimization (EM), in which each protein is allowed to rotate to minimize the intermolecular energy function. Then, both translation and rotation are allowed, and the proteins are docked by rigid body EM. Typically, 1000 complex conformations are calculated at this stage and the best 200 solutions in terms of intermolecular energies are subsequently refined. (ii) Three simulated annealing refinements: first, both proteins are considered as rigid bodies (1000 steps from 2000 to 50 K with 8 fs time steps); second, the side chains at the interface are allowed to move (4000 steps from 2000 to 50 K with 4 fs time steps); and third, the side chains and the backbone at the interface are allowed to move, allowing conformational rearrangements (1000 steps from 500 to 50 K with 2 fs time steps). (iii) A series of molecular dynamics simulations with explicit solvent, for a final refinement. The final docking solutions are clustered through the pairwise backbone root mean square deviation (I-RMSD) at the interface, and the cluster is defined as an ensemble that displays I-RMSD smaller than 1.0 Å. The resulting clusters ranked on HADDOCK score summarize the average interaction energies (electrostatic interaction energy; van der Waals interactions and restraints violation energy) and their average buried surface area. The interactions of protein-protein predicted surfaces were plotted using the DIMPLOT software (Laskowski and Swindells, 2011).

Results

Detection and characterization of putative sbβPLIs in Latin American pit vipers

The integrity of the starting RNAs was confirmed by the unique presence of two bands that corresponds to 18S and 28S bands of ribosomal RNA, by electrophoresis analysis (data not shown). RT-PCR products with the expected sizes for sbβPLIs (about 1000bp) confirmed the presence of sbβPLIs transcripts in snake tissues (liver, venom glands or both) of B. alternatus (Bal), B. jararaca (Bja), B. jararacussu (Bju), B. moojeni (Bmo), B. neuwiedi (Bne) and C. durissus terrificus (Cdt) (Fig. 1). The species abbreviation is followed by an asterisk when the source tissue was the venom gland.

Fig. 1

Electrophoresis on 1.0% agarose gel of the RT-PCR products after amplification of liver (left side) and venom glands (right side) from Latin American pit vipers, using specific primers for sbβPLIs. The amplification products have about 1000bp. Species abbreviations: Bal - B. alternatus, Cdt - C. durissus terrificus; Bju - B. jararacussu; Bja - B. jararaca; Bmo - B. moojeni; Bne - B. neuwiedi. Snake tissues: liver (left side); venom gland (right side with superscripted *). Act-actin. NC – negative control (no DNA). M - 1 Kb molecular marker (bp sizes are indicated by arrows).

Multiple nucleotide (nt) alignments of the cDNA sequences were performed and the primary structures of putative sbβPLIs were deduced from the consensus nt sequences, according to snake species and tissue source (provided as Supplementary data). Although the RNA and cDNA profiles on gel indicated the presence of transcripts in B. moojeni liver and B. jararacussu venom glands, the nucleotide sequences showed incomplete reads and were not included in our study.

All sbβPLI precursors from Latin American pit vipers displayed a signal peptide comprising 23 amino-acid residues. The mature proteins comprised 309 or 308 amino-acid residues. This difference was due to a proline deletion observed at the 221th position in seven out of the eleven sequences available for comparison (Fig. 2). The average molecular mass was 34667.40 ± 359.90 kDa (mean ± S.D.). The isoelectric points (pIs) varied from 6.4 to 7.6, therefore ensuring slightly acidic to basic character to these molecules. Leucine and proline residues accounted for 18.6 ± 1.0% and 7.3 ± 0.5% of the total amino acid content, respectively. Nine cysteines were present in most sequences, except for Bne and Lmm*, in which Ser210 was replaced by a tenth cysteine.

Fig. 2

Multiple alignment of the primary structures of putative sbβPLIs from Latin American pit vipers. Snake names are abbreviated by the first letter of the genus followed by the first letters of the species or subspecies. Primary sequences were deduced from transcripts of liver and venom glands (marked with * in the last case). The signal sequence is marked by a continuous line. Amino acid identities are represented by dots shaded in medium gray and similarities are shaded in light gray.

Comparison between Latin American and Asian sbβPLIs

The deduced primary structures of sbβPLIs from Latin American pit vipers were compared with known sbβPLIs from Asian species (G. brevicaudus: Gb; E. climacophora: Ec; E. quadrivirgata, subunits A and B: EqA and EqB, respectively) by calculating the identity (IS) and similarity (SS) scores. Both scores were higher for the sbβPLIs from colubrid species when compared to G. brevicaudus viperid. For the whole set of inhibitors, ISs varied from 61.9% (Lmm* vs Cdt) to 98.4% (EqB vs Ec), while SSs ranged from 75.8% (Bmo vs Cdt/Cdt*) to 99.4% (EqB vs Ec). A pair to pair comparison between ISs suggests three clusters: cluster 1 is formed by colubrid sbβPLIs, while cluster 2 encompasses sbβPLIs from viperids, except for the sbβPLIs from Lmm and Bne, which form a third cluster. A bar-coded representation of the primary structures, outlining cysteine, proline and leucine residues, expresses the consensus sequences in each cluster (Fig. 3).

Fig. 3

Top. Graphical representation of the identity scores (ISs) obtained by multiple alignments of the primary sequences of mature sbβPLIs from liver and venom gland (marked with * in the last case) from Asian snake species (white circles) and Latin American snake species (black circles). Suggested IS clusters are circled by dotted lines: c1 (on the left), c2 (centered), and c3 (on the right). Bottom: schematic bar-coded primary structure of representative amino acids of sbβPLIs, where leucine is coded in black, cysteine in red, and proline in blue. Each line represents the consensus sequences of sbβPLIs according to IS grouping.

Among the four Latin American pit viper species that had sbβPLI transcripts isolated from both liver and venom glands (Bal/Bal*, Bja/Bja*, Bne/Bne*, and Cdt/Cdt*), the putative inhibitor from C. d. terrificus (CdtsbβPLI) presented identical consensus sequences for both tissues (IS = 100%). CdtsbβPLI is part of cluster 2 and displayed ISs of 75.4%, 94.5%, and 62.5% with the consensus sequences of clusters 1, 2 and 3, respectively (data not shown). This inhibitor was chosen as a representative of the sbβPLIs from Latin American pit vipers for more detailed molecular studies.

Table 1 summarizes the sbβPLIs identified in species from Asian and Latin America, including the data of the present study.

Table 1

PLA2 inhibitors from snake blood identified in the structural beta class (sbβPLIs) identified, up to now, in species from Asia and Latin America.

Family	Species or subspecies	Source tissue	GenBank	Reference
Colubridae	Elaphe climacophora	Liver	AB462511	Shirai et al. (2009)
	Elaphe quadrivirgata	Liver	AB060637 AB060638	Okumura et al. (2002)
Viperidae	Bothrops alternatus	Liver	MH479016	Present study
		Venom gland	MH479017
	Bothrops jararaca	Liver	MH479019	Present study
		Venom gland	MH479020
	Bothrops jararacussu	Venom gland	MH479018	Present study
	Bothrops moojeni	Venom gland	MH479021	Present study
	Bothrops neuwiedi	Liver	MH479022	Present study
		Venom gland	MH479023
	Crotalus durissus terrificus	Liver	MH479024	Present study
		Venom gland	MH479025
	Gloydius brevicaudus	Liver	AB007198	Ohkura et al., 1997, Okumura et al., 1998
	Lachesis muta muta	Venom gland	–	Lima et al. (2011)

Structural characterization of the sbβPLI from C. durissus terrificus (CdtsbβPLI)

CdtsbβPLI is composed of 308 amino acids, containing 54 leucines (17.53%), 22 prolines (7.14%) and nine cysteines (2.92%). The calculated pI (7.2) is close to neutrality. The in silico prediction of the secondary structure of CdtsbβPLI points 27% of turns (T), 43% of sheets (S) and 17% (H) of helices, besides 13% of H/T, H/S and S/T segments. N102 (N102ASS105) and N209 (N209SLL212) were predicted as N-linked glycosylation sites, with respective scores of 0.55 and 0.59. These segments are highly conserved in clusters 1 and 2 of Latin American sbβPLIs (Fig. 3).

In order to gain insights about the tertiary structure of sbβPLIs and their interaction with basic phospholipases A2, we generated an in silico model for CdtsbβPLI. We searched for crystal structures of CdtsbβPLI-homologous proteins due to the absence of crystallographic data for sbβPLIs. The output structures were the LRR-containing G-protein coupled receptor 4 from Eptatretus burgeri (PDB ID 4QXE), the toll ectodomain from Drosophila melanogaster (PDB ID 4BV4) and the reticulon 4 receptor from Homo sapiens (PDB ID 1OZN). All these crystal structures presented e-values of −20 and IS of about 23% with CdtsbβPLI and were used for the generation of our template by molecular modeling and molecular dynamics simulations (Fig. 4). According to the Ramachandran plot, 93.8% of the residues are distributed in favored and allowed regions. Z- and QMEAN4 scores of −4.10 and −8.28, respectively, indicate an overall good quality for the proposed model.

Fig. 4

Leucine-rich repeat (LRR) motifs array and molecular modeling of the CdtsbβPLI. A. LRRs motifs array in the primary structure. The eight typical LRR motifs and the ninth C-terminal LRR motif associated with a C-terminal cysteine-rich cap are highlighted and the first residue of each LRR motif is indicated; B. Cartoon representation of CdtsbβPLI in silico structural model. Colors: yellow - eight conserved typical LRR motifs; white - ninth LRR motif flanked by C-terminal cysteine-rich cap; purple – cysteines; cyan - predicted N-linked glycosylation sites (N102 and N209). A detailed view of the N-terminal regions shows the enrichment of proline residues (yellow sticks); C and D. Root mean square fluctuation (RMSF) and root mean square deviation (RMSD), respectively, of CdtsbβPLI model during molecular dynamics simulations.

The search for consensus sequences of LRR motifs in the primary structure of CdtsbβPLI revealed eight typical motifs with full conservation of the ionic character of the Highly Conserved Segment (HCS) and high conservation of the Variable Segment (VS) (Table 2; Fig. 4). In five of the eight conserved HCS, the last leucine is replaced by amino acids with the same chemical character.

Table 2

Comparative analysis of the amino acid sequences of eight typical leucine-rich repeat (LRR) motifs and the ninth type 1 C-terminal LRR motif associated to a cysteine-rich cap in the CdtsbβPLI. HCS- highly conserved segment (in bold); VS- variable segment (regular font) Non-conserved residues are in red font. For the ninth LRR motif, the number of conserved residues of CdtsbβPLI is shown in comparison with type 1 C-terminal cysteine-rich region consensus sequence.

LRR motif	Sequence	No. of matching residues of HCS region with consensus sequence	No. of matching residues of VS region with consensus sequence
Consensus sequence	LxxLxLxxNxLxxLPxxoFxzLxx	–	–
1	58LQELHLSNNRL68KSLPSGLFRNLPQ81	5/5	5/5
2	82LHTLDLSRNFL92EDLPPEIFINASS105	5/5	4/5
3	106LTHLSLSENQL116AELRPSWFESLEK129	5/5	4/5
4	130LRILGLDHNQV140KEIPSCFDKLKE153	5/5	4/5
5	154LTSLDLSFNLI164HRLTTGMFSGLDN177	5/5	4/5
6	178LERLVLESNPI188QCIMGRTFHWRPK201	5/5	2/5
7	202LSVLSLKNSSL212THVIMGVFQ-LDQ224	5/5	4/5
8	225LELLDLSDNEF235TTLDPPVHK-PSA247	5/5	2/5
9	248NFSLDLDLSGNWACDCRLEENLLRW(13)FVC(19)CPC308	–	12/15*
Type 1 Cysteine-rich flanking consensus sequence	LxxLxLxxNPxCxCxoWLxxW(9–24)oxC(9–18)CxxP	–	–

The in silico model of CdtsbβPLI is in accordance to a horseshoe-shaped molecule, with the HCS regions located at its concave face, displaying the typical three-residue parallel β-strands of LRR motifs. A ninth LRR motif, flanked by a type I C-terminal cysteine-rich cap, starts at the 248th residue and displays a four-residue cysteine cluster with Cys260/Cys286 and Cys262/Cys306, forming disulfide bridges, besides a free thiol group Cys308 (Fig. 4). Two other cysteines (Cys147 and Cys190) are located at the convex face of the protein and may not be involved in disulfide bridges, since this would provide an unfavorable geometry for the LRR conformation. The LRR motifs and their flanking regions are conserved in every sbβPLIs (Fig. 2, Fig. 5). Remarkably, the N-terminal portion of CdtsbβPLI does not present the main characteristics of the consensus sequence of an N-terminal cysteine-rich flanking domain, despite the existence of a two-residue cysteine cluster with a Cys4/Cys18 disulfide bridge. An enrichment of eight proline residues in the N-terminal region is observed for all sbβPLIs (Fig. 2, Fig. 4, Fig. 5).

Fig. 5

Amino acid sequences alignment of CdtsbβPLI and known sbβPLIs from Asian species. Abbreviations: G. brevicaudus (GbsbβPLI), E. climacophora (EcsbβPLI), and E. quadrivirgata (EqsbβPLIA and B) species. The highly conserved segments (HCS) of the eight leucine-rich repeat (LRR) motifs from the sequences are highlighted in red boxes. The highly conserved core LxxLxLxxNxL of HCS from typical LRR motifs (1–8 in yellow) is indicated below the yellow boxes, where L is Leu, Ile, Val, or Phe; N is Asn, Thr, Ser, or Cys; C is Cys, Ser, or Asn, and x is any residue. The C-terminal LRR motif flanked by a cysteine-rich cap is also highlighted in a yellow box. The cysteines involved in disulfide bridges in the C-terminal (Cys260/Cys286 and Cys262/Cys306) and N-terminal (Cys4/Cys18) are highlighted in green boxes and by C symbols. Proline residues in the N-terminal region are highlighted in purple boxes and by P symbols, evidencing the proline-enrichment of this region. The acidic amino acids (Glu and Asp) and serine residues that form the negatively charged area located at the N-terminal region and LRR motifs 1–6 (see Fig. 6) in the concave face of the protein are highlighted in red boxes and by E, D and S symbols.

Fig. 6

Electrostatic potential surface analysis of the in silico structural model of CdtsbβPLI (A) and crotoxin B (isoform CBc; PDB ID 2QOG) (B). Red and blue-colored regions denote negative and positive charges, respectively. A. Detailed view of the negatively charged area located at the N-terminal region and LRR motifs 1–5 in the concave face of CdtsbβPLI in silico model, pointed as the binding region to basic PLA2s.The acidic amino acid (Glu and Asp) and serine residues located in this region of the protein are highlighted in sticks. B. Detailed view of the positively charged areas located in the β-wing (residues Arg74, Lys78 and Lys87), N (residues Lys7, Lys10, Arg14 and Lys16) and C-terminal regions (Lys114, Lys125 and Lys127), besides the Arg36, Arg38 and Arg43 residues (highlighted in sticks).

A cluster of acidic amino acids can be noted in the concave surface of the inhibitor, therefore ascribing a negative net charge to it. This charged surface is compatible with the previously reported specificity of sbβPLIs from Asian species to basic svPLA2 (Fig. 4).

Mapping the interface region between CdtβPLIs and basic svPLA2s

A large negatively-charged area mostly constituted by aspartic and glutamic acid residues was observed at the concave surface of the CdtsbβPLI in silico model in the N-terminal region, as well as in the six initial LRRs motifs (Fig. 6). This area also contains several serine residues, which are usually involved in hydrogen bonds as proton donors or acceptors. The residues that form this negatively-charged area in the in silico model display a large solvent-accessible surface area (SASA), exposing an acidic area of ∼730 Å2. This region is in the neighborhood of the above-mentioned proline-rich region at the N-terminal portion of CdtsbβPLI, probably offering structural stability for the binding of basic svPLA2. Remarkably, 10-50 N-terminal region presented the lowest root mean square fluctuation (RMSF) values of the overall structure during the molecular dynamic simulations (Fig. 4).

Comparative amino acid sequence analyses between CdtsbβPLI and other putative sbβPLIs from Latin American pit vipers showed that these acidic (aspartic and glutamic acids) and serine residues have great level of conservation in all sequences, except for a few residues (Fig. 2). Ser38, Glu40 and Ser64 are substituted by Tyr38, Gln40 and Thr64 residues in B. neuwiedi and L. muta muta sequences, respectively. However, these substitutions maintain the acidic character of the residue position. Other substitutions are Glu60/Gly60 and Glu184/Lys184 in B. neuwiedi and L. muta muta sequences, Ser112/Ile112 in B. alternatus venom gland sequence, and D136/Asn136 in B. moojeni venom gland sequence. Amino acid sequence comparisons between CdtsbβPLI and sbβPLIs from Asian snakes also revealed a great level of conservation of the residues present in this negatively-charged area, with only two substitutions: Ser110/Pro110 in GbsbβPLI and Ser154/Phe154 in non-venomous snake sequences (Fig. 5). Therefore, this negatively-charged area is present in all sbβPLIs that have been identified so far and possibly comprises the binding site region for positively-charged basic svPLA2s.

Unfortunately, no native CdtsbβPLI was available for experimental tests. Nevertheless, based on the expected selectivity of sbβPLIs for basic svPLA2, it is possible to suggest that CdtsbβPLI may target crotoxin B (CB), the major basic PLA2 in C. d. terrificus venom. Analyses of the electrostatic surface in the crystal structure of the monomeric CB (PDB ID 2QOG) showed that positively charged areas, composed of lysine and arginine residues, mainly concentrate in the β-wing (residues 75–84), N- and C-terminal regions (residues 1–20 and 115–133, respectively), besides Arg36, Arg38 and Arg43 residues (Fig. 6).

In order to gain insights on the interface region between sbβPLIs and basic svPLA2s, we performed docking predictions between CdtsbβPLI in silico model and the crystal structure of crotoxin B (CB) CBc isoform (PDB ID 2QOG) to obtain an in silico structural model for the CdtsbβPLI/CB complex (Fig. 7). For the in silico model of CdtsbβPLI, the previously identified negatively charged area formed by Glu15, Ser38, Glu40, Glu60, Ser64, Asp86, Thr84, Ser88, Ser110, Ser112, Glu113, Asp136, Ser156, Asp158, Ser160 and Glu184 residues was chosen as an active region to drive the docking (directly involved in the interaction) (Fig. 6). Regarding CB crystal structure, two different positively charged areas were chosen as active regions for dock driving comparison. Region 1 is formed by Lys7, Lys10, Arg14, Lys16, Arg74, Arg78 and Lys87 residues, located at the N-terminal and β-wing regions from CB; while region 2 is formed by Arg36, Lys38, Arg43 and Lys114, Lys125 and Arg127, located at the C-terminal portion of CB (Fig. 6). The passive residues (surrounding surface residues) are automatically defined around the active residues. Comparative analyses between the statistical values for the best cluster among the solutions found for docking predictions of regions 1 and 2 of CB with CdtsbβPLI in silico model showed that region 1 has the best HADDOCK score and the best free energy values for binding (Table 3). These data suggest that the N-terminal and β-wing regions from CB may constitute the region for interaction with CdtsbβPLI. A protein/protein interface analysis of the structural model of CdtsbβPLI/CB complex (Fig. 7) shows that this interface contains fourteen hydrogen bonds. Five residues from the β-wing, a neighboring residue (Lys87) and four residues from the N-terminal portion of CB establish nine, one and four hydrogen bonds, respectively, with serine and aspartic/glutamic acid residues located at the negatively charged area from CdtsbβPLI in LRR motif 1 (one residue/one hydrogen bond), LLR3 (one residue/one hydrogen bond), LRR4 (one residue/one hydrogen bond), LRR5 (two residues/three hydrogen bonds), LRR6 (two of these residues and the exception Arg180/six hydrogen bonds) and even in LRR8 (one residue/one hydrogen bond), which is located outside the negatively charged area. Besides hydrogen bonding, hydrophobic contacts are also established in the CdtsbβPLI/CB interface by three residues from the β-wing, two from its vicinities, six from the N-terminal region and three (Phe24, Trp31 and Trp70) from other regions of CB crystal structure; and two residues from LRR1, four from LRR2, two from LRR3, two from LRR4, three from LRR5, one from LRR6 and two from LRR7 of CdtsbβPLI in silico model (Fig. 7).

Fig. 7

Protein-protein docking predictions between the in silico structural model of CdtsbβPLI and crotoxin B (isoform CBc; PDB ID 2QOG). A. Cartoon representation of the best cluster of docking solutions for CdtsbβPLI/CB structural model. A negatively charged area from the structural model of CdtsbβPLI and two positively charged areas from the crystal structure of CB were chosen as the active residues (directly involved in binding) to drive the docking. These residues in CdtsbβPLI are in green; the N-terminal and β-wing regions from CB are in brown and cyan, respectively (see text for details). The eight conserved typical LRR motifs from CdtsbβPLI structure are highlighted in yellow. B. Hydrogen bonds (black dashes) established between serine and aspartic/glutamic acid residues (highlighted in yellow sticks; exception for Arg180) from CdtsbβPLI (in green) with residues located at N-terminal and β-wing regions (highlighted in brown sticks) from CB crystal structure (in cyan). The hydrogen bond distances (Å) and the eight typical conserved LRR motifs from CdtsbβPLI are indicated. C. Diagram of interactions established at CdtsbβPLI/CB interface. The residues involved in hydrogen bonds (green dashes with distances indicated in Å) from CdtsbβPLI and CB structures are shown in green and purple, respectively. The hydrophobic contacts are shown as arcs with radiation spokes, where the red and purple arcs indicate, respectively, the residues from CdtsbβPLI and CB structures.

Table 3

Statistics of the best cluster of solutions found by macromolecular docking predictions between the in silico model CdtsbβPLI and crotoxin B (CB). The crystal structure of CBc isoform (PDB ID 2QOG) was used in the prediction. RMSD- Root mean square deviation.

Proteins and parameters		CdtsbβPLI in silico model
		HADDOCK score (a.u.)	Cluster size	RMSD from the overall lowest-energy structure (Å)	Van der Waals energy (kcal/mol)	Electrostatic energy (kcal/mol)	Desolvation energy (kcal/mol)	Restraints violation energy (kcal/mol)	Buried Surface Area (Å2)	Z-score (a.u.)
CB crystal structure	R1	−215.7 ± 7.7	23	1.1 ± 0.4	−79.2 ± 4.1	−521.1 ± 52.8	−38.0 ± 6.2	57.7 ± 38.21	2637.3 ± 65.9	−1.6
	R2	−188.4 ± 4.0	18	0.9 ± 0.6	−49.1 ± 7.7	−687.3 ± 60.1	−10.9 ± 7.9	90.7 ± 33.15	2387.5 ± 107.7	−1.8

Discussion

GbsbβPLI was the first prototype for the structural model of blood plasma β class inhibitors from snakes, in spite of the absence of any structural model of its tertiary or quaternary structures. The native trimer (160 kDa) is composed of homogenous subunits of 50 kDa, as estimated by SDS-PAGE (Ohkura et al., 1997). It was reported to be heavily glycosylated with N-linked, since the estimated molecular mass of the monomer dropped from 50 to 39 kDa after enzymatic deglycosylation (Okumura et al., 1998). Two (N102 and N209) of the four potential N-linked glycosylations proposed for GbsbβPLI are present in CdtsbβPLI, therefore confirming the glycosylated nature of sbβPLIs. The average molecular mass of 34.7 kDa calculated for the sbβPLIs monomers described herein is compatible with that of deglycosylated GbsbβPLI monomer (34.6 kDa). Analysis of GbsbβPLI primary structure showed that this protein had LRR motifs homologous to leucine-rich α2-glycoprotein (Okumura et al., 1998).

The in silico model of CdtsbβPLI reveals eight typical LRR motifs in a horseshoe-shaped structure, and a ninth LRR motif flanked by a type I C-terminal cysteine-rich cap (Fig. 4, Fig. 5). Two cysteines (Cys147 and Cys190) located at the convex face of the protein probably occur as free thiols, in agreement with the GbsbβPLI description, since an intra-chain disulfide bond would cause an unfavorable geometry for the LRR conformation. Basic PLA2 molecules would supposedly bind to the concave face of the GbsbβPLI due to negatively charged residues at LLR1 (Okumura et al., 1998). Moreover, the LRR1 region might be responsible for the specific binding to G. brevicaudus basic PLA2, since it is fully conserved among three sbβPLIs from Asian species (Shirai et al., 2009). Curiously, GbsbβPLI gene expression in the liver was shown to be upregulated by acidic svPLA2 and not by basic svPLA2, as could be expected based on the selectivity of the inhibitor to the ionic character of svPLA2 (Kinkawa et al., 2010).

The in silico model of CdtsbβPLI highlights the relevance of a negatively charged area in its concave surface to the specificity of sbβPLIs to basic svPLA2s. Furthermore, the structural model shows that this negatively charged area is not restricted to LRR1 as previously suggested, but it spreads to the N-terminal region and LRR motifs 1 to 6 (Fig. 5). The clustering of negative amino acid residues results in the exposition to solvent of an acidic-charged area of ∼730 Å2. The availability of the concave surface of the inhibitor to svPLA2 binding is in accordance with the most typical profile of LRR-containing proteins. This negatively charged area in the CdtsbβPLI is in the neighborhood of a proline-rich region at the N-terminal portion of the molecule, which can confer structural stability to the binding of basic svPLA2 (Fig. 3). Remarkably, this area in our model is larger than the positively charged area pointed as the acidic svPLA2 binding region at the monomer of αBaltMIP, a sbαPLI from Bothrops alternatus blood (∼540 Å2) (Estevão-Costa et al., 2016).

CdtsbβPLI may target crotoxin B (CB), the major basic PLA2 in C. d. terrificus venom. CB has an isoelectric point of 8.4 and displays widely spread positively charged areas in its structure (Fernandes et al., 2017, Marchi-Salvador et al., 2008) mainly concentrated in the β-wing, N- and C-terminal regions, besides Arg36, Arg38 and Arg43 residues (Fig. 6). Both N- and C-terminal regions have been pointed out to be involved in the neurotoxic effect induced by CB (Čurin-Šerbec et al., 1994, Fernandes et al., 2017, Fortes-Dias et al., 2009). Docking predictions between CdtsbβPLI in silico model and CB crystal structure showed that, differently from the previous studies with GbsbβPLI that suggested LRR1 as the central point for the binding of basic svPLA2, the negatively charged area located at the N-terminal region and LRR motifs 1 to 6 of CdtsbβPLI establishes hydrogen bonds and hydrophobic contacts with CB, especially LRR motifs 2 (4 hydrogen bonds), 5 (3 hydrogen bonds and 3 residues involved in hydrophobic interactions) and 6 (1 hydrogen bond and 6 residues in hydrophobic interactions) (Fig. 7). It is important to point out that this binding region is similar to the binding site of MD-2 at the concave face of the ectodomain of Toll-like Receptor 4 (TLR4), at the junction of the N-terminal and central LRR motifs (Kim et al., 2007). Regarding CB, N-terminal and β-wing regions may constitute the region of interaction with CdtsbβPLI, displaying 14 hydrogen bonds with serine and aspartic/glutamic acid residues located at the negatively charged area of CdtsbβPLI, and 14 residues involved in hydrophobic interactions (Fig. 7).

The identification of sbPLIs from three different structural classes in non-venomous snakes raises questions about other physiological roles of these proteins beyond a natural resistance against the venom. In the particular case of sbβPLIs, Okumura et al. (1998) raised the hypothesis that GsbβPLI could correspond to snake LRG and that mammalian LRG might function as an svPLA2 inhibitor. Although the last assumption has not been proved, GbsbβPLI was shown to bind cytochrome c (Cyt c), an endogenous ligand of LRG, with higher affinity (Kd = 2.37 × 10−12 M) than to basic PLA2 from G. brevicaudus venom (Kd = 1.21 × 10−9 M) (Shirai et al., 2010). In fact, the tight binding of GbsbβPLI to Cyt c led the authors to suggest that the original function of sbβPLIs might be the clearance of autologous Cyt c from the circulatory system of snakes. They also proposed that the sbβPLIs may have acquired a new function as PLA2 inhibitory proteins in addition to the ancestral Cyt c binding, during snake evolution. This possibility indicates a great biotechnological potential of sbβPLIs, beyond the action of protecting against toxic svPLA2, since Cyt c is involved in apoptosis induction and proinflammatory mediation (Hiraoka et al., 2004, Pullerits et al., 2005). Cyt c has a positively-charged lysine-rich region on its surface involved in electrostatic interactions with subunit II of cytochrome c oxidase (Döpner et al., 1999). Thus, Cyt c may bind to sbβPLIs the same way as it binds to basic svPLA2s, by electrostatic interactions with the negatively-charged area. The structure identified in the CdtsbβPLI in silico model has high level of conservation among known sbβPLIs (Fig. 2, Fig. 5, Fig. 6).

CdtsbβPLI, as well as other previously described sbβPLIs, belong to the ‘typical’ subfamily of the structural class III of repetitive proteins, which contain an invariant eleven-residue core LxxLxLxxCxxL plus a variable segment (Matsushima et al., 2007). The whole repetitive segment is formed by the 24-25-residue sequence LxxLxLxxNxLxxLPxxOFxZxLxx, where, in the case of type C-terminal cysteine region identified in CdtsbβPLI, L is Leu, Ile, Val, or Phe; N is Asn, Thr, Ser; C is Cys, P is Pro but can be substituted by Trp or Phe residues in some sequences (Kobe and Kajava, 2001, Ng et al., 2011). SbβPLIs belong to type 2 α/β solenoid structures, which have been associated to the PFAM clan Cl022 of LRRs (Paladin and Tosatto, 2015). Since LRRs occur in a large number of proteins in plants, invertebrates and vertebrates, and participate in several biologically important processes (Ng et al., 2011), we believe that the present study contributes, not only to the study of sbPLIs, but to the general knowledge on the widely distributed LRR-containing proteins.

Conclusion remarks

In the last years, efforts have been made to determine the structure of endogenous phospholipase inhibitors and their binding mode with phospholipases A2 due to the potentiality of these endogenous molecules to be used as complements of conventional serum therapy and/or as inhibitors of secretory PLA2 in pathological processes in humans. Regarding sbβPLIs, few of these molecules have been identified so far and the reasons for their specific inhibition of basic svPLA2 are poorly understood. In the present study, we identified precursors of novel sbβPLIs in six species of Latin American pit vipers, and an extensive structural analysis of the transcribed proteins led to the construction of the first structural model of a sbβPLI, using the inhibitor from C. d. terrificus (CdtsbβPLI) as a prototype. We investigated the electrostatic surface of CdtsbβPLI in silico model and identified a conserved negatively charged area located at the N-terminal region and LRR motifs 1–6 (in the neighborhood of a proline-rich region). Docking predictions between CdtsbβPLI in silico model and CB crystal structure highlighted the role of this area in the binding interface between these two molecules, where LRR motifs 2, 5 and 6 from CdtsbβPLI establish the major number of contacts with CB. We also emphasized the positively charged areas at the surface of the basic svPLA2 from C d. terrificus venom crotoxin B (CB), and pointed that the N-terminal and β-wing portions of this toxin may constitute the interacting interface with CdtsbβPLI. Our data contribute to a better understanding of sbβPLIs as well as of LRR-containing proteins in general.

Conflict of interest

On behalf of the authors of the manuscript entitled “Identification, description and structural analysis of beta phospholipase A2 inhibitors (sbβPLIs) from Latin American pit vipers indicate a binding site region for basic phospholipases A2 from snake venoms”, I declare no conflicts of interest regarding the manuscript.

Financial support

We thank the Brazilian funding agencies: Fundação de Amparo à Pesquisa de Minas Gerais (), Fundação de Amparo à Pesquisa do Estado de São Paulo (), Conselho Nacional de Desenvolvimento Científico e Tecnológico (), Coordenação de Aperfeiçoamento de Pessoal de Ensino Superior ( Toxinology 1810/11; Programa Nacional de Pós-Doutorado) and Instituto Nacional de Ciência e Tecnologia em Toxinas (INCTTox/CNPq) for financial support.