Multiple Routes to Color Convergence in a Radiation of Neotropical Poison Frogs

1247 Syst. Biol. 72(6):1247–1261, 2023 https://doi.org/10.1093/sysbio/syad051 Advance Access Publication August 10, 2023 Spotlight Article © The Author(s) 2023. Published by Oxford University Press on behalf of the Society of Systematic Biologists. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Multiple Routes to Color Convergence in a Radiation of Neotropical Poison Frogs EvanTwomey1,*, PauloMelo-Sampaio2, Lisa M.Schulte1, FrankyBossuyt3, Jason L.Brown4, and SantiagoCastroviejo-Fisher5 1Department of Wildlife/Zoo Animal Biology and Systematics, Faculty of Biological Sciences, Goethe University Frankfurt, Max-von-Laue- Str. 13, Frankfurt am Main 60438, Germany, 2Departamento de Vertebrados, Museu Nacional, Universidade Federal do Rio de Janeiro, R. Gen. Herculano Gomes 41, Rio de Janeiro 20941- 360, Brazil, 3Amphibian Evolution Laboratory, Biology Department, Vrije Universiteit Brussel, Pleinlaan 2, Brussels 1050, Belgium, 4School of Biological Sciences, Southern Illinois University, 125 Lincoln Dr., Carbondale, IL 62901, USA, and 5Departamento de Zoología, Universidad de Sevilla, Av. de la Reina Mercedes, Seville 41012, Spain *Correspondence to be sent to: Department of Wildlife/Zoo Animal Biology and Systematics, Faculty of Biological Sciences, Goethe University Frankfurt, Frankfurt am Main, Germany; E-mail: [email protected]. Received 22 June 2022; reviews returned 3 August 2023; accepted 9 August 2023 Associate Editor: RaynaBell Abstract.—Convergent evolution is defined as the independent evolution of similar phenotypes in different lineages. Its existence underscores the importance of external selection pressures in evolutionary history, revealing how functionally similar adaptations can evolve in response to persistent ecological challenges through a diversity of evolutionary routes. However, many examples of convergence, particularly among closely related species, involve parallel changes in the same genes or developmental pathways, raising the possibility that homology at deeper mechanistic levels is an important facilitator of phenotypic convergence. Using the genus Ranitomeya, a young, color-diverse radiation of Neotropical poison frogs, we set out to 1) provide a phylogenetic framework for this group, 2) leverage this framework to determine if color phenotypes are convergent, and 3) to characterize the underlying coloration mechanisms to test whether color convergence occurred through the same or different physical mechanisms. We generated a phylogeny for Ranitomeya using ultraconserved elements and investigated the physical mechanisms underlying bright coloration, focusing on skin pigments. Using phylogenetic comparative methods, we identified several instances of color convergence, involving several gains and losses of carotenoid and pterin pigments. We also found a compelling example of nonparallel convergence, where, in one lineage, red coloration evolved through the red pterin pigment drosopterin, and in another lineage through red ketocarotenoids. Additionally, in another lineage, “reddish” coloration evolved predominantly through structural color mechanisms. Our study demonstrates that, even within a radiation of closely related species, convergent evolution can occur through both parallel and nonparallel mechanisms, challenging the assumption that similar phenotypes among close relatives evolve through the same mechanisms. [Carotenoids; coloration; convergent evolution; Dendrobatidae; pterins; ultraconserved elements.] The independent evolution of similar phenotypes in different lineages—convergent evolution—is a striking and ubiquitous feature of life (McGhee 2011). Well-known examples of convergent evolution include image-forming eyes in dinoflagellates, cephalopods, and vertebrates (Leander 2008), repeated pelvic reduction in sticklebacks (Bell et al. 1993), and recurrent ecomorph evolution in Anolis lizards (Losos et al. 1998). Because convergent evolution indicates that certain phenotypes evolve repeatedly and predictably, its existence supports the view that “replaying life’s tape” would result in similar, deterministic outcomes due to external forces such as natural selection. The best examples supporting this view are those of complex and clearly adaptive traits evolving across deep timescales (Powell and Mariscal 2015), which minimizes the extent to which homologous mechanisms (e.g., conserved genes, pathways) can explain the convergence (Leander 2008). On the other hand, because closely related species share more homologous traits than distantly related species, one expects (and frequently observes) that convergence among close relatives is due to homologous underlying mechanisms (frequently termed “parallel evolution”) (Shapiro et al. 2006; Rosenblum et al. 2010; Conte et al. 2012). Such cases support the view that internal forces (e.g., homology) are important facilitators of recurrent phenotypic evolution (Pearce 2012). However, there are a handful of examples of “nonparallel” convergence (i.e., convergence via different mechanisms, sensu Rosenblum et al. 2014) among close relatives, as well as parallel convergence among distant relatives (e.g., shared bone structure in tetrapod wings) demonstrating that convergence can occur through parallel or nonparallel pathways regardless of phylogenetic depth (Calboli et al. 2003; Wilkens and Strecker 2003; Leander 2008; Steiner et al. 2009; Kemppainen et al. 2021). As with distant relatives, nonparallel convergence among close relatives better demonstrates the importance of external selection pressures in that homologous mechanisms are not a prerequisite for convergent phenotypes to evolve. Thus, the extent to which convergence can inform us about evolution is in large part dependent on our understanding of the mechanistic basis of convergent phenotypes and determining whether convergence Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 SYSTEMATIC BIOLOGY 1248 among close relatives occurs through similar or different mechanisms. It is important to clarify two points when discussing parallel versus nonparallel evolution. The first is that phenotypes have hierarchical levels of causality (e.g., nucleotide, gene, pathway, physical mechanism), so convergence can occur through similarity at different causal levels. Independent changes to the same gene (or even the same nucleotide) can yield identical phenotypic effects; similarly, different genes can be co-opted for the same function (Elmer and Meyer 2011), as can different physical mechanisms (e.g., pigmentary vs. structural coloration, discussed below). Second, although phenotypic convergence is homoplasy, some of the causative mechanisms that underlie these phenotypes may be homologous. For example, mandibular (“true”) teeth were lost in the ancestor to frogs over 200 million years ago, and regained in a single frog species, Gastrotheca guentheri (Boulenger 1882; Wiens 2011). Given their anatomical similarity to other amphibian teeth, this suggests that true teeth in frogs arose from an ancient developmental pathway shared among all amphibians, perhaps through the deactivation of a suppressor regulatory gene (Paluh et al. 2021). In this example, it is easy to see how a conserved pathway for tooth development may have facilitated the re-evolution of this homoplasic trait. Animal coloration has played a prominent role in the study of convergence (Bates 1862; Hoekstra et al. 2006; Rosenblum et al. 2010; Friedman et al. 2014; Kratochwil et al. 2018). Bright, (i.e., non-melanic) coloration is produced from a wide variety of physical mechanisms, meaning that ample opportunity exists for convergence via different evolutionary paths. For example, blue coloration in animals is usually structural (but see Goda and Fujii 1995; Rüdiger et al. 1968; Taboada et al. 2020; Yu et al. 2008 for examples of blue pigments). However, a wide variety of blue-producing structures have evolved, including guanine-based iridophores (fish, reptiles, and amphibians; Bagnara et al. 2007), reflectin-based iridophores (cephalopods; Crookes et al. 2004), collagen arrays (bird and mammal skin; Prum and Torres 2003, 2004), the spongy medullary layer of avian feather barb rami (Prum 2006), and beta-keratin/ melanin/air matrices in avian feather barbules (Prum 2006). Other colors, such as yellow, orange, and red, are more commonly produced from pigments such as carotenoids and pterins. However, structural color can also play an important or even predominant role in generating these colors (Morrison 1995; Teyssier et al. 2015). Overall, convergent color evolution is not only common, but can occur through parallel and nonparallel mechanisms. Here, we study the physical mechanisms underlying a color radiation in the Neotropical poison frog genus Ranitomeya (family Dendrobatidae) in a comparative phylogenetic framework. This is a relatively young group (~10 million years; Santos et al. 2009) consisting of 16 species, and exhibits extreme color variation, spanning the visible spectrum from blue to red (Brown et al. 2011). Dendrobatids, like many other amphibians, produce color through a stacked arrangement of dermal cells often referred to as the “dermal chromatophore unit” (Bagnara et al. 1968; Kuriyama et al. 2006; Twomey et al. 2020b). The most superficial layer of the dermal chromatophore unit are the xanthophores/erythrophores, which contain brightly colored pigments such as carotenoids and pterins (Obika and Bagnara 1964; Bagnara 1966). These two classes of pigments originate from unrelated pathways, being of dietary and endogenous origin, respectively (Giordano et al. 2003; Ziegler 2003; Toews et al. 2017). Below are iridophores, which contain stacks of refractive guanine crystals interspersed with cytoplasm that produce structural coloration (Bagnara et al. 2007). Melanophores are the deepest cell type, containing the black pigment melanin. Besides producing black coloration, melanin contributes to structural color by absorbing incoherently scattered light (Shawkey and Hill 2006). Importantly, due to the broad overlap in terms of colors produced by each mechanism, this means that similar colors in different species can be produced by various combinations of carotenoids, pterins, and iridophores, or by each mechanism alone (Teyssier et al. 2015; Twomey et al. 2020b; Stuart-Fox et al. 2021a). Specifically, we ask the following question: Does color convergence among closely related species occur through the same or different mechanisms? To address this, we provided a phylogenetic framework for this group, including all 16 nominal species, using ultraconserved elements. Within this framework, we analyzed the evolution of bright dorsal coloration to examine whether similar phenotypes are convergent. Finally, we analyzed skin pigments and spectral reflectance data to characterize the underlying mechanisms of color convergence, allowing us to identify common or independent evolutionary causes of the same observed phenotypes. Materials and Methods Sample Collection and Data Use We aimed to sample the color diversity in Ranitomeya as completely as possible across nominal species and color morphs of polymorphic species. To do so, we sampled Ranitomeya populations throughout Peru and western Brazil. During our sampling in western Brazil, we also discovered several Ranitomeya populations with unique coloration but of unclear taxonomic status due to unresolved species limits. For this reason, the terminals in our phylogeny have no consistent taxonomic rank and are herein referred to as taxa. We collected pigment and reflectance data from 11 Ranitomeya taxa throughout the Amazonian regions of eastern Peru and western Brazil. For an additional 10 taxa, we took pigment and reflectance data from previous studies on color mechanisms in poison frogs (Twomey et al. 2020a, 2020b). Supplementary Table S1 provides information on sampling localities and sample sizes for the pigment and reflectance data sets. For VOL. 72 Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 TWOMEY ET AL. - COLOR CONVERGENCE IN POISON FROGS2023 1249 the metrics used in the phylogenetic comparative analyses (dorsal color metric PC1, reflectance λmax, absorbance-based pigment concentrations), we calculated each metric on a per-individual basis and then averaged at the level of taxon. Of the 38 samples included in the phylogenetic analysis, 15 were from Guillory et al. (2019), while the remaining 23 are new to this study. See Supplementary Table S2 for overview of sampling localities and accession numbers for samples in the phylogenetic analysis. Reflectance Spectrometry and Quantification of Dorsal Color Reflectance spectrometry was done with a Stellarnet Silver-Nova spectrometer (190–1100 nm range) and Spectrawiz software version 5.33. Reflectance measurements were taken on each frog at 4 points on the brightly colored parts of the dorsum (right/left head and right/ left mid-back, excluding black regions) and averaged to a single reflectance spectrum, except for frogs whose head and body were a different color (i.e., R. fantastica and R. benedicta, Fig. 1), in which case only the two head measurements were averaged. Reflectance measurements were taken relative to a halon white reflectance standard, with illumination and read angle at 45° relative to the skin, at a fixed distance of 3 mm using a custom 3D-printed probe tip. UV-VIS and VIS-NIR measurements were taken separately, using a Stellarnet SL3 deuterium lamp for the former and a Stellarnet SL1 halogen lamp for the latter. As the two readings overlapped from 400 to 700 nm, it was possible to splice the two spectra to a single spectrum with a clean signal from approximately 220–1100 nm. Reflectance spectra were plotted and analyzed in R version 3.6.3 (R Core Team 2020) using the Pavo2 package (Maia et al. 2019). As the wavelength of peak reflectance (λmax, commonly referred to as hue or spectral location) has been shown to be correlated to the thickness of iridophore platelets in poison frogs and other organisms (Morrison 1995; Twomey et al. 2020b), we used λmax of the reflectance spectrum as a proxy for structural color hue and included this trait in our comparative analysis. To investigate dorsal color evolution in Ranitomeya, we used two methods. First, for visualization purposes, we measured RGB values from dorsal photographs and mapped the individual color channels onto the phylogeny as shown on Figure 1 (additional details are given in the Supplementary Material). For a more rigorous analysis of dorsal color, we also derived a dorsal color metric from the spectral reflectance data. For this analysis, the main challenge was to reduce the dimensionality of the spectral reflectance data to one or two meaningful color metrics that could be analyzed. To do this, we used the vismodel and colspace functions of Pavo2 to calculate quantum catch data for each of 4 cone types in a tetrachromat visual system, reducing the dimensionality to 4 columns per individual (i.e., quantum catch values for u, s, m, and l cones). For each individual, these calculations were done separately on each dorsal reflectance measurement, and quantum catches were averaged per cone type per individual. Vismodel was calculated with the following parameters: visual = “avg.uv,” achromatic = “bt.dc,” relative = TRUE, illum = “ideal). This visual system was chosen as it represents an average avian UV visual system, and is thus representative of an ecologically important viewer (Maan and Cummings 2012). To further reduce the dimensionality of the quantum catch data, we performed a principal components analysis using the prcomp function in R, centering and scaling the data. We retained the first two principal components for further analysis. Pigment Analyses We analyzed skin pigments (carotenoids and pterins) using absorbance spectrometry (Supplementary Fig. S1) and thin-layer chromatography (TLC), following Twomey et al. (2020a, 2020b). These approaches yield complementary data sets with different strengths and weaknesses. Absorbance spectrometry of bulk pigment extracts gives a good measure of overall pigment concentrations and their color, as the processing of extracts is absolutely minimized; however, it provides limited information on pigment identity. TLC, on the other hand, allows for identifications of specific pigments but is far more labor intensive and involves much more sample processing, potentially leading to a loss of quantitative information. Therefore, for the purposes of this study, quantitative pigment information was taken from absorbance spectra while TLC was used only for qualitative determination of pigment presence/ absence. Details on TLC of carotenoids and pterins are given in the Supplementary Material. We identified individual carotenoid and pterin pigments following the methods described in Twomey et al. (2020a, 2020b). Absorbance spectra of carotenoid and pterin extracts were plotted with Pavo2 (Maia et al. 2019). Carotenoid and pterin concentrations were determined from absorbance spectra using the Beer-Lambert law (Hudon and Brush 1992) as detailed in the Supplementary Material. Our pigment analyses resulted in two quantitative traits derived from absorbance spectra (carotenoid concentration at λmax and pterin concentration at 500 nm, which corresponds to the drosopterin absorbance peak), and 7 discrete traits (5 carotenoids and two pterins) representing all skin pigments with a color in the visible range (presence/absence of the carotenoids beta-carotene, lutein, C-4 ketocarotenoids, canary xanthophyll, eschscholtzxanthin-like, and the pterins drosopterin and sepiapterin) derived from the TLC results. These were included in our phylogenetic comparative analysis. Phylogenetics We sampled tissues from as many species and color morphs as possible in Ranitomeya to infer a phylogeny for use in our comparative analysis of color evolution. Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 SYSTEMATIC BIOLOGY 1250 We used ultraconserved elements (UCEs), which are genome-wide molecular markers containing an “ultraconserved” core region with increasingly variable flanking regions, facilitating the identification and alignment of homologous sequences (Faircloth et al. 2012). Such markers have been used to resolve phylogenies at both deep and shallow timescales, including dendrobatid poison frogs (Crawford et al. 2012; Smith et al. 2013; Muell et al. 2022). We were able to include all 16 currently recognized species of Ranitomeya in the phylogeny, including R. defleri and R. ventrimaculata (for which we had no pigment or color data), plus two potentially undescribed species (Fig. 1, “sp. Catuaba” and “sp. Envira”). Three species of Ranitomeya were represented by multiple color morphs (R. variabilis, R. sirensis, and R. imitator; Fig. 1), resulting in a total of 23 terminals in our phylogenetic ingroup. Given the potentially strong effect of outgroup sampling on ingroup topology (Grant 2019), we included many outgroup samples from the subfamily Dendrobatinae, with at least one species per genus, except for “Colostethus” ruthveni (Supplementary Fig. S2). Genomic DNA was extracted from Brazilian tissue samples with the Promega Wizard extraction kit and from Peruvian tissue samples with the Qiagen DNeasy extraction kit. We sent extracted DNA to RAPiD Genomics (Gainesville, FL, USA), where library preparation and Illumina sequencing of UCEs was done following Faircloth et al. (2012). The Tetrapod UCE 5Kv1 probe set was used to target 5060 UCE loci. Roughly 1–2 million paired-end 150 bp reads were sequenced per sample. We used Phyluce 1.6 (Faircloth 2016) to clean sequencing reads, assemble reads into contigs, and extract and align UCE loci for phylogenetic analyses as detailed in the Supplementary Material. We also used Phyluce to create a 95% complete data matrix, that is, one where only loci present in 95% or more of taxa were retained. Aligned loci were concatenated into a single alignment for phylogenetic analysis. We searched for the best partition scheme and model of nucleotide evolution using ModelFinder (Kalyaanamoorthy et al. 2017), with the commands -m testnewmergeonly and -rcluster 30 and the Bayesian information criterion (Schwarz 1978). We performed 200 independent tree searches under the maximum likelihood (ML) criterion in IQ-TREE 2.0.6 (Minh et al. 2020b), with the command -nstop set to 500, the option -allnni activated, and using the models and partitions resulting from ModelFinder (Chernomor et al. 2016). We used default settings for other parameters. We also calculated 1000 replicates of the ultrafast bootstrap approximation (UFBOOT) using IQ-TREE with default settings and the partitions selected by ModelFinder (Minh et al. 2013; Hoang et al. 2018). We used the SumTrees package in DendroPy v4.4.0 (Sukumaran and Holder 2010) to add UFBOOT frequencies to the best ML tree. To quantify genealogical concordance, we calculated gene and site concordance factors (gCF and sGF, respectively) using IQ-TREE (Minh et al. 2020a). To assess phylogenetic resolution, we collapsed nodes incompatible with those from slightly suboptimal trees, defined as the trees from our 200 independent searches that differ from the best reported tree by < 0.1 log-likeli- hood unit. This value is typically an order of magnitude less than the change in log-likelihood after adding a single autapomorphy to the data matrix (Simmons and Kessenich 2020). The UCE sequence alignment, partition scheme, and phylogenetic trees are available on Dryad at https:// doi.org/10.5061/dryad.gf1vhhmsn Comparative Analyses To examine how frog color and coloration mechanisms evolved within Ranitomeya, we used phylogenetic comparative methods implemented in the R package phytools version 0.6-99 (Revell 2012). We pruned the full tree (Supplementary Fig. S2) using the drop.tip function of phytools to remove terminals lacking color/pigment data (mostly outgroup taxa) and made our phylogeny ultrametric with an arbitrary calibration (total tree age = 1) using the chronos function of the R package ape version 0.54 (Paradis and Schliep 2019). We compared ultrametric phylogenies created with each of the 3 available models of substitution rate variation (correlated, relaxed, and discrete), two levels of the lambda parameter (0 and 1), and 1–5 rate categories for the discrete model (nb.rate. cat parameter). For further analyses, we selected the ultrametric tree with the lowest phylogenetic information criterion (PHIIC or ΦIC; Paradis 2013), which was the discrete model with a single rate, equating to a strict-clock. All subsequent analyses were based on this tree. We explored the mode of trait evolution by fitting Brownian motion (BM), Ornstein-Uhlenbeck (OU), and white-noise models to continuous color traits using the fitContinuous function of the R package geiger version 2.0.6.4 (Pennell et al. 2014), implemented in phytools. For each trait, we evaluated model support by comparing AIC weights among models. However, given the criticisms of AIC for phylogenetic model selection, we also performed Phylogenetic Monte Carlo (PMC) analyses, implemented in the R package pmc, to assess model support as well as statistical power to discriminate among competing models (Boettiger et al. 2012). PMC analyses were run with 500 bootstrap replicates and by testing all 3 possible pairwise model comparisons. Ancestral state reconstruction for continuous traits was done with the contMap function of phytools using likelihood estimation (anc.ML), according to their bestfit model of trait evolution. For the dorsal color metric and pterin concentration, this was Brownian motion (Supplementary Table S3; Supplementary Fig. S3). For carotenoid concentration and reflectance λmax, this was OU, which was marginally favored over white noise in the PMC analyses (Supplementary Fig. S3). VOL. 72 Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 TWOMEY ET AL. - COLOR CONVERGENCE IN POISON FROGS2023 1251 Ancestral state reconstruction for discrete traits was done with stochastic character mapping (make. simmap function in phytools), an equal rates model (model = ER), and 1000 simulations. We found that the ER model gave much more parsimonious reconstructions (minimizing the number of state changes) than an all-rates-different model, and thus used ER for each trait. The traits used in these analyses were presence/ absence of different pigments based on the TLC data. When reporting numbers of state changes, we give the mode (most frequent number of state changes), in addition to the mean and standard errors, which were calculated from the 1000 simulations of the stochastic character mapping. We calculated phylogenetic signal for continuous traits using Blomberg’s K (Blomberg et al. 2003), implemented in phytools with the phylosig function, with 10,000 simulations. For discrete traits, we used the method of Fritz and Purvis (2010). This method calculates a metric, D, which can take negative values when traits are clumped (e.g., a trait evolving once and being retained in all descendants), values around zero when consistent with Brownian motion, and values above 1 when random or overdispersed. We calculated D using the phylo.d function of caper version 1.0.1 with P-values based on 10,000 permutations (Orme et al. 2013). Results Phylogeny Across all 38 samples included in the phylogenetic analysis, we recovered a total of 3258 UCE loci. After sequence trimming and filtering for completeness (keeping only loci present in 95% or more of taxa), the data set was reduced to 1550 UCE loci. Our final alignment consisted of 786,579 characters, of which 57,236 were parsimony-informative. Average missing data per sample was 4.95%. ModelFinder selected 49 partitions, each including from 1 to 80 loci. Partitions and their corresponding models of nucleotide evolution are given in Supplementary Table S4. The optimal phylogenetic tree had a log-likelihood of −2268095.948. We found 26 slightly suboptimal trees (i.e., a difference < 0.1 log-likelihood unit with the optimal tree) but their topologies were identical to the optimal tree (only the branch lengths differed). As expected with our ratio of informative characters to nodes, UFBOOT values are high. In general, concordance factors are medium to high, but with low values associated with some of the more recent splits such as the most exclusive clade including R. defleri and R. summersi. Consistent with previous studies (Grant et al. 2006, 2017; Santos et al. 2009; Brown et al. 2011; Muell et al. 2022), our phylogenetic analysis of Ranitomeya based on UCEs recovered two major clades, herein labeled the vanzolinii clade and the fantastica clade (Fig. 1). For an updated infrageneric taxonomy, we refer the reader to Muell et al. (2022). We also note that, based on their phylogenetic position, sampling location, and morphology, it is currently not possible to assign the Catuaba and Envira populations to any known species, and are thus referred to as “sp.” herein. Ranitomeya cf. yavaricola bears morphological similarity to bona fide R. yavaricola, but as this population is ~300 km S from the type locality, we use the label cf. to express uncertainty regarding the species designation. Despite our inclusion of a single sample per taxon, our phylogenetic results are similar to those of Muell et al. (2022). Most differences are associated with different sampling within species, such as the inclusion of additional sampling localities (e.g., R. sirensis and R. uakarii). However, the relative position of R. defleri and R. toraro deserves further comparisons. When the UCE data set was analyzed allowing for different models of evolution for different parts of the data, these two species were nested within the fantastica clade and sister to the reticulata group sensu Muell et al. (2022) (BS = 100%, Fig. 1, Supplementary Fig. S2). On the other hand, analyses applying a single model of evolution for all the data (Supplementary Fig. S4) or a quartet-based approach (Muell et al. 2022) recovered R. toraro as sister of all the other species of the fantastica clade (BS = 100%). We present and discuss our results on the basis of the former topology as it represents a more intensive tree search and considers that different parts of the genome can change under different models of nucleotide evolution. However, we also performed all downstream analyses using the topology shown in Supplementary Figure S4 and the main results of our study did not change. Dorsal Color Evolution Our PCA of dorsal reflectance spectra yielded two main principal components (eigenvalues > 1). Principal component 1 explained 58.2% of the variation, loading positively on quantum catches for S (short-wavelength cone, loading = 0.549) and M (medium-wavelength cone, loading = 0.517), and negatively on L (long-wave- length cone, loading = −0.653). Biologically, this means that most color variation involves correlated changes in the short to medium wavelength range, with opposite changes in the long wavelength range (e.g., blue frogs with high reflectance in the blue-green range and low in the red; or red frogs with high reflectance in the red range and low in the blue-green). The second principal component explained 27.6% of the variation, loading heavily on U (ultraviolet cone, loading 0.908). For further analysis, we analyzed only PC1, which we refer to as the dorsal color metric (Fig. 2a). More generally, most taxa within Ranitomeya were found to have yellowish-orange dorsal coloration; our analyses further indicate that yellowish-orange dorsal coloration was the ancestral state both for Ranitomeya overall as well as the fantastica and vanzolinii clades Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 SYSTEMATIC BIOLOGY 1252 (Figs 1 and 2a). In this context, color transitions and convergence can be more clearly defined as they generally involve the transition from yellowish-orange to some other color. Bluish-green dorsal coloration evolved once within Ranitomeya in the flavovittata clade (Fig. 1). Within this clade, R. flavovittata regained yellow coloration from a blue/green ancestral state. Other gains of yellow or greenish coloration (e.g., the spotted morphs of R. variabilis and R. imitator and the yellow morph of R. sirensis) involve a color shift in the opposite direction, that is, yellow or greenish coloration evolving from an orange ancestral state. Our analyses reveal at least two independent gains of red coloration: once in the reticulata clade, and once in the red morph of R. sirensis. Other gains of reddish coloration are more subtle in that the species did not evolve red coloration per se, but rather evolved redder coloration compared to the inferred ancestral state. For example, in R. amazonica and the banded and varadero morphs of R. imitator, each has a deeper reddish-orange coloration compared to their sister taxa as well as their ancestral state (Fig. 1). Carotenoids Chloroform skin extracts, which contain carotenoid pigments, were generally yellow across most samples. Absorbance spectrometry of these extracts confirmed the presence of carotenoids, which usually have characteristic triple-peak absorbance spectra and absorbance maxima typically in the range of 450–460 nm in the case of chloroform extracts (Supplementary Fig. Figure 1. Color evolution in Ranitomeya. On the left, a maximum-likelihood ultrametric phylogeny of Ranitomeya based on 1550 ultraconserved element loci. Branches are colored according to ancestral state reconstructions of dorsal coloration. On the right, boxes show dorsal color, pigment colors, and reflectance spectra. Dorsal color: Red-Green-Blue (RGB) color measured from dorsal photographs (i.e., the observed frog color that was mapped onto the phylogeny). Carotenoids: RGB rendering of the chloroform extract, modeled onto a white reflector. Pterins: RGB rendering of the ammonia extract, modeled onto a white reflector. Refl: Dorsal reflectance spectrum; numbers indicate the wavelength of peak reflectance (reflectance λmax). RGB renderings, reflectance spectra, and reflectance λmax values represent population averages; in the case of reflectance spectra the solid red line is the mean and red shaded area is the standard deviation. RGB renderings of pigment absorbance spectra were done with the spec2rgb function of pavo2 (see Supplementary Material for more details). For complete phylogeny with outgroups see Supplementary Figure S1. VOL. 72 Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 TWOMEY ET AL. - COLOR CONVERGENCE IN POISON FROGS2023 1253 Figure 2. Ancestral state reconstructions for 4 continuous traits. (a) Dorsal color metric PC1 (58.2% of variation), derived from reflectance spectra, (b) Carotenoid concentration, derived from chloroform absorbance spectra, (c) Pterin concentration at 500 nm, derived from ammonia absorbance spectra, and (d) Wavelength of peak reflectance (reflectance λmax), derived from dorsal reflectance spectra. Note that the pterin concentration is based on absorbance at 500 nm and is thus indicative of drosopterin concentration. For b, c, and d, the color scales are representative of the color gamut controlled by that trait; for a, the color scale is a standard blue-to-red color ramp which approximately coincides with the observed colors of these frogs. The rectangles at the tips show the RGB dorsal color of that taxon. Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 SYSTEMATIC BIOLOGY 1254 S1). Carotenoids were further confirmed with TLC (see below). Although most samples contained carotenoids based on the chloroform absorbance spectra, there were some notable exceptions. Within the vanzolinii clade, carotenoid reduction or loss was observed in taxa with blue or bluish-green coloration, that is, Ranitomeya cyanovittata, R. cf. yavaricola, and R. sp. Envira (Fig. 2b, Supplementary Fig. S5). In the fantastica clade, carotenoid loss or reduction was found in taxa with red coloration (Fig. 2b, Supplementary Fig. S5). Ranitomeya benedicta, which has vivid red coloration on the head, had no detectable carotenoids from the chloroform absorbance spectra (although trace beta-carotene was found with TLC; see below), while R. reticulata and R. uakarii, which are also red, had very low amounts of carotenoids (Supplementary Fig. S5). By contrast, the red morph of R. sirensis, which is the only red taxon in the vanzolinii clade, had a distinctly red chloroform extract due to its high carotenoid concentrations, including red C-4 ketocarotenoids (Twomey et al. 2020a). Ancestral state reconstruction of carotenoid concentration (Fig. 2b) suggests that the ancestral state in Ranitomeya was an intermediate carotenoid concentration, with two subsequent reductions: once in the reticulata clade and once in the flavovittata clade. Carotenoid concentrations showed independent increases from a lower ancestral state in R. summersi, R. flavovittata, the banded morph of R. imitator, and, to a lesser extent, in the most recent common ancestor of R. vanzolinii and R. sirensis. Using TLC, we were able to identify a number of the skin carotenoids. The most widespread were the dietary carotenoids beta-carotene and lutein (Fig. 3). Beta-carotene was present in all taxa except for R. cf. yavaricola, R. cyanovittata, and R. sp. Envira, and stochastic character mapping indicated presence at the root (posterior probability = 0.94), followed by a single loss (mean = 1.4 ± 0.03) and single gain (mean = 1.4 ± 0.03). Lutein was also absent in these 3 taxa, as well as R. benedicta and R. reticulata, and most likely present at the root (posterior probability = 0.63), followed by 4 losses (mean = 4.6 ± 0.05) and 4 gains (mean = 4.6 ± 0.07), making it the most evolutionarily labile skin pigment (Fig. 3). Red C-4 ketocarotenoids have been found previously in R. sirensis and R. summersi (Twomey et al. 2020a, 2020b). We did not find C-4 ketocarotenoids in any other sample and thus their presence remains restricted to these two species. Character mapping indicates absence at the root (posterior probability = 0.97), with two independent gains (mean = 2.1 ± 0.01) in R. summersi and the red morph of R. sirensis, and no losses (mean = 0.3 ± 0.02). Results from other carotenoids are given in Supplementary Fig. S4. Phylogenetic signal of categorical traits (D) was negative for beta-carotene (D = −0.54), canary xanthophyll (D = −1.4), and drosopterin (D = −0.14), consistent with “clumped” trait distributions where traits evolve few times and are retained in descendants; similarly, P-values for these 3 pigments rejected the possibility of random trait evolution. In the other pigments, D was positive and consistent with a random phylogenetic pattern. Pterins Ammonia-based skin extracts, which contain pterin pigments, were mostly colorless or pale yellow in species of the vanzolinii clade (Fig. 1). However, all the species in the fantastica clade that had orange or red dorsal coloration also had orange/red pterin extracts (Fig. 1) with absorbance peaks or shoulders around 500 nm (Supplementary Fig. S1b), characteristic of the pterin pigment drosopterin. TLC of pterin extracts confirmed the presence of large quantities of drosopterin in these species, in particular R. reticulata and R. benedicta (Supplementary Fig. S5). In the vanzolinii clade, drosopterin is known from R. imitator (Twomey et al. 2020b), and we did not detect it in any of the other species in this clade. Stochastic character mapping of drosopterin reveals a slightly higher probability of absence at the root (posterior probability 0.58), with 3 gains (mean = 2.9 ± 0.04) and 3 losses (mean = 3.3 ± 0.05); notable putative gains involve at the base of the reticulata clade, and in the species R. imitator and R. amazonica (Fig. 3). For results on other pterin pigments, see Supplementary Material. Reflectance Spectra and Structural Coloration All dorsal reflectance spectra, except for Ranitomeya sp. Catuaba, had a hump-shaped reflectance curve (Fig. 1) that is suggestive of a color contribution from structurally colored iridophores. Coherently scattering iridophores are predicted to yield similar hump-shaped reflectance curves where the location of the spectral peak is controlled in large part by the thickness of the iridophore platelets (Morrison 1995; Saenko et al. 2013). Reflectance curves based on purely structural coloration are expected to be roughly symmetrical around the peak, although here we observed that reflectance spectra typically have a steeper slope left of the peak, likely due to carotenoid and pterin pigments absorbing these shorter wavelengths. In some taxa, such as R. cf. yavaricola, R. sp. Envira, and R. cyanovittata, which largely lack carotenoid and pterin pigments (Supplementary Fig. S1), structurally colored iridophores represent the most plausible mechanism generating the observed colors. Similarly, R. toraro was found to have a relatively small color contribution from carotenoid and pterin pigments, yet this species has a bright yellow dorsal coloration and peaked reflectance spectrum (Fig. 1); Twomey et al. (2020b) also reported that some individuals of R. imitator also had bright yellow or orange dorsal coloration despite low concentrations of carotenoid and pterin pigments. Ranitomeya sp. Catuaba is the lone species whose dorsal reflectance spectrum is not distinctly humpshaped, instead sloping up continuously past the edge of the measurement range, unlike any of the other species (Fig. 1). This raises the possibility that structural VOL. 72 Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024 TWOMEY ET AL. - COLOR CONVERGENCE IN POISON FROGS2023 1255 Figure 3. Ancestral state reconstructions for presence/absence of carotenoid and pterin pigments in Ranitomeya. (A) Beta carotene (carotenoid), (B) Lutein (carotenoid), (C) C-4 ketocarotenoids (carotenoid), and (D) Drosopterin (pterin). Pigments were scored as present/ absent based on TLC results and absorbance spectra. Pigment color (red/yellow) is indicated with the label color. Observed states are shown at the tips; pie charts at internal nodes represent posterior probabilities of each state. Boxes report the mean ± standard error for the numbers of total changes, gains, and losses calculated from the 1000 stochastic mapping iterations. For 3 taxa (R. cf. yavaricola, R. sp. Envira, and R. cyanovittata), pigments were scored as absent based on absorbance spectra that indicated a lack of carotenoids and drosopterin. Other species that had skin carotenoids but for which we did not have TLC data (R. vanzolinii, R. toraro, and R. sp. Catuaba) were omitted from this analysis. Results from 3 additional pigments are shown on Supplementary Figure S6. Downloaded from https://academic.oup.com/sysbio/article/72/6/1247/7240510 by Sevilla Psych Lib B21 user on 08 May 2024