Abstract

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 identified 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 benefit from the resemblance with other Micrurus species as perfect and imperfect mimics.

Key words
Co-occurrence; ecological niche modeling; models; mimics; color variation

INTRODUCTION

The selection of adaptive similarities involving mimics and models represents a classic problem in evolutionary biology (Kikuchi & Pfennig 2012KIKUCHI DW & PFENNIG DW. 2012. A Batesian mimic and its model share color production mechanisms. Curr Zool 58: 657-666.). The term Batesian mimicry refers to a phenomenon in which a non-noxious species replicates a visual, acoustic or chemical signal emitted by a venomous species (model), resulting in adaptive advantages against predation (Ruxton et al. 2004RUXTON GD, SHERRATT TN & SPEED MP. 2004. Avoiding Attack: the evolutionary ecology of crypsis, warning signals, and mimicry. Oxford University Press, 249., Hossie & Sherratt 2013HOSSIE TJ, SHERRATT TN. 2013. Defensive posture and eyespots deter avian predators from attacking caterpillar models. Anim Behav, 86: 383-389., Rabosky et al. 2016RABOSKY ARD, COX CL, RABOSKY DL, TITLE PO, HOLMES IA, FELDMAN A & MCGUIRE JA. 2016. Coral snakes predict the evolution of mimicry across New World snakes. Nat Comm 7: 11484.). In such systems – Batesian mimicry complexes – Ruxton et al. (2004)RUXTON GD, SHERRATT TN & SPEED MP. 2004. Avoiding Attack: the evolutionary ecology of crypsis, warning signals, and mimicry. Oxford University Press, 249. claim that imitators should occur only in sympatry with their models. However, Pfennig & Mullen (2010)PFENNIG DW & MULLEN SP. 2010. Mimics without models: causes and consequences of allopatry in Batesian mimicry complexes. Proc R Soc B Biol 277: 2577-2585. 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 co-occurrence of mimics and models may affect the abundance of perfect and imperfect mimics (Huheey 1976HUHEEY JE. 1976. Studies of warning coloration and mimicry. VII. Evolutionary Consequences of a Batesian-Müllerian Spectrum: a model for Müllerian Mimicry. Evolution 30: 86-93., Brodie III & Janzen 1995BRODIE EDIII & JANZEN FJ. 1995. Experimental studies of coralsnakes mimicry: Generalized avoidance of ringed snake patterns by free-ranging avian predators. Funct Ecol 9: 186-190., Pfennig et al. 2001PFENNIG DW, HARCOMBE WR & PFENNIG KS. 2001. Frequency-dependent Batesian mimicry. Nature 410(6826): 323., Caley & Schluter 2003CALEY MJ & SCHLUTER D. 2003. Predators favour mimicry in a tropical reef fish. Proc R Soc B Biol Sci 270: 667-672., Harper & Pfennig 2007HARPER GR JR & PFENNIG DW. 2007. Mimicry on the edge: why the mimics vary in resemblance to their model in different parts of their geographical range? Proc R Soc B Biol Sci, 274: 1955-1961., Kikuchi & Pfennig 2013KIKUCHI DW, PFENNIG DW. 2013. Imperfect mimicry and the limits of the natural selection. Univ Chic Press J 88(4): 297-315.). In areas where models are scarce or absent, the chances of predation on mimics would tend to increase (Huheey 1976HUHEEY JE. 1976. Studies of warning coloration and mimicry. VII. Evolutionary Consequences of a Batesian-Müllerian Spectrum: a model for Müllerian Mimicry. Evolution 30: 86-93., Pfennig et al. 2001PFENNIG DW, HARCOMBE WR & PFENNIG KS. 2001. Frequency-dependent Batesian mimicry. Nature 410(6826): 323., 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 2003CALEY MJ & SCHLUTER D. 2003. Predators favour mimicry in a tropical reef fish. Proc R Soc B Biol Sci 270: 667-672., Kikuchi & Pfennig 2013KIKUCHI DW, PFENNIG DW. 2013. Imperfect mimicry and the limits of the natural selection. Univ Chic Press J 88(4): 297-315.).

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 1996ROZE JA. 1996. Coral snakes of America: Biology, Identification, and Venoms. Krieger Publishing Company, p. 1-328.). 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 1981GREENE HW & MCDIARMID RW. 1981. Coral Snake Mimicry: does it occur? Science 213, 4513: 1207-1212., Campbell & Lamar 1989CAMPBELL JA & LAMAR WW. 1989. The Venomous Reptiles of Latin America. Comst Pub Assoc. Cornell University Press. Ithaca. New York., 2004CAMPBELL JA & LAMAR WW. 2004. The Venomous Reptiles of the Western Hemisphere. Volume I. Comstock Publishing Associates. Cornell University Press. Ithaca. New York., Savage & Slowinski 1992SAVAGE JM & SLOWINSKI JB. 1992. The colouration of the venomous coral snakes (Family Elapidae) and their mimics (families Aniliidae and Colubridae). Biol J Linn Soc 45: 235-254., Roze 1996ROZE JA. 1996. Coral snakes of America: Biology, Identification, and Venoms. Krieger Publishing Company, p. 1-328., Almeida et al. 2016ALMEIDA PCR, PRUDENTE ALC, CURCIO FF & RODRIGUES MTU. 2016. Biologia e História Natural das Cobras-Corais. In: Silva NJ Jr (Ed), As Cobras Corais do Brasil: Biologia, Taxonomia, Venenos e Envenenamentos. PUC Goiás, Goiânia, p. 168-215.).

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 1981GREENE HW & MCDIARMID RW. 1981. Coral Snake Mimicry: does it occur? Science 213, 4513: 1207-1212., Marques & Puorto 1991MARQUES OAV & PUORTO G. 1991. Padrões cromáticos, distribuição e possível mimetismo em Erythrolamprus aesculapii (Serpentes: Colubridae). Mem Inst But 53(1): 127-134., Brodie III & Janzen 1995BRODIE EDIII & MOORE A. 1995 Experimental studies of coralsnakes mimicry: do snakes mimic millipedes? Anim Behav 49: 534-536., Torre-Loranza et al. 2006TORRE-LORANZA MA, LÉON GA & LUNA MAL. 2006. Coralillos verdaderos (Serpentes: Elapidae) y coralillos falsos (Serpentes: Colubridae) de Veracruz, México. Acta Zoológica Mex (nueva serie) 22: 3: 1-13., Almeida et al. 2014ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80., Bosque et al. 2016BOSQUE RJ, NOONAN BP & COLLI G R. 2016. Geographical coincidence and mimicry between harmless snakes (Colubridae: Oxyrhopus) and harmful models (Elapidae: Micrurus). Glob Ecol Biogeogr 25: 218-226.). Combining phylogenetic techniques, distributional comparisons, and general similarity, Rabosky et al. (2016)RABOSKY ARD, COX CL, RABOSKY DL, TITLE PO, HOLMES IA, FELDMAN A & MCGUIRE JA. 2016. Coral snakes predict the evolution of mimicry across New World snakes. Nat Comm 7: 11484. demonstrated a high correlation between the occurrence of red and black band patterns (RBB, sensu Rabosky et al. 2016RABOSKY ARD, COX CL, RABOSKY DL, TITLE PO, HOLMES IA, FELDMAN A & MCGUIRE JA. 2016. Coral snakes predict the evolution of mimicry across New World snakes. Nat Comm 7: 11484.) 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)BRATTSTROM BH. 1955. The coral snake “mimic” problem and protective coloration. Evolution 9: 217-219., 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.

Recently, Akcali & Pfennig (2017)AKCALI CK & PFENNIG DW. 2017. Geographic variation in mimetic precision among different species of coral snake mimics. J Evol Biol 30: 1420-1428. showed that resemblance level of mimics depends on their degree of sympatry (deep, edge, and allopatry) with the model, since each species exhibits different responses to selective pressures. For instance, the polychromatism of the false coral snake Atractus latifrons and its geographic congruence with species of Micrurus in Amazonia (M. albicinctus Amaral 1925AMARAL A. 1925. South American Snakes in the collection of the United States National Museum. Proc Unit Stat Nat Mus. 67, 24: 1-30., M. averyi Schmidt 1939SCHMIDT KP. 1939. A new coral snake from British Guyana. Field Mus Nat Hist Publ Zool Ser 24: 45-47., M. filiformis (Günther 1859GÜNTHER A. 1859. On the genus Elaps of Wagner. Proc Zool Soc London 1859: 79-89.), M. langsdorffi Wagler 1824WAGLER J. 1824. Serpentum Brasiliensium species novae, ou historie naturelle des espèces nouvelles de serpens. In: Jean de Spix, Animalia nova sive species novae. [NAtrix bahiensis: 27. Monaco Typis Franc. Seraph. Hübschmanni, vii + 75 p. , M. lemniscatus (Linnaeus 1758LINNAEUS C. 1758. Systema naturae per regna tria naturae, secundum classes, ordines, genera species, cum characteribus, differentiis, synonymis, locis. Tomus I. Editio decima, reformata Laurentii Salvii Holmae, 10th ed., 824 p.) and M. diutius Burger 1955BURGER WL. 1955. A new subspecies of the coral snake Micrurus lemniscatus from Vénézuela, British Guiana and Trinidad, and a key for the identification of associated species of coral snakes. Bol Mus Cienc Nat Caracas 1: 35-40.) indicated a high distributional correspondence between mimics and models (Almeida et al. 2014ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80.). Furthermore some authors have already been suggested a possible mimetic relationships between A. latifrons and Micrurus in the Amazon (Cunha & Nascimento 1983CUNHA OR & NASCIMENTO FP. 1983. Ofídios da Amazônia XX - As espécies de Atractus Wagler, 1828, na Amazônia Oriental e Maranhão (Ophidia, Colubridae). Bol Mus Paraense Emílio Goeldi, Nova Série Zoologia 123: 1-40., Savage & Slowinski 1992SAVAGE JM & SLOWINSKI JB. 1992. The colouration of the venomous coral snakes (Family Elapidae) and their mimics (families Aniliidae and Colubridae). Biol J Linn Soc 45: 235-254., Martins & Oliveira 1993MARTINS M & OLIVEIRA ME. 1993. The snakes of the genus Atractus Wagler (Reptilia: Squamata: Colubridae) from the Manaus region, central Amazonia, Brazil Zool Med 67: 1-26, 21-40., Silva 2004SILVA JJH. 2004. Las Serpientes del género Atractus Wagler, 1828 (Colubridae, Xenodontinae) en la Amazonia Colombiana. Rev Acad Colomb Cienc XXVIII 108: 409-446., Almeida et al. 2014ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80.).

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 2012RANGEL TF & LOYOLA RD. 2012. Labeling Ecological Niche Models. Nat Cons 10(2): 119-126.). From such associations, one can identify and map new locations where viable populations of species are likely to occur (Peterson et al. 2011PETERSON AT, SOBERÓN J, PEARSON RG, ANDERSON RP, MARTINÉZ-MEYER E, NAKAMURA M & ARAÚJO MB. 2011. Ecological Niches and Geographic Distributions. Princeton University Press, New Jersey.). 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. 2011PETERSON AT, SOBERÓN J, PEARSON RG, ANDERSON RP, MARTINÉZ-MEYER E, NAKAMURA M & ARAÚJO MB. 2011. Ecological Niches and Geographic Distributions. Princeton University Press, New Jersey., Uribe-Rivera et al. 2017URIBE-RIBERA DE, SOTO-AZAT C, VALENZUELA-SANCHEZ A, BIZAMA G, SIMONETTI JA & PLISCOFF P. 2017. Dispersal and extrapolation on the accuracy of temporal predictions from distribution models for the Darwin’s frog. Ecol Soc Am 27: 1633-1645.). 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. 2008JIMÉNEZ-VALVERDE A, LOBO JM & HORTAL J. 2008. Not as good as they seem: the importance of concepts in species distribution modelling. Biodiv Let 14: 885-890., Peterson et al. 2011PETERSON AT, SOBERÓN J, PEARSON RG, ANDERSON RP, MARTINÉZ-MEYER E, NAKAMURA M & ARAÚJO MB. 2011. Ecological Niches and Geographic Distributions. Princeton University Press, New Jersey.), the distributional area estimated by ENM could provide a good approximation of the occupied distributional area at large scales (Peterson et al. 2011PETERSON AT, SOBERÓN J, PEARSON RG, ANDERSON RP, MARTINÉZ-MEYER E, NAKAMURA M & ARAÚJO MB. 2011. Ecological Niches and Geographic Distributions. Princeton University Press, New Jersey.). Therefore, ENM represents an important tool for overcoming Wallacean shortfalls (Lomolino 2004LOMOLINO MV. 2004. Frontiers of Biogeography. In: Lomolino M & Heaney L (Eds), Frontiers of Biogeography: New Directions in the Geography of Nature. Sunderland, MA: Sinauer., Hortal et al. 2015HORTAL J, BELLO F, DINIZ-FILHO JAF, LEWINSOHN TM, LOBO JM & LADLE RJ. 2015. Seven shortfalls that beset large-scale knowledge of biodiversity. Annu Rev Ecol Evol Syst. 46: 523-549.) 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.

MATERIALS AND METHODS

We examined 157 specimens of Atractus latifrons and 402 specimens of coral snakes of the genus Micrurus: M. albicinctus (n= 7), M. averyi (n= 19), M. filiformis (n= 74); M. langsdorffi (n= 24), M. lemniscatus (n= 325), and M. diutius (n= 109) with an Amazonian distribution, housed in 21 institutions (Supplementary Material – Appendix S1); institutional abbreviations followed Sabaj (2016)SABAJ MH. 2016. Standard symbolic codes for institutional resource collections in herpetology and ichthyology: an Online Reference. Version 6.5 (16 August 2016). Electronically accessible at http://www.asih.org/, American Society of Ichthyologists and Herpetologists, Washington, DC.
http://www.asih.org/, American Society o... , except for Centro de Estudos e Pesquisas Biológicas (CEPB) and Coleção Herpetológica da Universidade Federal do Ceará (CHUFC), not listed there). Nomenclature of the color patterns of A. latifrons followed Almeida et al. (2014)ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80.. 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)ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80. 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).

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 species-level identification of most of the Micrurus taxa followed Roze (1996)ROZE JA. 1996. Coral snakes of America: Biology, Identification, and Venoms. Krieger Publishing Company, p. 1-328., except the for M. lemniscatus species group, in which we followed Pires et al. (2014)PIRES MG, SILVA-JR NJ, FEITOSA DT, FILHO GAP, PRUDENTE ALC & ZAHER H. 2014. A new species of triadal coral snake of the genus Micrurus Wagler, 1824 (Serpentes: Elapidae) from northeastern Brazil. Zootaxa 3811(4): 569-584..

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)ARAÚJO MB & NEW M. 2007. Ensemble forecasting of species distributions. Trends Ecol Evol 22: 42-7. to generate consensus predictions by combining outputs from different niche modeling methods (weighted by individual model accuracy) and climatic models (Rangel & Loyola 2012RANGEL TF & LOYOLA RD. 2012. Labeling Ecological Niche Models. Nat Cons 10(2): 119-126.). We obtained a set of bioclimatic variables at 0.5^o resolution from the ecoClimate database (http://ecoclimate.org, Lima-Ribeiro et al. 2015LIMA-RIBEIRO MS, VARELA S, GONZÁLEZ-HERNÁNDEZ J, OLIVEIRA G, DINIZ-FILHO JAF & TERRIBILE LC. 2015. Ecoclimate: a database of climate data from multiple models for past, present, and future for macroecologists and biogeographers. Bio Inf 10: 1-21.), 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. 2012TERRIBILE LC, LIMA-RIBEIRO MS, ARAÚJO MB, BIZÃO N, COLLEVATTI RG & DOBROVOLSKI R. 2012. Areas of climate stability of species ranges in the Brazilian Cerrado: disentangling uncertainties through time. Nat Conserv 10: 152-159.). 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)FRANKLIN J. 2009. Mapping species distributions: spatial inference and prediction. Cambridge: Camb Univ Press. and Peterson et al. (2011)PETERSON AT, SOBERÓN J, PEARSON RG, ANDERSON RP, MARTINÉZ-MEYER E, NAKAMURA M & ARAÚJO MB. 2011. Ecological Niches and Geographic Distributions. Princeton University Press, New Jersey..

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 pseudo-absences 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. 2006ALLOUCHE O, TSOAR A & KADMON R. 2006. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J Appl Ecol 43: 1223-1232., Diniz-Filho et al. 2009DINIZ-FILHO JAF, BINI ML, RANGEL FT, LOYOLA RD, HOF C & NOGUÉS-BRAVO D. 2009. Partitioning and mapping uncertainties in ensembles of forecasts of species turnover under climate change. Ecography (Cop). 32: 897-906.). 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. 2009DINIZ-FILHO JAF, BINI ML, RANGEL FT, LOYOLA RD, HOF C & NOGUÉS-BRAVO D. 2009. Partitioning and mapping uncertainties in ensembles of forecasts of species turnover under climate change. Ecography (Cop). 32: 897-906., Terribile et al. 2012TERRIBILE LC, LIMA-RIBEIRO MS, ARAÚJO MB, BIZÃO N, COLLEVATTI RG & DOBROVOLSKI R. 2012. Areas of climate stability of species ranges in the Brazilian Cerrado: disentangling uncertainties through time. Nat Conserv 10: 152-159.).

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 2004GOTELLI NJ & ENTSMINGER GL. 2004. EcoSim: Null models software for ecology. Version 7. Acquired Intelligence Inc. & Kesey-Bear.). 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 1990STONE L & ROBERTS A. 1990. The checkboard score and species distribution. Oecologia 85: 74-79.) 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.

RESULTS

Morphological analyses

Based on the correspondence of color pattern among Atractus latifrons and Micrurus species, we recognized six species of Micrurus as potential mimetic models for A. latifrons: M. albicinctus, M. averyi, M. filiformis, M. langsdorffi, M. lemniscatus and M. diutius (Figure 1). Therefore, based on similarities of general color patterns, we assumed the following associations: the bicolor monads pattern of Atractus latifrons resembles Micrurus albicinctus (Figure 2a–b); the tricolor monads pattern 1 of A. latifrons resembles M. averyi (Figure 2c–d); the tricolor monads pattern 2 resembles M. langsdorffi (Figure 2e–f); the tricolor dyads pattern of A. latifrons resembles the tricolor triads pattern of M. lemniscatus and M. diutius species (Figure 3a–c) and the tricolor tetrads pattern of A. latifrons resembles the tricolor triads pattern of M. filiformis (Figure 4a–b), although with a few number of triads along the body M. diutius and M. lemniscatus can also serve a models for A. latifrons (Figure 3a-c).

Figure 1
Color patterns of possible mimetic models of genus Micrurus: a) Micrurus albicinctus; b) Micrurus averyi; c) Micrurus langsdorffi; d) Micrurus filiformis; e) Micrurus lemniscatus; and f) Micrurus diutius.

Figure 2
Chromatic bicolor monads patterns of mimic Atractus latifrons (a) and model Micrurus albicinctus (b); tricolor monads pattern morph 1 of A. latifrons (c) and model Micrurus averyi (d); and the tricolor monads pattern morph 2 of A. latifrons (e) and model Micrurus langsdorffi (f). Scale: 10mm.

Figure 3
Chromatic tricolor dyads patterns of mimic Atractus latifrons (a) and model Micrurus lemniscatus (b); and Micrurus diutius (c). Scale: 10mm.

Figure 4
Chromatic tricolor tetrads patterns of mimic Atractus latifrons (a) and model Micrurus filiformis (b) Scale: 10mm.

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).

Figure 5
Geographic distribution from ENM for chromatic patterns of Atractus latifrons and correspondent coralsnake mimetic model.

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) the co-occurrence of the bicolor monads pattern of A. latifrons and M. albicinctus in a restricted region of Amazonia, specifically the states of Rondônia and Amazonas; (2) the tricolor dyads pattern of A. latifrons co-occurs with M. lemniscatus and M. diutius as potential models in Amazonia, specifically in the states of Pará and Amazonas; (3) the tricolor tetrads pattern of A. latifrons co-occurs with M. filiformis and M. lemniscatus in central Amazonia, and (4) the tricolor monads pattern of A. latifrons co-occurs with the potential models M. langsdorffi 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 co-occurrence test from the polygons obtained by convex hulls, except for the pattern A. latifrons TM vs M. averyi, which was non-significant (i.e., 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).

Figure 6
Histograms of simulated proportions of range overlapping derived from the 1000 resampled geographic distributions for each chromatic pattern. The red line indicates the observed proportion of range overlapping between each A. latifrons chromatic pattern and the respective Micrurus species.

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. 2014ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80.) 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 1991MARQUES OAV & PUORTO G. 1991. Padrões cromáticos, distribuição e possível mimetismo em Erythrolamprus aesculapii (Serpentes: Colubridae). Mem Inst But 53(1): 127-134.).

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 associations on a more recent study work by Akcali & Pfennig (2017)AKCALI CK & PFENNIG DW. 2017. Geographic variation in mimetic precision among different species of coral snake mimics. J Evol Biol 30: 1420-1428. 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 2006HARPER GR. 2006. Evolution of a Snake Mimicry Complex. PhD thesis, Univ North Carolina., Harper & Pfennig 2007, 2008HARPER GR & PFENNIG DW. 2008. Selection overrides gene flow to break down maladaptive mimicry. Nature 451:1103-1106., Pfennig et al. 2007PFENNIG DW, HARPER GR, BRUMO AF, HARCOMBE WR & PFENNIG KS. 2007. Population differences in predation on Batesian mimics in allopatry with their model: Selection against mimics is strongest when they are common. Behav Ecol Sociobiol 61: 505-511., Kikuchi & Pfennig 2010KIKUCHI D & PFENNIG DW. 2010. Predator Cognition Permits Imperfect Coral Snake Mimicry. Am Nat 176(6): 830-834., Akcali & Pfennig 2017AKCALI CK & PFENNIG DW. 2017. Geographic variation in mimetic precision among different species of coral snake mimics. J Evol Biol 30: 1420-1428.). 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. 2001PFENNIG DW, HARCOMBE WR & PFENNIG KS. 2001. Frequency-dependent Batesian mimicry. Nature 410(6826): 323., Ruxton et al. 2004RUXTON GD, SHERRATT TN & SPEED MP. 2004. Avoiding Attack: the evolutionary ecology of crypsis, warning signals, and mimicry. Oxford University Press, 249., Harper & Pfennig 2007, Iserbyt et al. 2011ISERBYT A, BOTS J, VAN DONGEN S, TING JJ, VAN GOSSUM H & SHERRATT TN. 2011. Frequency-dependent variation in mimetic fidelity in an intraspecific mimicry system. Proc Biol Sci 278: 3116-3122.). 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 2000EDMUNDS M. 2000. Why are there good and poor mimics? Biol J Linn Soc 70: 459-466., Kikuchi & Pfennig 2013KIKUCHI DW, PFENNIG DW. 2013. Imperfect mimicry and the limits of the natural selection. Univ Chic Press J 88(4): 297-315.). 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 1983CUNHA OR & NASCIMENTO FP. 1983. Ofídios da Amazônia XX - As espécies de Atractus Wagler, 1828, na Amazônia Oriental e Maranhão (Ophidia, Colubridae). Bol Mus Paraense Emílio Goeldi, Nova Série Zoologia 123: 1-40., Savage & Slowinski 1992SAVAGE JM & SLOWINSKI JB. 1992. The colouration of the venomous coral snakes (Family Elapidae) and their mimics (families Aniliidae and Colubridae). Biol J Linn Soc 45: 235-254., Martins & Oliveira 1993MARTINS M & OLIVEIRA ME. 1993. The snakes of the genus Atractus Wagler (Reptilia: Squamata: Colubridae) from the Manaus region, central Amazonia, Brazil Zool Med 67: 1-26, 21-40., Silva 2004SILVA JJH. 2004. Las Serpientes del género Atractus Wagler, 1828 (Colubridae, Xenodontinae) en la Amazonia Colombiana. Rev Acad Colomb Cienc XXVIII 108: 409-446., Almeida et al. 2014ALMEIDA PC, FEITOSA DT, PASSOS P & PRUDENTE ALC. 2014. Morphological variation and taxonomy of Atractus latifrons (Günther, 1868) (Serpentes: Dipsadidae). Zootaxa 3860: 64-80.), we were unable to assess the presumed mimetic relationship with M. hemprichii ortoni suggested by Silva (2004)SILVA JJH. 2004. Las Serpientes del género Atractus Wagler, 1828 (Colubridae, Xenodontinae) en la Amazonia Colombiana. Rev Acad Colomb Cienc XXVIII 108: 409-446. 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 non-random 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 2002STOCKWELL DRB & PETERSON AT. 2002. Effects of sample size on accuracy of species distribution models. Ecol Model 148: 1-13., Wisz et al. 2008WISZ MS, HIJMANS RJ, LI J, PETERSON AT, GRAHAM CH & GUISAN A & ‘NCEAS, PREDICTING DISTRIBUTIONS WORKING GROUP’. 2008. Effects of sample size on the performance of species distribution models. Biodiver Res 14: 763-773.), however this does not necessarily mean that the predicted distribution is inaccurate (e.g. Kadmon et al. 2003KADMON R, FARBER O & DANIN A. 2003. A systematic analysis of factors affecting the performance of climatic envelope models. Ecol Appl 13: 853-867.). For instance, a few records may suffice to characterize distributions of species with narrow environmental tolerances (Kadmon et al. 2003KADMON R, FARBER O & DANIN A. 2003. A systematic analysis of factors affecting the performance of climatic envelope models. Ecol Appl 13: 853-867.), 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 co-occurrence. 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.

ACKNOWLEDGMENTS

We thank the curators of the following collections for their assistance in loaning study material: Zoologische Staatssammlung München (ZSM); Zoologische Museum Hamburg (ZMH); Senckenberg Museum Frankfurt (SMF); Zoologisches Museum zu Berlin (ZMB); Naturhistorisches Museum Wien (NHMW); Museo Argentino de Ciencias Naturales Bernardino Rivadavia (MACN); Museo de Ciencias Naturales de La Plata (MCNLP); Museo Noel Kempff Mercado (MNKR); Centro de Estudos e Pesquisas Biológicas da Pontifícia Universidade Católica de Goiás (CEPB); Fundação Medicina Tropical (FMT); Instituto Nacional de Pesquisas da Amazônia (INPA); Instituto Butantan de São Paulo (IBSP); Museu de Zoologia da Universidade de São Paulo (MZUSP); Museu Paraense Emílio Goeldi (MPEG); Universidade Federal do Mato Grosso (UFMT); Museu Nacional do Rio de Janeiro (MNRJ); Instituto Alexander von Humboldt (IAvH); Museo de la Universidad La Salle (MLS); Universidad Nacional de Colombia (ICN); Museo de Historia Natural de la Universidad del Cauca (MHNUC); British Museum of Natural History (BMNH); Museo de Historia Natural de la Universidad Mayor de San Marcos (MHNSM); Naturhistoriska Riksmuseet (NRM). We are grateful to the Museu Paraense Emílio Goeldi (MPEG) for logistical support. Financial support for ALCP was provided by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) productivity scholarship grant number 302611/2018-5 and PROTAX 44.0413/2015-0; PA – PROTAX grant 562171/2010-0); for PCRA was provided by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq; processes PROTAX 562171/2010-0). The research by LCT is supported CNPq productivity grants and by National Institutes for Science and Technology (INCT) in Ecology, Evolution and Biodiversity Conservation, supported by MCTIC/ CNPq (proc. 465610/2014-5) and FAPEG. We thank Guilherme Oliveira for helping us with niche modeling. We also thank Glenn Shepard (native speaker) for English revisions, and three anonymous reviewers for their helpful comments and suggestions.

SUPPLEMENTARY MATERIAL

REFERENCES

Publication Dates

History