Co-occurrence patterns between false coral snake Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae) and venomous coral snakes from the Amazon

: Batesian mimicry may result in remarkable cases of phenotypic convergence that represent classic examples of evolution through natural selection. The existence of mimicry systems among coral snakes, however, remains controversial because of contradictions between the predictions of mimetic theory and the empirical patterns of co-occurrence and species abundance. Here, we analyze the geographic distribution of coral snake species of the genus Micrurus and populations of the false coral snake Atractus latifrons in Amazonia, and perform ecological niche modeling (ENM) analyzes to generate potential geographic distributions of species of Micrurus and A. latifrons , identify patterns of co-occurrence and assess whether the distribution of A. latifrons coincides with the distribution of Micrurus species, which could suggest the existence of a possible mimetic relationship between the species. We identifi ed six Micrurus species that may represent mimetic models for A. latifrons. The results of the co-occurrence analysis corroborates the results from ENM, indicating that chromatic patterns of A. latifrons and their respective model species are aggregated. Our study suggests that all color patterns of A. latifrons – including the tricolor monads, and the more common tricolor dyads and tricolor tetrads – may benefi t from the resemblance with other Micrurus species as perfect and imperfect mimics.


INTRODUCTION
The selection of adaptive similarities involving mimics and models represents a classic problem in evolutionary biology (Kikuchi & Pfennig 2012). The term Batesian mimicry refers to a phenomenon in which a nonnoxious species replicates a visual, acoustic or chemical signal emitted by a venomous species (model), resulting in adaptive advantages against predation (Ruxton et al. 2004, Hossie & Sherratt 2013, Rabosky et al. 2016. In such systems -Batesian mimicry complexes -Ruxton et al. (2004) claim that imitators should occur only in sympatry with their models. However, Pfennig & Mullen (2010) mentioned exceptions in which Batesian mimics exist in areas lacking their putative models, or where distributions of mimics are not totally included in the models' range. Nonetheless, the precise distributional relationships (full sympatry or allopatry) of potential models and their respective presumed mimics are largely unknown for many groups potentially involved in Batesian mimetic associations.
Several studies suggest that the cooccurrence of mimics and models may affect the abundance of perfect and imperfect mimics (Huheey 1976, Brodie III & Janzen 1995, Pfennig et al. 2001, Caley & Schluter 2003, Harper & Pfennig 2007, Kikuchi & Pfennig 2013. In areas where models are scarce or absent, the chances of predation on mimics would tend to increase (Huheey 1976, Pfennig et al. 2001, Harper & Pfennig 2007; alternatively, in areas where the model is relatively common, the odds of the predator finding it will be higher, promoting avoidance (i.e., lower predation rates) of any perfect or imperfect mimics (Brodie III & Janzen 1995, Caley & Schluter 2003, Kikuchi & Pfennig 2013.
Mimicry evolved in different snake lineages, with vipers (Viperidae) and coral snakes (Elapidae) frequently serving as models. Coral snakes advertise their danger with aposematic red, black, and white (or yellow) rings (Roze 1996). In the Neotropical region, coral snakes of the genera Micrurus, Leptomicrurus, and Micruroides are often associated with mimetic complexes involving non-venomous (aglyphous) or mildly-venomous (opisthoglyphous) species, known as false coral snakes, belonging mostly to the genera Atractus, Lampropeltis, Tantilla, Oxyrhopus, Erythrolamprus, Pliocercus and Simophis (Greene & McDiarmid 1981, Campbell & Lamar 1989, 2004, Savage & Slowinski 1992, Roze 1996, Almeida et al. 2016. The vast majority of assumptions about mimicry systems involving coral snakes rely on sympatric distribution and similarities in chromatic patterns between venomous and non-venomous species (Greene & McDiarmid 1981, Marques & Puorto 1991, Brodie III & Janzen 1995, Torre-Loranza et al. 2006, Almeida et al. 2014, Bosque et al. 2016. Combining phylogenetic techniques, distributional comparisons, and general similarity, Rabosky et al. (2016) demonstrated a high correlation between the occurrence of red and black band patterns (RBB, sensu Rabosky et al. 2016) of venomous and non-venomous coral snakes in space and time. Thus, RBB patterns would only have appeared in non-venomous coral snakes after its emergence in venomous coral snakes, and only in areas of sympatry between both forms. However, according to Brattstrom (1955), the color similarities between such species may have resulted from evolutionary convergence around a cryptic pattern that provides protection, while the fixation of the color pattern resulted from the selection of a disruptive effect.
The assessment of congruent geographic patterns among mimics and models may be limited by low rates of abundance, precluding trustworthy analyses about patterns of geographic overlap. In this sense, ecological niche modeling (ENM) has been broadly used in recent decades to infer the potential geographic distribution of species from incomplete empirical information about their distributions. ENM associates known occurrences of species throughout the geographic space with environmental conditions of evaluated sites (Rangel & Loyola 2012). From such associations, one can identify and map new locations where viable populations of species are likely to occur (Peterson et al. 2011). An important point of using ENMs to estimate geographic distributions is that they often predict broader ranges because they do not consider biotic interaction and dispersal constraints which limits the species distributions at local scales (Peterson et al. 2011, Uribe-Rivera et al. 2017. Despite this, if some aspects are taken into account in the modeling process, such as the accessible area for species occurrence, the quality of species occurrence records, and the niche-modeling method used (Jiménez-Valverde et al. 2008, Peterson et al. 2011, the distributional area estimated by ENM could provide a good approximation of the occupied distributional area at large scales (Peterson et al. 2011). Therefore, ENM represents an important tool for overcoming Wallacean shortfalls (Lomolino 2004, Hortal et al. 2015 and testing ecological and biogeographical hypotheses. In this study, we evaluated the congruence between the geographic distributions of different color patterns of the false coral snake A. latifrons and the true coral snakes of the genus Micrurus in Amazonia to suggest, based on the co-occurrence between the species and the similarity between the color patterns, the existence of a possible mimetic relationship between the tested species. In order to reduce the effects of sampling gaps in the accurate assessment of distributions, we use ENM to generate potential geographic distributions of Micrurus and A. latifrons species to test species correlation, assuming sympatry as evidence of the Batesian mimetic complexes. However, to perform the co-occurrence test of chromatic patterns of A. latifrons and Micrurus species, we recognize here the variations described as "tricolor dyads" defined by Almeida et al. (2014) as actually variations of the pattern "Tricolor monads", which we designate as "Tricolor monads 1" (similar to Micrurus averyi) and "Tricolor monads 2" (similar to Micrurus langsdorffi).

MATERIALS AND METHODS
To select Micrurus species acting as potential models for Atractus latifrons, we followed two parameters: (1) the potential model of a mimetic species should occur in sympatry with A. latifrons; (2) the potential model should exhibit a typical coral pattern (black, white and/ or red rings), forming monads or triads on the body and tail, as well as a black dorsal cephalic cap similar to the mimic patterns. The specieslevel identification of most of the Micrurus taxa followed Roze (1996), except the for M. lemniscatus species group, in which we followed Pires et al. (2014).

Ecological niche modeling and geographic distributions
We used ecological niche models to predict geographic distributions for species from known occurrence records. For niche modeling analyses, we used longitude and latitude data associated with specimens having clear and unproblematic occurrence records (Appendix S2). The models were generated considering the extent of the Neotropical region as the historically accessible area since both genera, Atractus and Micrurus, have a Neotropical distribution. We followed the recommendation by Araújo & New (2007) to generate consensus predictions by combining outputs from different niche modeling methods (weighted by individual model accuracy) and climatic models (Rangel & Loyola 2012). We obtained a set of bioclimatic variables at 0.5 o resolution from the ecoClimate database (http://ecoclimate.org, Lima-Ribeiro et al. 2015), representing pre-industrial climatic conditions, for five coupled atmosphere-ocean general circulation models (AOGCM): CCSM4, CNRM, GISS-E2-R, MIROC-ESM, and MRI-CGCM3. To minimize collinearity among bioclimatic variables, we used a varimax-rotated factor analysis to select the set of variables with the highest loadings in the first five factors (Terribile et al. 2012). Thus, we chose the following five orthogonal variables: annual mean temperature, annual temperature range, precipitation of the wettest month, precipitation of the driest month, and precipitation of the warmest quarter. We performed factor analysis based only on CCSM bioclimatic variables, as these variables were highly correlated across all AOGCMs (i.e., the loadings across all AOGCMs provided quite similar results, with shifts usually below 0.05 of difference in the highest loadings).
We used 11 modeling algorithms, including presence-only and presence-absence methods, ranging from simple bioclimatic envelope models (BIOCLIM) and distance-based methods (Euclidian distance and Gower distance), multivariate methods (ecological niche factor analysis -ENFA and Mahalanobis distance), as well as more complex methods such as regression models (Generalized linear model -GLM, Multivariate adaptive regression splines -MARS and Generalized additive model) and machine learning-based methods (Maxent, Artificial Neural Networks, and Randon forest). For a general description of niche models, see Franklin (2009) andPeterson et al. (2011).
For the modeling process, we mapped occurrence records at the same resolution as climate data. To evaluate the congruence between model and mimic species distribution, we modeled each chromatic pattern of Atractus latifrons separately ( Figure S1, Table SI). The species presence records (and pseudoabsences for those methods that use presence and pseudo-absence) were allocated randomly into 75% for calibration and 25% for evaluation, and we repeated this process 50 times. Pseudo-absences were randomly selected in the background region (excluding cells with occurrence records) with the same proportion as species records (prevalence of 0.5). The 50 repetitions in each method were converted into presence-absence maps based on thresholds established by the area under the ROC curve (AUC). All 300 models (50 repetitions x 6 modeling methods) were included in the consensus maps of each AOGCM, weighted by their model fit according to the True Skill Statistics (TSS, Allouche et al. 2006, Diniz-Filho et al. 2009). The average across the ensemble outputs of each AOGCM resulted in the final consensus map of habitat suitability (varying from 0 to 1) for each Micrurus species and chromatic pattern of A. latifrons. Final accuracy of ensemble models is presented in the Appendix S2. We generated all models using the computational platform BioEnsembles (Diniz-Filho et al. 2009, Terribile et al. 2012. We used habitat suitability predictions to map the potential geographic distribution of each Micrurus species and chromatic patterns of A. latifrons. Because ENM analyses do not consider biotic interactions and other local factors limiting species ranges (e.g. geographic barriers), their transferability (extrapolation in geographic space) may include areas that are inaccessible for the species. Since we are dealing with species with poor dispersal, we adopted conservative parameters to delimitate more realistic, easily interpreted distributions. Thus, we applied a decision threshold of 0.5, and considered areas (or cells) suitable if they had suitability values higher than 0.5. We assessed the final maps based on expert background knowledge (P.A. and D.T.) to ensure the predictions reflected the real potential distribution of each form.

Testing for co-occurrence patterns
From the geographic distributions obtained from ENM, we tested for non-random patterns of co-occurrence as evidence of positive relationships between Micrurus species and A. latifrons chromatic population coexistence. We performed the co-occurrence analysis using the co-occurrence module of EcoSim (Gotelli & Entsminger 2004). We used the matrix of presence (1) and absence (0) generated by ENM, and set each species and chromatic pattern in rows and the localities (or cells) in columns. We used the C-score index (number of "checkboard units") (Stone & Roberts 1990) as a quantitative co-occurrence index (the higher the C-score, the less co-occurrence, on average, between all of the species' pairs in the matrix). We estimate C-score through 5,000 simulations by reshuffling elements within each row of the matrix to test the probability that each species pair is significantly aggregated (observed C-score lower than simulations) or segregated (observed C-score higher than simulations). We also generate polygonal geographic distributions for each Micrurus model and mimetic Atractus by generating a convex hull around the occurrence locations ( Figure S2), while removing raster cells within these polygons that the ENMs identified as unsuitable. Then, we repeated the co-occurrence analysis with these polygons and compare it with that from ENM to evaluate if the results could be influenced by overprediction in ENM distributions.
To assess a positive relationship between the distributions of mimics and models, we replicated the procedure of niche modeling and included the vector of habitat suitability of each model species of Micrurus as a predictor variable (together with the five climatic predictors described above) in the modeling process of each corresponding Atractus latifrons mimic chromatic pattern. Then, we adjusted a multiple regression model between habitat suitability of each chromatic pattern of A. latifrons (the dependent response variable) and the set of six independent variables (habitat suitability of the corresponding Micrurus sp. Model plus the five climatic variables), and noted the standardized regression coefficients for each independent variable. We used these standardized regression coefficients to evaluate the correlation between habitat suitability (i.e., the geographic distribution) of each mimic with the habitat suitability of its corresponding model. Although we are aware that the significance of these coefficients cannot be considered due to the implicit collinearity between suitability and the climatic variables used in the ENM process, this approach was useful to assess a possible association of the geographic distribution of mimics with the distribution of the corresponding model (each standardized coefficient has partial importance regarding the other predictor variables in the regression model).
Finally, to evaluate if the geographic distribution of mimics could be constrained by the distribution of Micrurus models, we modeled the geographic distribution of A. latifrons as a whole (i.e., using all records of all mimics) following the same protocol as above. From this modeled distribution (see Figure S3), we re-sampled 1000 times geographic ranges (i.e. null ranges) for each chromatic pattern and calculated the proportion that each of these null ranges overlapped with the respective Micrurus model distribution. The size of these null ranges for each A. latifrons chromatic pattern was equal to the ENM distribution, and we also set for these null ranges to be continuous (thus approximating the modeled ranges). We then compare if the observed overlapping proportion for each model Micrurus and respective A. latifrons chromatic pattern was frequently higher (or not) than the proportion from the null ranges.

Geographical analyses
Considering only the chromatic patterns of mimics and models, we observed a strong association between venomous coral snakes and the different forms already described for Atractus latifrons. However, in addition to this chromatic association, the distribution analyzes also indicated the existence of sympatry between supposed models and mimics ( Figure 5).
The distributions generated by niche modeling showed a clear pattern of range overlap between the tricolor dyads (TD), tricolor tetrads (TT), tricolor monads (TM) and bicolor monads (BM) chromatic patterns of Atractus latifrons. The analysis of potential species distribution indicated (1)   and M. averyi in part of the Colombian Amazon, in the Peruvian Amazon and in the Brazilian Amazon, specifically in the states of Amazonas and Roraima ( Figure 5).
The results of the co-occurrence test for Atractus latifrons color pattern and their respective models indicated positive relationships (i.e., the probability of the observed C-score index ≤ the simulated index was 0.999) for all associations, indicating that chromatic patterns of A. latifrons and their respective model species are aggregated. The co-occurrence test corroborated the association involving patterns of bicolor monads, tricolor  monads, tricolor dyads, and tricolor tetrads, suggesting that the distributions of chromatic patterns of the mimic A. latifrons and of the different models of Micrurus do not occur randomly. The results were the same for the cooccurrence test from the polygons obtained by convex hulls, except for the pattern A. latifrons TM vs M. averyi, which was non-significant (i.e., Figure 5. Geographic distribution from ENM for chromatic patterns of Atractus latifrons and correspondent coralsnake mimetic model. the probability of the observed C-score index ≤ the simulated index was 0.389). This probably occurred because the convex hull for both A. latifrons TM and M. averyi were smaller than the distributions estimated by ENM and thus, the overlap between them was minimal ( Figure S2).
We found that the spatial distribution of A. latifrons exhibiting bicolor monads and M. albicinctus (measured as habitat suitability) had the highest standardized coefficient, indicating that the distribution of this chromatic pattern of A. latifrons has a strong correlation with the distribution of M. albicinctus (Table I; Figure 6), corroborating the co-occurrence test. For the chromatic tricolor dyads pattern, the models M. lemniscatus and M. diutius exhibited the lowest correlations. For the other chromatic patterns, correlations were high for the associations between the tricolor monads pattern vs. M. langsdorffi and intermediary correlations for the tricolor monads pattern vs. M. averyi. Similar results were obtained in the comparison of the observed and null proportions of ranges overlapping (Table II, Figure 6), in which the observed proportion for A. latifrons TD vs. M. diutius and M. lemniscatus were lower than most of the null proportions. For all other mimetic relationship, but A. latifrons TT vs. M. filiformis, the observed overlapping proportion was higher than null proportions. This result for A. latifrons TT vs. M. filiformis is expected since the observed proportion was rather low (only 11%, Table II; see also Figure 5).

DISCUSSION
Reports of mimicry involving snakes in the Amazonia usually suggest the relationship between a species of non-venomous coral snake imitating a single venomous coral snake species (1 vs. 1). Here we report on a mimetic complex involving a single polychromatic species of non-venomous coral snake (Almeida et al. 2014) and several sympatric species of Micrurus, with associations that include different color patterns and different levels of similarity, representing perfect and imperfect instances of mimicry. In Brazil, the only comparable instance of potential mimicry already tested in the literature involves the false coral snake Erythrolamprus aesculapii, for which populations from the Atlantic Forest appear to represent mimics of the monad patterned species M. corallinus, while populations from the Cerrado biome resemble the triad patterned species Micrurus frontalis and M. lemniscatus (Marques & Puorto 1991).
The mimicry assumptions in this study rely strictly on similarities of color patterns and distributional congruence based on objective data of localities of occurrence, without any additional analytical support. Experimental studies by Smith (1975) and Brodie III & Moore (1995) revealed an innate avoidance by birds of all ring color patterns, regardless of width and color. Nevertheless, we based our mimetic Table I. Beta coefficients (β) from multiple regression analysis between habitat suitability of chromatic patterns of Atractus latifrons (response variable) and bioclimatic variables plus habitat suitability of each Micrurus species (predictor variables). TD -tricolor dyads; TT -tricolor tetrads; TM -tricolor monads, and BM -bicolor monads. Bio1 -annual mean temperature, Bio7 -annual temperature range, Bio13 -precipitation of the wettest month, Bio14precipitation of the driest month, Bio18 -precipitation of the warmest quarter. associations on a more recent study work by Akcali & Pfennig (2017) who suggest that the greater the geographic overlap between models and mimics, the greater their similarity. Several studies have tested the similarity of chromatic patterns between models and mimics with respect to levels of geographic congruence (Harper 2006, Harper & Pfennig 2007, 2008, Pfennig et al. 2007, Kikuchi & Pfennig 2010, Akcali & Pfennig 2017. All of them agree that the relationship between geographical variation and mimetic precision supports the assumption that the mimic's accuracy increases as the model becomes rare (Pfennig et al. 2001, Ruxton et al. 2004, Harper & Pfennig 2007, Iserbyt et al. 2011. Our results regarding a positive mimetic relationship between A. latifrons as a perfect mimic of the apparently rare model M. albicinctus corroborate this tendency. On the other hand, the abundance and wide distribution of imperfect chromatic patterns, such as tricolor dyads in A. latifrons, may favor mimics through co-occurrence with a greater number of models, such as M. lemniscatus, M. diutius and M. filiformis (see Edmunds 2000, Kikuchi & Pfennig 2013. Indeed, we observed that the distribution of the pattern tricolor dyads of A. latifrons largely overlaps with areas of greatest Micrurus species richness (see Figure S4).
Among the proposed mimetic relationships involving Atractus latifrons and sympatric species of Micrurus in Amazonia (Cunha & Nascimento 1983, Savage & Slowinski 1992, Martins & Oliveira 1993, Silva 2004, Almeida et al. 2014), we were unable to assess the presumed mimetic relationship with M. hemprichii ortoni suggested by Silva (2004) due to the absence of the tricolor triads pattern in our sample of A. latifrons. Our distributional data corroborate previous evidence of perfect mimicry among some populations of A. latifrons and M. albicinctus by ecological niche modeling, the co-occurrence test and the strong similarity of color. Some populations of A. latifrons may also represent imperfect mimics, considering the relationships established among tricolor monads, tricolor dyads and tricolor tetrads of mimics and their models M. averyi, M. langsdorffi, M. lemniscatus, M. diutius, and M. filiformis, respectively. The combination of ENM methods promotes additional support to the assumption of nonrandom co-occurrence of color patterns for at least one mimetic system involving species with bicolor monads. Additional effects combining niche modeling techniques and distributional patterns of model and mimics of other systems may improve the accuracy of mimicry studies based on indirect data such as sympatry of color patterns and behavioral repertoire.
Finally, we note that for some Micrurus species and particular chromatic patterns of A. latifrons, the niche models that estimated potential geographic distributions relied on a limited number of occurrence records (e.g. M. albicinctus, M. averyi, bicolor patterns in monads and tetrads of Atractus latifrons), despite our exhaustive efforts at gleaning all information available in museum collections. Several studies point out that model performance generally decreases with sample size (Stockwell & Peterson 2002, Wisz et al. 2008, however this does not necessarily mean that the predicted distribution is inaccurate (e.g. Kadmon et al. 2003). For instance, a few records may suffice to characterize distributions of species with narrow environmental tolerances (Kadmon et al. 2003), certainly the case for the species analyzed here. Also, although niche modeling is expected to overpredict geographic distributions, our approach indicated it can be useful to investigate the patterns of species cooccurrence. Nevertheless, our results can only be considered as indications of the existence of a mimicry complex between these species, since it was not possible to indicate precisely which species acted as model or mimic. Further studies focused on these relationships would be required to confirm the mimicry complexes. We hope that new records from future sampling will provide additional support for our results, especially regarding the chromatic patterns of A. latifrons.

SUPPLEMENTARY MATERIAL
Appendix S1. Examined specimens of the Atractus latifrons and Micrurus.Countries are identified with capital letters and bold, states are only capitalized, counties in italics, and locations in simple text.
Appendix S2. Number of spatially unique occurrence records used in the niche modelling for each Micrurus species and Atractus latifrons chromatic pattern. Table SI. Spatially unique occurrence records used in the niche modelling for each Micrurus species and Atractus latifrons chromatic pattern. Table SII. Accuracy of ENM models for the ensemble suitability among modeling methods for each climatic model, and overall mean accuracy for the ensemble among modeling methods and climatic models. Figure S1. Modeled geographic distributions and occurrence records for Atractus latifrons chromatic patterns. Figure S2. Polygons generated by a convex hull for each combination of mimic A. latifrons and Micrurus model species. Figure S3. Modeled geographic distribution of Atractus latifrons, based on all occurrence records for this species. Figure S4. Richness pattern of Micrurus mimetic models in the Amazonia, resulting from stacking the geographic ranges from ENM.