*2.1. Rarity Patterns*

We used literature and unpublished information to characterize the geographical range (GR), habitat breadth (HB), and local abundance (LA) of 143 species of NW pitvipers (see Tables S1–S3). We follow the taxonomy of the Reptile Database as of August, 2015 [28], except for some island populations of *Crotalus* (see [43,44]). Each variable considered was dichotomized (see details below). As a result, we built a table with the eight combinations of these variables in the table cells (cf. Rabinowitz [10]). We then categorized species in four rarity categories (see [14]): Not rare (NR), which includes those species that are not rare in any of the three aspects of rarity; low intermediate (LI), with species that are rare in one aspect; high intermediate (HI), with species that are rare in two aspects; and rarest (RT), with species that are rare in all three aspects.

Locality data used to calculate GR of each species were gathered mostly from literature (scientific papers, theses, dissertations, books, reports etc.), but also from museum databases (including VertNet [45]), and from an unpublished database for species that occur in Brazil [46] (see Table S1). Each resulting map was compared to the maps provided by Campbell and Lamar [29] to exclude possible erroneous localities. In a few cases in which locality data were not available, we obtained data directly from published distribution maps.

We estimated GR using a map depicting the point occurrences of a given species and drawing the smallest convex polygon connecting external points and that encompasses the remaining points (cf. extent of occurrence, EOO [47]). We then excluded areas that were clearly unsuitable habitat (using literature and unpublished information on degree of habitat specialization) and that represented over 1% of the total polygon. In order to make very small distributions more realistic, we considered a bu ffer of five km radius centered at each point of occurrence. For island endemics, we used the area of the island as the species' GR. For most species, the area of our GR was the same as the EOO [47], but in some, the GR obtained was up to one fourth the EOO (e.g., *Bothrops asper*; see Table S4). We then dichotomized the distribution of GRs using its median: GR > 82,900 km<sup>2</sup> = wide distribution; GR < 82,900 km<sup>2</sup> = restricted distribution). All spatial analyses were performed in the software Quantum GIS [48].

We estimated HB of each species using literature data to obtain the number of major habitat types in which the species can be found, considering gross vegetation types: Desert, non-desert open habitats (grassland, shrubland, and savanna), and forests (Table S2). We then dichotomized the distribution of number of habitat types using its median: > one habitat type = wide habitat breadth; one habitat type = narrow habitat breadth (Table S2). To estimate LA of each species, we used a simplification of the ACFOR (abundant, common, frequent, occasional, rare) scale [49] with three categories (abundant, common, and rare), which refer only to the apparent abundance found in the literature (see Table S3). For each species we assigned a value from one to three (rare to abundant). For dichotomization, we used the mean of these values (1.98; one = rare; two and three = not rare) instead of its median (2), as this would make dichotomization impossible. We failed to find LA data for 16 species. Thus, of the 143 species of New World pitvipers considered herein, we had enough data to analyze rarity for 127 species (88% of the species).

We analyzed the associations between the three rarity variables (GR, HB, and LA) using log-linear models [50,51]. We fitted alternative log-models to the frequencies of species in each combination of the three rarity variables. The models expressed all combinations of associations among variables in a multi-way contingency table. We started with a full model (all possible associations between two variables plus the association of the three variables) and then we tested sequentially the loss of fit caused by the removal of each one of these associations [51]. The loss of fit in these sequential comparisons was gauged by the log-likelihood ratio test (LRT) with the Chi-square test. A significant loss of fit when an association was removed was considered support for such association. The first step in this forward selection procedure was to compare through LRT the model with all interactions

with the model without the three-way interaction. If the test showed a non-significant loss of fit, we inferred that the data did not support the three-way association between the three forms of rarity. In this case, we proceeded by testing the significance of the three two-way associations (that is, the association between GR and HB, GR and LA, HB and LA). In each of these tests we compared through LRT the model with the target two-way interaction (but with the other remaining interactions) with the model with all two-way interactions. Each of these three tests then allowed us to identify if there was evidence for one of the associations among pairs of rarity forms (for instance, if species with small geographic range and small populations were more frequent than expected if these two states were independent, and so on). We performed these analyses in R 3.1.3 (R Core Team, Vienna, Austria) [52].

## *2.2. Predicting Rarity*

We explored if body size and latitude can predict rarity in NW pitvipers. Body size is frequently related to several life history aspects in animals [53], and might affect abundance (see [54]) and home range size (e.g., [55–57]). In snakes, body size can also be tied to aspects often predicting generation time, such as time of sexual maturity, fecundity, and lifespan (e.g., [58,59]). Thus, body size could potentially predict rarity in snakes. We collected data for body size (total length, i.e., snout-vent length plus tail length) through an extensive literature review (see Table S5).

Latitude could also predict rarity patterns among NW pitvipers. According to Rapoport's rule, in temperate zones species tend to have larger geographical ranges than in the tropics [60–64] and species occurring in tropical regions tend to have narrower habitat breadth than those in temperate zones [65–67]. Additionally, there is evidence that species abundance is positively related to latitude [68,69]. To evaluate if latitude predicts rarity patterns across NW pitvipers, we used the latitude of the centroid of the geographical range of each species.

To investigate if body size and latitude predict rarity, we used Bayesian phylogenetic generalized mixed models implemented in the R package MCMCglmm [70]. Under this approach phylogenetic relatedness is incorporated as a covariance matrix. We followed the approach of Verde Arregoitia et al. [71] and used body size and latitude as fixed effects and the rarity ranking categories as an ordinal dependent variable. We assigned a rank ranging from 1 to 4 to our rarity categories (cf. [14]), in which the lower value represents the rarest category and the higher value indicates the not rare category. We implemented our models using the codes by Verde Arregoitia et al. [71] and the phylogeny generated in Alencar et al. [31]. We ran the analyses for at least 10,000,000 steps to achieve convergence, discarding the first 10% iterations as burn-in. We verified convergence using the Heidelberger and Welch's convergence diagnostic [72], checking for the autocorrelation of successive samples and through visual inspection of trace plots. We considered associations as significant when the 95% credible intervals did not include zero.

#### *2.3. Rarity and Other Assessments of Extinction Risk or Prioritization*

Categorizations under the IUCN criteria [47] were gathered from the IUCN Red List [42], from assessments by the IUCN Snake and Lizard Red List Authority, and from the recent assessment of Brazilian snakes for species endemic to Brazil (see [43]). We assessed the remaining species for which previous IUCN categorizations were not available using the IUCN categories and criteria [47] (see Table S4). Continuing declines in habitat quality under criterion B were inferred visually using maps of human influence index [73], superimposed on distribution maps, combined with information on the ability of species to persist in disturbed habitats (obtained in the literature or by specialist consultation).

We also compared our rarity categorization of New World pitvipers of this study (only the rarest, RT category) with those of the Mexican red list [74], the Threat Index of Maritz et al. ([43] for the World vipers), and conservation status assessments using the Environmental Vulnerability Score (EVS; [75]) for Central American [76] and Mexican species [77].
