Individual Variability in Bothrops atrox Snakes Collected from Different Habitats in the Brazilian Amazon: New Findings on Venom Composition and Functionality

Differences in snake venom composition occur across all taxonomic levels and it has been argued that this variation represents an adaptation that has evolved to facilitate the capture and digestion of prey and evasion of predators. Bothrops atrox is a terrestrial pitviper that is distributed across the Amazon region, where it occupies different habitats. Using statistical analyses and functional assays that incorporate individual variation, we analyzed the individual venom variability in B. atrox snakes from four different habitats (forest, pasture, degraded area, and floodplain) in and around the Amazon River in Brazil. We observed venom differentiation between spatially distinct B. atrox individuals from the different habitats, with venom variation due to both common (high abundance) and rare (low abundance) proteins. Moreover, differences in the composition of the venoms resulted in individual variability in functionality and heterogeneity in the lethality to mammals and birds, particularly among the floodplain snakes. Taken together, the data obtained from individual venoms of B. atrox snakes, captured in different habitats from the Brazilian Amazon, support the hypothesis that the differential distribution of protein isoforms results in functional distinctiveness and the ability of snakes with different venoms to have variable toxic effects on different prey.


Introduction
Snake venoms are complex mixtures of toxic components belonging to multipleprotein families [1], each of which expresses several isoforms that are present in the venoms in different proportions [2]. The concentration of each isoform present in the venoms is highly variable and modulates venom function [3] and, as a consequence, the complexity groups and suggest that the functional diversity of the venoms appears to be relevant to the ability of these snakes to persist in highly variable unsettled environments, such as the floodplain habitat.

Individual Variation of Venom Composition Is Associated with Habitats
Chromatographic analyses by RP-HPLC were our first approach to investigate the individual variability of B. atrox venoms among the specimens collected in each habitat: forest, pasture, degraded area, and floodplain. Overall, similar profiles can be seen in all chromatograms, although the relative heights or areas of peaks varied greatly between venoms even among snakes collected at the same habitat ( Figure 1 and Supplementary Figure S1, for individual chromatograms). In all venom chromatographies, the highest peaks were eluted after 80 min, while peaks eluted between 50 and 70 min showed the most variable area percentages, being higher in floodplain venoms and lower in pasture venoms. According to our previous study [15], Snake Venom Metalloproteinases (SVMPs) are the main components eluted after 80 min, while Phospholipases A 2 (PLA 2 s), C-type Lectin-like proteins (CTLs), and Snake Venom Serine proteinases (SVSPs) are the prevalent protein families in the fractions eluted in the intermediate steps (50-70 min). supporting our initial assumption. Moreover, our findings provide a deeper view of the main toxins and biological activities related to the individual venom variability within these B. atrox groups and suggest that the functional diversity of the venoms appears to be relevant to the ability of these snakes to persist in highly variable unsettled environments, such as the floodplain habitat.

Individual Variation of Venom Composition Is Associated with Habitats
Chromatographic analyses by RP-HPLC were our first approach to investigate the individual variability of B. atrox venoms among the specimens collected in each habitat: forest, pasture, degraded area, and floodplain. Overall, similar profiles can be seen in all chromatograms, although the relative heights or areas of peaks varied greatly between venoms even among snakes collected at the same habitat ( Figure 1 and Supplementary Figure S1, for individual chromatograms). In all venom chromatographies, the highest peaks were eluted after 80 min, while peaks eluted between 50 and 70 min showed the most variable area percentages, being higher in floodplain venoms and lower in pasture venoms. According to our previous study [15], Snake Venom Metalloproteinases (SVMPs) are the main components eluted after 80 min, while Phospholipases A2 (PLA2s), C-type Lectin-like proteins (CTLs), and Snake Venom Serine proteinases (SVSPs) are the prevalent protein families in the fractions eluted in the intermediate steps (50-70 min).  atrox snakes collected at the forest, pasture, degraded area, or floodplain were applied to a Vydac C-18 column. Mobile phases used were 0.1% TFA in water (solution A) or 0.1% TFA in acetonitrile (solution B). Proteins were gradient-eluted at 2 mL/min (5% B for 5 min, 5-15% B over 10 min, 15-45% B over 60 min, 45-70% B over 10 min, 70-100% over 5 min, and 100% B over 10 min). Separation was monitored at 214 nm.

Population-Level Differentiation between Spatially Disparate B. atrox Groups
First, we identified the population-level differentiation among venoms of snakes collected at the different habitats, based on the individual venom composition of each group of snakes. The PERMANOVA analysis (Table 1) showed that the environment explained 22.5% of the total variation in venom peak abundances in the sampled snakes (df = 3.33, F = 3.2, p < 0.0001). Further exploration of these differences with compositionally robust PCA (Figure 2, score plots) showed that PC1 largely differentiates floodplain from pasture venoms. Forest and degraded venoms occupying an intermediate space on the PC1 axis. Components with the most prominent loadings are fraction 10 and the fractions eluted after 85 min. The highest negative values of PC1, responsible for clustering the pasture venoms, were fraction 23 (PC1 = −0.7308, PC2 = −0.6035) and fraction 21 (PC1 = −0.5028, PC2 = −0.6410). These fractions were characterized in our previous study [15] as PIII-class and PI-class SVMP isoforms, respectively. The highest positive values of PC1, responsible for clustering the floodplain venoms, were fraction 20 (PC1 = 0.2958, PC2 = −0.3732), fraction 24 (PC1 = 0.1241, PC2 = 0.2006), and fraction 10 (PC1 = 0.2284, PC2 = −0.0353), which in our previous study have been characterized as different isoforms of the PIII-class and PI-class SVMPs, with minor proportions of CTL isoforms, and an SVSP, respectively ( Figure 2, loading plots). pasture, degraded area, or floodplain were applied to a Vydac C-18 column. Mobile phases used were 0.1% TFA in water (solution A) or 0.1% TFA in acetonitrile (solution B). Proteins were gradienteluted at 2 mL/min (5% B for 5 min, 5-15% B over 10 min, 15-45% B over 60 min, 45-70% B over 10 min, 70-100% over 5 min, and 100% B over 10 min). Separation was monitored at 214 nm.

Population-Level Differentiation between Spatially Disparate B. atrox Groups
First, we identified the population-level differentiation among venoms of snakes collected at the different habitats, based on the individual venom composition of each group of snakes. The PERMANOVA analysis (Table 1) showed that the environment explained 22.5% of the total variation in venom peak abundances in the sampled snakes (df = 3.33, F = 3.2, p < 0.0001). Further exploration of these differences with compositionally robust PCA (Figure 2, score plots) showed that PC1 largely differentiates floodplain from pasture venoms. Forest and degraded venoms occupying an intermediate space on the PC1 axis. Components with the most prominent loadings are fraction 10 and the fractions eluted after 85 min. The highest negative values of PC1, responsible for clustering the pasture venoms, were fraction 23 (PC1 = −0.7308, PC2 = −0.6035) and fraction 21 (PC1 = −0.5028, PC2 = −0.6410). These fractions were characterized in our previous study [15] as PIII-class and PI-class SVMP isoforms, respectively. The highest positive values of PC1, responsible for clustering the floodplain venoms, were fraction 20 (PC1 = 0.2958, PC2 = −0.3732), fraction 24 (PC1 = 0.1241, PC2 = 0.2006), and fraction 10 (PC1 = 0.2284, PC2 = −0.0353), which in our previous study have been characterized as different isoforms of the PIII-class and PI-class SVMPs, with minor proportions of CTL isoforms, and an SVSP, respectively ( Figure 2, loading plots).  In addition, posthoc comparisons of individual population pairs using PERMANOVA showed that the spatially proximate degraded and forest habitats were the only locations between which venom composition did not significantly differ. All other population pairs differed significantly, with the greatest degree of differentiation between the spatially disparate floodplain and pasture populations.
Next, peak-by-peak comparisons showed that several chromatographic peaks contributed to population-level differentiation among venoms. After FDR-correction, nine ), pasture ( In addition, posthoc comparisons of individual population pairs using PERMANOVA showed that the spatially proximate degraded and forest habitats were the only locations between which venom composition did not significantly differ. All other population pairs differed significantly, with the greatest degree of differentiation between the spatially disparate floodplain and pasture populations.
Next, peak-by-peak comparisons showed that several chromatographic peaks contributed to population-level differentiation among venoms. After FDR-correction, nine ), floodplain ( In addition, posthoc comparisons of individual population pairs using PERMANOVA showed that the spatially proximate degraded and forest habitats were the only locations between which venom composition did not significantly differ. All other population pairs differed significantly, with the greatest degree of differentiation between the spatially disparate floodplain and pasture populations. Next, peak-by-peak comparisons showed that several chromatographic peaks contributed to population-level differentiation among venoms. After FDR-correction, nine ), and degraded ( In addition, posthoc PERMANOVA showed that the only locations between w other population pairs differ between the spatially dispara Next, peak-by-peak co contributed to population-lev peaks showed significant var ) habitats.
In addition, posthoc comparisons of individual population pairs using PERMANOVA showed that the spatially proximate degraded and forest habitats were the only locations between which venom composition did not significantly differ. All other population pairs Toxins 2021, 13, 814 5 of 18 differed significantly, with the greatest degree of differentiation between the spatially disparate floodplain and pasture populations.
In addition, posthoc comparisons of individual population pairs using PERMANOVA showed that the spatially proximate degraded and forest habitats were the only locations between which venom composition did not significantly differ. All other population pairs differed significantly, with the greatest degree of differentiation between the spatially disparate floodplain and pasture populations.
Next, peak-by-peak comparisons showed that several chromatographic peaks contributed to population-level differentiation among venoms. After FDR-correction, nine peaks showed significant variation: Peaks 1, 3,5,7,8,12,20,21, and 23 ( Figure 3). The centered log-ratio mean abundance for each reversed-phase high-performance liquid chromatography peak was plotted. Negative numbers correspond to low-abundance peaks, whereas positive numbers correspond to high abundance peaks. (*) Asterisks indicate peaks that show significant population-level variation. X-axis labels correspond to RP-HPLC peak numbers (1 to 26).
In terms of comparisons between habitats, venoms from the forest or recently degraded area were similar in relative peak abundances. Pasture venom was distinct with a higher abundance of Peaks 5 and 7, which correspond to the acidic PLA2s, Peak 12, with CTL as the major toxin, and Peaks 21 and 23, which include the SVMPs. Significantly lower abundances were observed in pasture venoms for Peaks 8 and 20, which have as major toxins PLA2s and PI-SVMPs, respectively. Peak 20 was practically absent in pasture venoms ( Figure 1). Floodplain venoms were distinct in peaks related to SVMPs and PLA2s: Peaks 1, 21, and 23, which contain SVMPs and disintegrins, were present at a lower abundance, while Peak 20, which contains mostly PI-SVMPs, is at a comparatively higher abundance. In peaks containing PLA2s, Peak 3, which contains K-49 basic PLA2s, is less abundant, while Peak 7, which contains acidic PLA2s, and Peak 8, rich in D-49 basic PLA2s, were proportionally higher in the venoms from the floodplain.

Venom Differentiation Is Not Limited to Rare (Low Abundance) Proteins
We also investigated whether the venom variability in B. atrox groups would be associated with the degree to which the venom proteins were expressed in the different habitats. RP-HPLC peaks for each habitat were classified into two sets: low or high abundance, based on the clr-transformed mean of each peak and a nonparametric MANOVA analysis, as previously described [23], where rare proteins = the mean of an  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 clr- The centered log-ratio mean abundance for each reversed-phase high-performance liquid chromatography peak was plotted. Negative numbers correspond to low-abundance peaks, whereas positive numbers correspond to high abundance peaks. (*) Asterisks indicate peaks that show significant populationlevel variation. X-axis labels correspond to RP-HPLC peak numbers (1 to 26).
In terms of comparisons between habitats, venoms from the forest or recently degraded area were similar in relative peak abundances. Pasture venom was distinct with a higher abundance of Peaks 5 and 7, which correspond to the acidic PLA 2 s, Peak 12, with CTL as the major toxin, and Peaks 21 and 23, which include the SVMPs. Significantly lower abundances were observed in pasture venoms for Peaks 8 and 20, which have as major toxins PLA 2 s and PI-SVMPs, respectively. Peak 20 was practically absent in pasture venoms ( Figure 1). Floodplain venoms were distinct in peaks related to SVMPs and PLA 2 s: Peaks 1, 21, and 23, which contain SVMPs and disintegrins, were present at a lower abundance, while Peak 20, which contains mostly PI-SVMPs, is at a comparatively higher abundance. In peaks containing PLA 2 s, Peak 3, which contains K-49 basic PLA 2 s, is less abundant, while Peak 7, which contains acidic PLA 2 s, and Peak 8, rich in D-49 basic PLA 2 s, were proportionally higher in the venoms from the floodplain.

Venom Differentiation Is Not Limited to Rare (Low Abundance) Proteins
We also investigated whether the venom variability in B. atrox groups would be associated with the degree to which the venom proteins were expressed in the different habitats. RP-HPLC peaks for each habitat were classified into two sets: low or high abundance, based on the clr-transformed mean of each peak and a nonparametric MANOVA analysis, as previously described [23], where rare proteins = the mean of an individual peak < the geometric mean, and abundant proteins = mean of an individual peak > the geometric mean ( Figure 4). Proteins showing deviations off the middle line are those in which the expression levels are different between the compared populations whereas those that are on the line are similar.
individual peak < the geometric mean, and abundant proteins = mean of an individual peak > the geometric mean ( Figure 4). Proteins showing deviations off the middle line are those in which the expression levels are different between the compared populations whereas those that are on the line are similar. Figure 4. Comparisons of the degree of differentiation in low vs. high abundance proteins based on chromatographic peak areas from B. atrox individual venoms. The clr mean was plotted for each RP-HPLC peak across the different habitat (x-axis) and (y-axis) populations for B. atrox snakes. Bars indicate the SE, a solid line indicates a perfect agreement, dashed lines indicate the origin (i.e., the geometric mean), and proteins less than these values were considered low-expression proteins.
In our PERMANOVA analyses, only in the spatially proximate degraded and forest habitats the venom composition did not significantly differ (Table 1). Accordingly, minor Comparisons of the degree of differentiation in low vs. high abundance proteins based on chromatographic peak areas from B. atrox individual venoms. The clr mean was plotted for each RP-HPLC peak across the different habitat (x-axis) and (y-axis) populations for B. atrox snakes. Bars indicate the SE, a solid line indicates a perfect agreement, dashed lines indicate the origin (i.e., the geometric mean), and proteins less than these values were considered low-expression proteins.
In our PERMANOVA analyses, only in the spatially proximate degraded and forest habitats the venom composition did not significantly differ (Table 1). Accordingly, minor differences in both rare and abundant proteins were observed between these populations, which showed only slight differences in the rare proteins. However, in the populations from the other habitats, the venom differentiation was not constrained to rare proteins. On the contrary, in some cases, these proteins showed more similarity across populations from different habitats (for example floodplain vs. forest; pasture vs. forest) than proteins Toxins 2021, 13, 814 7 of 18 classified as abundant. In fact, except for peak 23 (rich in PIII-class SVMPs), abundant proteins showed differences across all habitats. Interestingly, the abundant proteins would be closely related to the functional activity of the venom toward their prey. As shown before [15], the conserved fraction 23 contains Batroxrhagin, a multifunctional PIII-class SVMP extremely conserved in different samples of B. atrox venom [2]. However, in Figure 4, we show that other abundant toxins are differentially expressed between the snakes from different habitats: fraction 21, which predominantly contains a PI-class SVMP, had lower expression in floodplain venoms; fraction 10, which contains mostly SVSPs, had higher expression in floodplain venoms; and fraction 20, which predominantly contains a PI-class SVMP and CTLs, was relatively over-expressed in venoms from pasture compared to the other areas.

Specific-Level Differentiation among Venoms of Snakes Collected in the Same Environment
Once the variability among the groups was defined, we quantified intrapopulation variability in venom composition, among specimens collected in close proximity in the same areas. As shown in Figure 1, HPLC chromatograms show variable profiles within venoms from snakes collected in the same environment (for more details see Supplementary Figure S1). Consistent with this pattern, we found great heterogeneity in the percentage area of each peak in the groups ( Table 2).
In every habitat, individual venoms showed important differences in both their percent area and the presence/absence of some peaks. Venoms from pasture and floodplain snakes were again the most distinct among the individuals from the same group (Table 2; Supplementary Figure S1), with venoms from the floodplain being the most heterogeneous. From Table 2 we emphasize three regions: In the PLA 2 -eluting fractions (Peaks 3 and 5), Peak 3 was absent in several chromatograms, but its absence is offset by the homogeneous increase of Peak 5 in pasture specimens. Venoms from floodplain specimens showed the most heterogeneous distribution of both PLA 2 fractions. CTL fractions (Peaks 15-17) are expressed at low but homogeneous levels in pasture venoms with increased variation in venoms from the forest and degraded area, reaching the highest heterogeneity in venoms from the floodplain. A notable observation is the high level of variation in the expression of the peaks containing SVMPs (peaks from 21 to 23) in floodplain specimens. These are the most abundant and homogeneous fractions in venoms from the other areas, but in venoms from the floodplain, fraction 21 is very low or even absent in venoms of five specimens while fraction 22, almost not present in venoms from other areas, is detectable in high levels in venoms from four specimens from the floodplain, indicating a high variety of SVMP isoforms venoms from snakes in this environment. It is important to note that the higher variability observed in venoms from floodplain snakes is not due to the two distant floodplain spots of snake collection. The differences appointed here can be observed within the snakes collected at Santarém (ATXV 5, 7, 8, 9, and 16) or Oriximiná (ATXV 10, 11, 12, and 13), which are approximately 300 km apart.

Differences in the Composition of Individual Venoms Resulting in Functional Variability
The most variable components among the venoms included SVMPs, SVSPs, and PLA 2 s. Our next step was to evaluate some of the main biological activities related to these protein families. With the goal of reducing the number of experimental animals for toxicity tests for ethical reasons and due to the limited amounts of some individual venom samples, venoms from only 16 specimens were evaluated per functional test, comprising four from each habitat. As shown in Figure 5, individual variation was observed in the catalytic activities of the major enzymatic components (SVMPs, SVSPs, and PLA 2 s) among the venom samples from each habitat ( Figure 5A,C,E), except in the venoms from the floodplain, in which the SVMP catalytic activity was similar and low. There were significant differences among the venoms collected at the same habitat in hemorrhagic ( Figure 5B) and myotoxic ( Figure 5F) activities. The four venoms from the floodplain snakes (V5, V8, V13, and V16) induced hemorrhagic spots comparable to those induced by snake venoms from the other habitats, with an emphasis on the V8 snake, whose venom had the highest hemorrhagic activity in the floodplain group. For myotoxic activity, the greatest variation was found among the venoms from floodplain snakes: this group was the only one in which all tested venoms showed statistically significant differences in activity ( Figure 5F). same habitat. Venoms from five specimens from each habitat were injected into groups of six mice or six chicks, and the time of death of the animals was monitored over 48 h ( Figure  6).  The greatest within-group difference was found in procoagulant activity among the floodplain snakes, in which the DC 50 values varied from 0.0035 to 1.965 µg for V8 and V10 snakes, respectively ( Figure 5D). No other group showed such a large variation in DC 50 values. However, the wide range of activity in this group was mostly due to the very high DC 50 value of only one specimen (V10). On the other hand, pasture venoms were similar in procoagulant activity with only a small amount of variation observed in the venoms from the forest and degraded areas.

Individual Heterogeneity in Venom-Induced Lethality
Finally, to link functional variation to toxicity to specific prey, we decided to investigate the ability of the venoms from each habitat to kill rodents (mice) or birds (chicks). We also quantify individual variation in venom-induced lethality within the same habitat. Venoms from five specimens from each habitat were injected into groups of six mice or six chicks, and the time of death of the animals was monitored over 48 h ( Figure 6). We observed a great deal of variation in the lethal activity of venoms from individual venoms from all groups analyzed. In venoms from the forest, degraded area, and pasture areas, mortality of experimental animals was observed from 1 to 2 h after venom injection, while in a number of cases, several animals survived through 48 h of observation. Similar results were observed for both the mice and birds, which matches the generalist diet preferences of B. atrox. There were some quite interesting patterns of toxicity in floodplain We observed a great deal of variation in the lethal activity of venoms from individual venoms from all groups analyzed. In venoms from the forest, degraded area, and pasture areas, mortality of experimental animals was observed from 1 to 2 h after venom injection, while in a number of cases, several animals survived through 48 h of observation. Similar results were observed for both the mice and birds, which matches the generalist diet preferences of B. atrox. There were some quite interesting patterns of toxicity in floodplain venoms. These venoms had similar activity to venoms from the forest, pasture, and degraded area when injected in mice. However, some venoms from the floodplain habitat (V8 and V16) induced deaths only a few minutes after its inoculation in chicks, a pattern that was not observed in mice. In addition, one venom (V12) was noticeably more lethal to chicks, killing half the animals in the first two hours and all animals of the group within 24 h, whereas in mice, the same venom induced the death of only two animals, 12 and 24 h after inoculation. However, not all floodplain venoms were highly lethal to birds. The V10 snake venom was much more lethal to rodents than to birds, killing all mice within 24 h, and only two chicks until the end of the experiment. The venom of snake V13 killed the same number of chicks and mice (only 2), although at different time intervals. Similar to the variations observed in the venom composition, the differential lethality to chicks observed with floodplain venoms was not due to the two distant floodplain spots of snake collection. Venoms that preferentially killed chicks were from snakes collected at Santarém (ATXV 8 and 16) or Oriximiná (ATXV 12 and 13).

Discussion
In this study, we analyzed the influence of different environments on the individual variability of venom samples obtained from B. atrox snakes captured in four different habitats of the Brazilian Amazon: forest, pasture, degraded area, and floodplain. Our analyses revealed clear differences in both venom composition and its biological activities among snakes from different habitats. However, major differences were also observed among venom samples collected from snakes in the floodplain habitat, a dynamic environment, subject to periods of annual drought and floods.
In a previous report involving B. atrox snakes, the variability in venom composition was attributed mainly to low-abundance proteins that would be a genetic reservoir for quick adaptive changes [2], while the most abundant isoforms from each toxin family were considered as "core toxins", conserved in venoms of most specimens and responsible for the major activities of the venom [2]. A similar observation was reported by Margres and collaborators [23], showing that, in Crotalus adamanteus, Sistrurus miliarius, and Agkistrodon piscivorus snakes, the differences in the expression levels were present mostly in the low-expression proteins. However, in the present study, both rare (low expression) and abundant (high expression) proteins show contribute to venom variation within and between snakes from different habitats.
For example, our data also showed that Peak 23 was composed mostly of Batroxrhagin, a P-III class SVMP from B. atrox venom isolated by our group [24], was widely conserved in all the 37 venoms analyzed, within and across all groups; this confirms a previous observation [2] that Batroxrhagin would act as a "core toxin" in B. atrox venoms, acting on important physiological targets of diverse prey types. In contrast, other abundant venom proteins showed less conservation and could represent more specific or even "adaptive" variation in these proteins. Clear examples of abundant fractions differentially expressed between the groups are fraction 21, which contains predominantly a PI-class SVMP and had lower expression in floodplain venoms; fraction 10, which contains mostly SVSPs and had higher expression in floodplain venoms; and fraction 20, which contains predominantly a PI-class SVMP and CTLs and is higher expressed in venoms from pasture compared to the other areas. Other fractions rich in SVSPs and CLTs were also were variable. These proteins participate to the hemorrhagic and coagulant activity of B. atrox venom and the observed differences in expression of the isoforms may alter procoagulant activity and to promote differential lethality in rodents and birds, suggesting an important role of the environment in selecting for venom variation.
Individual variation in the abundance of SVMPs, SVSPs, and PLA 2 s could have arisen as a result of environmental differences between habitats. In our previous study [15], we compared the functional profiles of pools of venoms from B. atrox snakes from the same habitats and observed a less hemorrhagic and more procoagulant phenotype in the pool of venoms from floodplain snakes. The pool of venoms from the floodplain snakes was also significantly faster to induce clotting in several types of plasmas, including avian plasma [25]. Here, we used similar experimental approaches to investigate the individual variability within each group of snakes, introducing modifications in the tests of DC 50 , which were performed in the presence of calcium, to ensure the detection of coagulotoxins, dependent or not on this cofactor. In addition, our tests were performed with avian plasma as previously reported [26,27], which allowed us to construct better dose-response curves to assess the procoagulant effects of the individual venoms.
Venoms from B. atrox specimens from the floodplain showed low SVMP activity and higher SVSP catalytic activities. However, the hemorrhagic and procoagulant phenotypes observed in the pool of venoms were not consistent across all individual venoms collected from this environment. Some individual venoms such as from V5 and V8 snakes, showed both a potent procoagulant action on avian plasma and high hemorrhagic activities in mice. Moreover, the higher DC 50 value observed, which corresponds to the less procoagulant venom, was precisely from a floodplain snake venom (V10), demonstrating that not all snakes from the floodplain have venoms that induce potent clotting of avian plasma, but the range of clotting activity was the highest in this environment.
Bernardoni and collaborators [3] analyzed Bothrops neuwiedi venom and found SVMPs isoforms that are functionally different and capable of affecting different targets in the hemostatic systems of birds, rodents, and humans. They suggested that some SVMPs are less selective, guaranteeing the action of the venom on different targets, while other isoforms are more selective, modulating the action of the venom for specific prey. More recently, we also investigated the effects of these SVMP isoforms on amphibian plasma (Rinella marina), demonstrating the coagulotoxic effects of these toxins on different types of animal plasmas [18]. In the present study, some fractions rich in SVMPs had a very variable distribution in the B. atrox venoms of all habitats and may be involved in the variability observed in the coagulotoxic activity observed in these individual venoms. In contrast, the major SVMP eluted in Peak 23 was very conserved in the venoms of all 37 specimens analyzed, confirming the previous assumption [2] that Batroxhagin is the core SVMP of B. atrox venom. On the other hand, the strong hemorrhagic action observed in some floodplain venoms (V5, V8, and V13), is compatible with a prominent increase of the fractions represent by Peaks 21 and 22 that shows height/area comparable to the ones observed in venoms from the other habitats. These peaks contain Atroxlisin-Ia, a PI-class SVMP that induces hemorrhage comparable to the P-III class enzymes [28].
During envenomation, various components present in snake venoms can act synergistically to cause the prey's organism to collapse, which usually results in rapid death. The speed to kill/immobilize prey is crucial for food acquisition of terrestrial viperids during hunting [29,30]. Previous works also report that lethality profiles of snake venoms can be variable in different animal models, such as mammals, birds, reptiles, and amphibians [5,31]. For this reason, we included the functional characterization of the individual venoms in two different animal models, birds (chicks) and mammals (mice). We found substantial differences in venom activities in both, the number of test animals dying, and the time of deaths in the two models. Forest venoms showed similar lethality profiles in rodents and birds across individual venom samples, while venoms from pasture and degraded areas were slightly more lethal to rodents. Of particular note, venoms from floodplain snakes were highly variable according to their ability to prey on mammals or birds. For example, the venoms of the snakes V8 and V16 were more toxic to chicks, inducing deaths within few minutes in birds, even more quickly than the venom of B. insularis, a snake with a diet specialized in birds [32]. The venom of the snake V10 however was more lethal to rodents. Interestingly, the venoms with the highest procoagulant activity are more lethal to birds while V10 venom, the weakest procoagulant, killed predominantly mice. The limitations imposed by the small available quantities of venom prevented using the same individual venoms in all the functional tests performed or to carry out doseresponse analyses, hindering a more direct assessment of the relationship between lethality and the other toxic activities evaluated.
Nevertheless, our findings strongly suggest that the procoagulant activity of the floodplain venoms contributes significantly to their lethal effects in rodents and birds. The most striking and differential characteristics of the floodplain venoms are its potent effects on coagulation in different types of plasma. The venom of the snake V10 (from floodplain habitat) was more efficient in killing rodents than birds, and this same venom had the highest EC 50 value on avian plasma. A possible explanation for the remarkable action of some floodplain venoms on hemostasis and velocity to induce death in birds could be directly linked to the floodplain environment, where these snakes were captured. Floodplains are periodically flooded by the lateral overflow of waters rich in sediments from the Amazon River. The floodplain regions alternate cyclical periods of drought and flood, remaining flooded for a few months during the year, this seasonal variation being driven by flood and precipitation seasons [33]. Balancing selection pressures could result in greater venom heterogeneity, as evidenced by the differential lethality to rodents or birds observed in some floodplain venoms, while the snake venoms from the other habitats presented a more homogeneous lethality pattern.
In line with our findings, Smiley-Walters et al. [31] showed that variation in the venom of Sistrurus miliarius at the population level has an adaptive role in terms of toxicity to prey. Testing adaptive hypotheses requires careful analysis of phenotypic characters whose variation has clear functional consequences. Venoms directly affect the ability of an individual snake to immobilize and kill its prey [31,[34][35][36][37]. For B. atrox, a snake with a generalist diet that includes arthropods, frogs, lizards, birds, and small mammals [38,39], the functional diversity of the venoms could represent a local adaptation of individual venoms in each population to different sets of prey. Within an ecological context, the nature of the interactions between species, including prey-predator relationships, may change among populations [40,41], especially for species with a wide geographical distribution. of the functional versatility and diversity of B. atrox venoms may explain its wide distribution throughout all Amazon regions [7], reflectinghow this species adapts to prey communities in different habitats.

Conclusions
Individual venom variation of B. atrox snakes, captured in different habitats within the Brazilian Amazon, support the idea that dynamic environments may select for more variable venoms. Such differential distribution of protein isoforms leads to differential function and toxic effects on different prey. The fact that the greater heterogeneity in terms of composition/toxicity of the venom was found in snakes from the floodplain habitat suggests that balancing selection for expression of different isoforms could be enacted by the drastic seasonal changes (drought/flood) present in this type of environment. Variable selective pressures in different habitats may likely exert localized impacts on the venom phenotype in a species with a wide geographic distribution, such as B. atrox.

Snakes and Venoms
Thirty-seven male and female adult specimens of B. atrox snakes, with sizes ranging from 71.2 to 124.5 cm, were captured in four environments at Santarem and Oriximiná, in the western region of the state of Pará, Brazil (Supplementary Figure S2) under ICM-Bio/SISBio license 32098-1 and SISGEN number A78BD88: (1) Pasture: ten adult snakes were collected from a pasture area in the municipality of Oriximiná ( ATXO 1,2,3,5,6,7,9,15,16,19) on the north shore of the Amazon River. This site was historically Terra Firme forest characterized by large trees that were cleared for pasture~20 years earlier.
(2) Floodplain (ATXV): five adults were collected from a seasonally flooded island in the main course of the Amazon river near Santarem (at Urucurituba: ATXV 5,7,8,and 9; at Tapará: ATXV 16) or Oriximiná (ATXV 10,11,12,13). These areas are typical of floodplain habitats subject to periodic flooding by the Amazon River and formed by the deposition of sediments that have led to the formation of many islands. The typical vegetation consists of grasses that grow on highly fertile alluvial soils. (3) Forest: ten snakes were captured at the Floresta Nacional do Tapajós, a protected area located in the municipality of Belterra (ATXF 24,25,26,28,29,30,31,33,34,and 35) next to the Tapajos River about 50 km south of the Amazon River. This site also represents the upland Terra Firme forest. (4) Degraded: eight snakes were collected at a recently degraded (ATXD 3, 4, 6, 7, 8, 9, and 10) area contiguous to the forest area at Floresta Nacional do Tapajós, which was cleared for pasture. After capture, the snakes were transferred to the Herpetarium of Laboratório de Pesquisas Zoológicas, Universidade da Amazônia (UNAMA), in Santarém, PA. Venom samples were collected using manual extraction techniques, freeze-dried, and kept at −20 • C until use. Animal care and procedures used in the handling of snakes were undertaken according to the guidelines and permits (CEUAIB 1244/14, Instituto Butantan, São Paulo, Brazil).

Chromatographic Analysis
The B. atrox venom samples were individually fractionated by reverse-phase highperformance liquid chromatography (RF-HPLC) in the Vydac C18 column as previously described [15]. The fractions had their toxin composition predicted according to peak shape and retention time comparing to a previous standard chromatography of B. atrox venom from snakes collected in the same areas from which components present in each fraction were identified and quantified by free-label mass spectrometry [15].

Enzymatic Assays on Synthetic Substrates
SVMP, SVSP, and PLA 2 enzymatic activities of individual venoms from B. atrox were assayed using synthetic substrates according to the procedures previously standardized in our lab [42]. Briefly, the substrate used for SVMPs was Abz-AGLA-EDDnp substrate (Peptide International, Gardner, MA, USA), and the results are expressed as Relative Fluorescence Units-RFU/min/µg. The PLA 2 activity was assayed using the substrate 4-nitro-3-[octanoyloxy] benzoic acid (Enzo Life Sciences, New York, NY, USA) and activity were determined according to the absorbance at 425 nm and expressed as Absorbance/min/µg of venom. For SVSP catalytic activity, the chromogenic synthetic substrate benzoyl-arginylp-nitroanilide (L-BAPNA) (Sigma-Aldrich ® , St. Louis, MO, USA) was used and the hydrolysis of the substrate was expressed as the increase in Absorbance/min/µg of venom. The results represent the mean ± SD of three independent experiments, undertaken in triplicate.

Procoagulant Activity
The procoagulant activity was evaluated in recalcified chicken (White leghorn) plasma, obtained under license CEUAIB n • 13,710-14, Instituto Butantan, by thromboelastometry using a four-channel ROTEM ® system (Pentapharm, Munich, Germany), as previously described [26]. The plasma was prepared as 3.2% citrated stock and stored at −80 • C until needed in 1 mL aliquots. For the experiments, an aliquot was defrosted by placing it into a 37 • C water bath for 10 min. Venom samples and controls were dissolved in PBS (60 µL) and added in the specific cups with 20 µL of CaCl 2 (0.02 M) and 260 µL of citrated plasma (final volume = 340 µL), pipette-mixed, and the clotting time was immediately recorded. PBS and ellagic acid were used as negative and positive controls, respectively. Under such conditions, recalcified chicken plasma treated with PBS only is unclottable for at least 30 min (CT value = 1800 s). The results are shown as Coagulant Dose 50% (DC 50 ), defined as the amount of venom (µg) that reduces the CT parameter to 900 s.

Venom Activities Using Animal Models
Hemorrhagic and myotoxic activities were carried out as previously described [43] using male Swiss mice (18-20 g) under the approval of the Butantan Institute Ethics Committee on Animal Use (Protocol Number: 13710-14). Briefly, for hemorrhagic activity, samples containing 10 µg of each venom pool, diluted in 50 µL of PBS, were injected intradermally into the dorsal skin of mice. After 3 h, the animals were euthanized in a CO 2 chamber, the dorsal skin was removed, and the hemorrhagic halos were measured. Groups of 5 animals were tested, and control group animals were injected with PBS only or B. jararaca venom (10 µg). Results represent the mean ± SE of the area of hemorrhagic spots (cm 2 ) from mice tested in at least 2 independent experiments. The myotoxic activity was assayed using 50 µg of venom pools in 30 µL of PBS, injected intramuscularly into the gastrocnemius muscle in Swiss mice. After 3 h, the animals were bled via ophthalmic plexus and the sera were assayed for creatine-kinase activity with a commercial kit CK-UV (Bioclin), according to the manufacturer's instructions. Groups of 5 animals were tested, and animals from the control groups were injected with PBS only or B. jararacussu venom (50 µg). Results represent the mean ± SE of the CK activity in mice sera (U/mL) from mice tested in at least 2 independent experiments.
The lethality induced by the venoms was tested in murine (Swiss mice 18-20 g body weight) and avian (three days old Bovan chicks 20-30 g body weight) models with protocols approved by the Animal Ethical Use Committees of the Instituto Butantan (CEUAIB n • 13710-14). Groups of 6 mice or chicks were injected intraperitoneally with a single dose (200 µg) of venom samples, diluted in 500 µL of PBS. Following injection, the times of death for each group were recorded hourly until the sixth hour, and then at 12, 24, and 48 h after inoculation of the samples. After 48 h, the surviving animals were euthanized, using a CO 2 chamber for the mice and injection of sodium pentobarbital overdose for the chicks. The data obtained were graphed on a survival curve plot. The tests were undertaken in two independent experiments.

Statistical Analysis
We tested for significant compositional differentiation using a hierarchical multivariate approach that is robust to the compositional nature of the data [44]. First, we brought the relative abundance data out of the compositional simplex space using the isometric log-ratio (ilr) transformation, which also avoids the zero-sum constraint of the centered-log ratio (clr) transformations but at the cost of one dimension of the data. The ilr-transformed peak data were then subjected to permutational multivariate analysis of variance based on distances using the adonis function from the vegan R package [45], with all peak dimensions as a response and environment as a predictor. Significance and R 2 were estimated using 10,000 permutations of the raw data. We then visualized the venom phenotype of each snake in multivariate space using robust principal components analysis (PCA) as implemented in the pcaCoDA function from the robCompositions R package [46].
After detection of a significant global association between environment and venom composition, we conducted pairwise posthoc PERMANOVA, as described above, between the venom compositions of all pairs of environments to determine which environments differed significantly. Lastly, we preserved the full dimensionality of the data but removed the data from the simplex using the clr-transformation, to determine which peaks specifically varied across environments. To do this, we use the lm function from R stats, with the clr-transformed abundance of an individual peak as the response and environment of origin as the predictor. We repeated this test for all 26 peaks and used the false-discovery rate (fdr) approach to correct the p-values for multiple tests.
For the functional analyses (in vitro and in vivo assays), firstly, the data were evaluated for a normal distribution (Shapiro-Wilk), and then differences among the means were evaluated by one-way ANOVA, followed by a Tukey post-test (for multiple comparisons). Data that did not meet the normality criteria were analyzed using a non-parametric test (Kruskal-Wallis). Results represent the mean and standard deviation or standard error, as appropriate, and the level of significance was set at p ≤ 0.05. Data were analyzed using the GraphPad Prism statistical program (version 7.00 for Windows, GraphPad Software, San Diego, CA, USA).

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/toxins13110814/s1, Figure S1: Chromatographic profiles of B. atrox venoms of different habitats from west of Pará State, Brazilian Amazon; Figure S2: Characteristics and localization of areas of snake collection.