Radiation and hybridization of the Little Devil poison frog (Oophaga sylvatica) in Ecuador

Geographic variation of color pattern in the South American poison frogs (Dendrobatidae) is an intriguing evolutionary phenomenon. These chemically defended anurans use bright aposematic colors to warn potential predators of their unpalatibility. However, aposematic signals are frequency-dependent and individuals deviating from a local model are at a higher risk of predation. The well-known examples of Batesian and Müllerian mimics, hymenopterans (wasps and bees) and Heliconius butterflies, both support the benefits of unique models with relatively high frequencies. However, extreme diversity in the aposematic signal has been documented in the poison frogs of the genus Dendrobates, especially in the Oophaga subgenus. Here we investigate the phylogenetic and genomic differentiations among populations of Oophaga sylvatica, which exhibit one of the highest phenotypic diversification among poison frogs. Using a combination of PCR amplicons (mitochondrial and nuclear markers) and genome wide markers from a double-digested RAD data set, we characterize 13 populations (12 monotypic and 1 polytypic) across the O. sylvatica distribution. These populations are mostly separated in two lineages distributed in the Northern and the Southern part of their range in Ecuador. We found relatively low genetic differentiation within each lineage, despite considerable phenotypic variation, and evidence suggesting ongoing gene flow and genetic admixture among some populations of the Northern lineage. Overall these data suggest that phenotypic diversification and novelty in aposematic coloration can arise in secondary contact zones even in systems where phenotypes are subject to strong stabilizing selection.


INTRODUCTION
Aposematism is an antipredator adaptation that has evolved in many animals as a defense strategy. This adaptation combines warning signals (e.g., vivid coloration) with diverse predator deterrents such as toxins, venoms and other noxious substances. Most groups of animals include at least one aposematic lineage such as insects (e.g., Heliconius (Jiggins & McMillan 1997) and Monarch (Reichstein et al. 1968) butterflies), marine gastropods (e.g., nudibranchs (Tullrot & Sundberg 1991;Haber et al. 2010), amphibians (Howard & Brodie 1973;Brodie 2008), reptiles (Brodie III 1993;Brodie III & Janzen 1995) and birds (Dumbacher et al. 1992(Dumbacher et al. , 2000. The aposematism strategy is dependent on the predictability of the warning signal for effective recognition as well as learning and avoidance by predators (Benson 1971;Mallet et al. 1989;Kapan 2001;Pinheiro 2003;Ruxton et al. 2004;Chouteau et al. 2016). However, some aposematic species present high levels of polymorphic coloration and patterning (Przeczek et al. 2008). The causes and consequences of intraspecific variation in warning signals remains a fundamental question in the evolution and ecology of aposematism.
Dendrobatid poison frogs are endemic to Central and South America and have evolved chemical defenses coupled with warning coloration 3-4 times within the clade (Santos et al. 2003(Santos et al. , 2009(Santos et al. , 2014. Interestingly, some poison frog genera/subgenera display a wide range of inter-and intra-specific color and pattern variability, including Dendrobates sensu lato (Noonan & Wray 2006;Noonan & Gaucher 2006;Comeault & Noonan 2011), Adelphobates (Hoogmoed & Avila-Pires 2012), Ranitomeya (Symula et al. 2001;Twomey et al. 2013Twomey et al. , 2014Twomey et al. , 2016 and Oophaga (Wollenberg et al. 2008;Wang & Summers 2010;Brusa et al. 2013;Medina et al. 2013). For example, some species like R. imitator have evolved variation in coloration and patterning in a unique vertebrate example of Müllerian mimicry (Symula et al. 2001). Another example is the Strawberry poison frog, Oophaga pumilio, which is extremely polymorphic within the geographically small archipelago of Bocas del Toro (Panamá), but has more continuous and less variable coloration in its mainland distribution from Nicaragua to western Panamá. Many factors might account for the origin and maintenance of O. pumilio color polymorphism including selective pressure from multiple predators, mate choice based on visual cues, and genetic drift among the islands of the archipelago (Gehara et al. 2013). However, it is unclear how the subgenus Oophaga as a whole has evolved such extreme diversity in coloration and patterning without evident geographic barriers.
Molecular phylogenies show that Oophaga as lineage includes many cryptic species that underlie its diversity. For example, O. pumilio populations include two distinctive mitochondrial lineages, each of which contains at least one closely related congeneric species O. speciosa, O. arborea or O. vicentei (Hagemann & Prohl 2007;Wang & Schaffer 2008;Hauswaldt et al, 2011). Mitochondrial haplotypes from both lineages are found co-occurring across a wide zone along the Panamá and Costa Rica border. These populations present variable levels of admixture that also might account for their phenotypic diversity. This complex phylogeographic pattern suggests a series of dispersals and isolations leading to allopatric divergence and then subsequent admixture and introgression among O. pumilio and other species of Oophaga (Hagemann & Pröhl 2007, Wang & Shaffer 2008, Hauswaldt et al. 2011 (Myers & Daly 1976). Oophaga histrionica has long been suspected to be a species complex (Silverstone 1975;Lotters et al. 1999), and has recently been proposed to include three new species evolving along a diversifying landscape (Posso-Terranova & Andrés 2016a). However, these new species descriptions need further morphological, behavioral, toxicological, and genomic characterization to be clearly validated. The Oophaga subgenus is composed of nine species, which have extraordinary morphological and chemical diversity (Daly & Myers 1967;Daly et al. 1978;Daly 1995;Saporito et al. 2007a). Most research has been conducted in O. pumilio, including a number of studies in population genetics, phylogeography, behavior, diet specialization, and chemical defenses (Saporito et al. 2007b;a;Richards-Zawacki et al. 2012;Gehara et al. 2013;Stynoski et al. 2014;Dreher et al. 2015). Comparative work in other Oophaga species would lend insight into interesting speciation and diversification processes.
We explored population structure of the highly polymorphic Little Devil (or Diablito) poison frog, O. sylvatica, which has an inland distribution in the Chocoan rainforest and ranges from southwestern Colombia to northwestern Ecuador (Funkhouser 1956). Two types of population structures exist within this taxon: one composed of many phenotypically distinctive populations with low within-population color and pattern variability, and another structure containing a unique population with extreme withinpopulation variability. The color and pattern observed in this latter population may represent either an admixture or hybridization zone, with individuals showing intermediate color patterns that combine traits from the surrounding monotypic populations. Previous research in other Oophaga (i.e., O. histrionica and O. lehmanni) suggests that hybridization may promote color polymorphism (Medina et al. 2013). However, whether such similar admixture or hybridization mechanisms also promote color polymorphisms within O. sylvatica is unknown, as genetic studies have never been conducted with this species.
To investigate the population structure of O. sylvatica, we used a set of PCR amplicons including mitochondrial and nuclear markers, and a collection of genome wide single nucleotide polymorphisms (double digested RAD sequencing) from 13 geographically distinct populations (12 low phenotypic diversity and 1 highly polymorphic). The aims of our research were: 1) explore the relationship between individuals and their geographic localization, 2) estimate population structure and differentiation and 3) evaluate if admixture or hybridization could account for color diversity within and among populations.

Sample collection
We sampled 13 populations of Oophaga sylvatica (N = 200 individuals) in July 2013 and July 2014 throughout its Ecuadorian distribution ( Figure 1). Frogs were collected during the day with the aid of a plastic cup and stored individually in plastic bags with air and leaf litter for 3-8 h. In the evening the same day of capture, individual frogs were photographed in a transparent plastic box over a white background. Frogs were anesthetized with a topical application of 20% benzocaine to the ventral belly and euthanized by cervical transection. Tissues were preserved in either RNAlater (Life Technologies, Carlsbad, CA, USA) or 100% ethanol. Muscle and skeletal tissue were deposited in the amphibian collection of Centro Jambatu de Investigación y Conservación de Anfibios in Quito, Ecuador (Suppl. were also treated using the same protocol. In order to protect the vulnerable O. sylvatica populations that are highly targeted by illegal poaching, specific GPS coordinates of frog collection sites can be obtained from the corresponding author.

DNA extraction and amplification
Liver or skin tissue stored in RNAlater was homogenized in using Trizol (Life Technologies) in tubes with 1.5 mm TriplePure Zirconium Beads (Bioexpress, Kaysville, UT, USA).
After tissue homogenization, DNA was purified by a standard phenol/chloroform procedure followed by ethanol precipitation according to the Trizol manufacturer's instructions.
Mitochondrial genes (12S-tRNA Val , 16S, CO1) were considered as a single unit and concatenated in a matrix of 2084 bp for each individual of the different populations of O. sylvatica (N=199, one sample was excluded as it failed to amplify for 12S-tRNA Val ), and for O. histrionica (N=2) and O. pumilio (N=6) individuals. Diversity indices were calculated using DnaSP v5.10.01 and Arlequin 3.5 (Excoffier & Lischer 2010). Estimation of genetic distances between mtDNA sequences was tested under different models for nucleotide substitution and the HKY+I model was selected according to the Bayesian information criterion implemented in jModeltest2 (Guindon & Gascuel 2003;Darriba et al. 2012). We retained the Tajima and Nei model (Tajima & Nei 1984) to infer the haplotype network in Arlequin 3.5, which was a representative consensus of the results obtained with the other models tested. Minimum spanning networks (MSN) were generated using Gephi (Bastian et al. O. sylvatica were found in lowland and foothill rainforest (0 to 1,020m above sea level) in northwestern Ecuador. Most frogs were phenotypically variable among geographical localities (populations), while relatively stable within populations (Suppl. Figure 1). Color diversity is particularly dramatic, ranging from yellow to red to brown and greenish and can be combined with either markings or spots of different colors. B. A striking example of diversity within the population of Otokiki, located in the center of the northern range, with phenotypes similar to the surrounding monomorphic populations as well as intermediate phenotypes.
We assessed population differentiation by calculating the conventional F-statistics from haplotype frequencies (F ST ) and the statistics from haplotype genetic distances based on pairwise difference (Φ ST ) using Arlequin 3.5. P values for F ST and Φ ST were estimated after 10,000 permutations and significance threshold level was fixed at p = 0.05.
We performed a Bayesian assignment test on the nuclear dataset to estimate the number of genetic clusters as implemented in the program STRUCTURE 2.3.4. The software assumes a model with K populations (where K is initially unknown), and individuals are then assigned probabilistically to one or more populations (Pritchard et al. 2000). We ran the admixture model with correlated allele frequencies (Falush et al. 2003) for 50,000 burn-in and 100,000 sampling generations for K ranging from one to the number of sampled population plus three (K = 1-16) with 20 iterations for each value of K. We determined the number of clusters (K) that best described the data using the delta K method (Evanno et al. 2005) as implemented in Structure Harvester (Earl & vonHoldt 2012), which was inferred at K=5. Then we ran the same model for 500,000 burn-in generations and 750,000 MCMC, from two to six clusters (K = 2-6) and 20 iterations and analyzed the results using CLUMPAK (Kopelman et al. 2015).

ddRADseq library generation and sequencing
We constructed double-digested restrictionsite-associated DNA sequencing (ddRAD) libraries on a subset of 125 samples following the protocol in Peterson et al. (2012). Samples include three specimens randomly drawn within each sampling site from the monomorphic populations and all the specimens from the putative admixture zone (Otokiki locality, see Figure 1). DNA was extracted from skin tissues preserved at -20°C in RNAlater using the Nucleospin DNA kit (Macherey-Nagel, Bethlehem, PA). Genomic DNA of each sample (1 µg) was digested using 1 µL of EcoRI-HF (20,000 U/mL) and 1 µL SphI-HF (20,000 U/mL) (New England Biolabs, Ipswitch, MA) following the manufacturer's protocol. Digested samples were then cleaned with Agencourt Ampure XP beads (Beckman Coulter, Danvers, MA). Purified digested DNA (100ng) was ligated to double stranded adapters (biotin-labeled on P2 adapter) with a unique inline barcode using T4 DNA ligase (New England Biolabs, Ipswitch, MA) and purified with Agencourt Ampure XP beads. Barcoded  (SRP078453).

ddRAD sequence analysis
Raw fastq reads were demultiplexed, quality filtered to discard reads with Phred quality score < 20 within a sliding window of 15% of read length and trimmed to 120bp using the process_radtags.pl command from the Stacks v1.35 pipeline (Catchen et al. 2013). ddRAD loci were constructed de novo with the Stacks denovo_map pipeline (parameters: m = 5, M = 3, n = 3) and then corrected for misassembled loci from sequencing errors using the corrections module (rxstacks, cstacks and sstacks), which applies population-based corrections. In order to minimize the effect of allele dropout that generally leads to over estimation of genetic variation (Gautier et al. 2013), we selected loci present in at least 75% of the individuals of each of the 13 sampled populations and generated population statistics and output files using Stacks population pipeline (parameters: r = 0.75, p = 13, m = 5).
To test for admixture, we estimated the number of genetic clusters as implemented in the program STRUCTURE 2.3.4 (Pritchard et al. 2000). The data set was filtered using the Stacks population pipeline (parameters: r = 0.75, p = 13, m = 5) to contain only loci present in at least 75% of each of the 13 sampled localities. We retained only one SNP per RAD locus to minimize within-locus linkage. The data set contains 3,785 SNPs from concatenated pairended reads. We ran the admixture model with correlated allele frequencies (Falush et al. 2003) for 100,000 burn-in generations and 1,000,000 MCMC, for K = 2-16 with 10 iterations for each cluster. We determined the number of clusters K that best described the data (Evanno et al. 2005) as implemented in Structure Harvester (Earl & vonHoldt 2012) and analyzed the results using CLUMPAK (Kopelman et al. 2015). We used a Discriminant Analysis of Principal Components (DAPC) to identify clusters of genetically related individuals, implemented in the R package adegenet v2.0.1 (Jombart et al. 2010;Jombart & Ahmed 2011). This multivariate method is suitable for analyzing large number of SNPs, providing assignment of individuals to groups and a visual assessment of between-population differentiation. It does not rely on any particular population genetics model and DAPC is free of assumption about Hardy-Weinberg equilibrium or linkage equilibrium (Jombart et al. 2010). The data set was built with the Stacks populations pipeline and is composed of 3,785 SNPs (5955 alleles) present in at least 75 % of the individuals from each of the 13 populations. To avoid over-fitting of the discriminant functions we performed stratified cross-validation of DAPC using the function xvalDapc from adegenet and retained 20 principal components, which gave the highest mean success and the lowest Root Mean Squared Error (MSE). We also generated consensus sequences of the 125 individuals grouped by geographical sampling localities to assess their relationship at the population level. We used Bayesian inferences as implemented in Mr Bayes 3.2.6 (Ronquist et al. 2012) under Geneious to build a phylogenetic tree from a matrix of 41,779 SNPs (fixed within but variable among populations), with O. pumilio and O. histrionica as an outgroup. We used the default settings under the HKY85 + gamma model and ran a single MCMC with four chains (0.02 heated chain temp) for 1,100,000 generations, out of which the first 100,000 were discarded as burn-in.

Mitochondrial diversity and structure
After alignment of all 207 individuals, sequences show 163 segregating sites (S) that compose 66 haplotypes (H) ( Table 1 and Figure 2). Among the 13 populations of O. sylvatica sampled, genetic diversity is relatively high, with 83 segregating sites, 61 haplotypes, and a mean haplotype diversity (Hd) of 0.948. Both O. histrionica and O. pumilio are well separated from O. sylvatica and link the network by two distinct branches (Figure 2). These branches constitute two remote clusters, one composed of haplotypes from Otokiki and Felfa populations that link to O. histrionica, and one composed of haplotypes from the San Antonio population that link to O. pumilio. Both branches link to the same node, a haplotype of one individual from Durango. Six short branches (2 to 4 inferred mutations) radiate from this central node: two are composed of private haplotypes from the Durango and San Antonio populations, and the other four have mixed origins belonging to the northern populations of Otokiki, Durango, Lita, Alto Tambo and San Antonio. These populations show the highest sequence (K) and nucleotide (π) diversity, as well as the highest haplotype diversity (except San Antonio population for Hd) ( Table 1). One of the mixed clusters includes haplotypes from Durango, Lita, Alto Tambo and Otokiki populations; another includes haplotypes from Alto Tambo and Otokiki, and two different clusters are composed of haplotypes from the Otokiki and Lita populations. All of these clusters have inter-haplotype distances ranging from 1 to 4 inferred mutational steps. With 29 haplotypes, Otokiki is the most genetically diverse population (Table 1 and Figure 2), in addition to exhibiting the most variable color patterning ( Figure 1B). Most haplotypes are unique to the Otokiki population, except for 6 that are shared with the populations of Alto Tambo, Lita and Durango. These are also the geographically the closest populations to Otokiki ( Figure 1A). Genetically, Lita and Alto Tambo seem weakly separated from Otokiki, with the Lita population having no unique haplotype and only two out of five unique haplotypes in the Alto Tambo population. The southern and northern populations are separated by a unique haplotype from the Felfa population. Southern populations show less genetic diversity than northern populations, where 67 individuals collapsed in 15 haplotypes (Figure 2). One branch is composed of a small group of 3 haplotypes from Puerto Quito, Simón Bolívar and Cristóbal Colón. The second branch leads to one haplotype including 38 individuals from 6 populations: Quingüe, Cube, Simón Bolívar, Puerto Quito, Santo Domingo and La Maná. Radiating from this haplotype are 11 unique haplotypes (1 to 4 individuals each) with mostly 1 inferred mutational step. While closely related to the more southern populations, frogs from the Cristóbal Colón population do not share any haplotype with the rest of the southern populations.

Nuclear DNA diversity and population structure
Nuclear sequences (RAG-1, TYR, and NCX) were phased and analyzed separately. Across all three Oophaga species in this study, the NCX gene presents 10 segregating sites and 15 haplotypes, TYR has 9 segregating sites and 9 haplotypes, and RAG-1 has 4 segregating sites and 5 haplotypes (Suppl. Table 3). Among the nuclear genes, NCX is the most variable in O. sylvatica, with 7 segregating sites and 12 haplotypes, whereas RAG-1 has 2 segregating sites and 3 haplotypes, and TYR has only 1 segregating sites and 2 haplotypes. We then combined and phased the three sets of nuclear markers in a unique matrix for every individual (Table 2). This slightly increases the number of segregation sites (S = 10) and haplotypes (H = 16), as well as diversity indices for our O. sylvatica group. Similar to mitochondrial genes, sequence (K) and nucleotide (π) diversity are the highest for northern populations (Durango, San Antonio, Lita, Otokiki, Felfa, and Alto Tambo), but haplotype diversity (Hd) is higher in a subset of the southern group (Quingüe, Puerto Quito, and Santo Domingo). Within the minimum spanning network of haplotypes, inferred distances between nodes are very short, ranging from 1 to 5 mutation steps (Figure 4)

Population structure based on ddRAD markers
A structure plot was drawn using Bayesian inferences with our ddRAD data (a subset of 125 individuals of O. sylvatica). The optimal number of clusters inferred by Evanno's method was K = 3 (see Figure 5 for colors assigned to clusters). In O. sylvatica we can observe a clear genetic structure. One genetic cluster (blue) represents mostly populations from the northernmost part of the range (San Antonio, Lita, Alto Tambo, Durango and Otokiki). A second genetic cluster (orange) appears mostly in populations from the Southern part of the range (Felfa, Cristóbal Colón, Simón Bolívar, Quingüe, Cube,

Discriminant Analysis of Principal Components (DAPC) on the ddRAD data
We used 3,785 SNPs from the ddRAD dataset to conduct a DAPC analysis ( Figure 6). The optimal number of clusters inferred under the Bayesian Information Criterion was K = 2. The two clusters follow the north/south split observed previously. In the north, individuals from the populations of Lita, Alto Tambo and Durango overlap with individuals from Otokiki, but San Antonio is well separated from them. In the south, populations are separated from each other and individuals did not overlap with any other population. By plotting the membership probability of each individual using the compoplot function with two discriminant analysis (DA), we can see that overlapping individuals from Lita, Alto Tambo, Durango and Otokiki, suggestive of admixture ( Figure  6). In addition, some individuals from the southern populations of Cube and La Maná have been assigned with mixed proportions to both populations.

Phylogenetic relationships derived from ddRAD data
Using consensus sequences for each geographical locality, we built a population-based tree using Bayesian inferences and setting O. pumilio as the outgroup species (Figure 7). The San Antonio population, which is located in the northern-most part of the Ecuadorian range, has diverged before the rest of the O. sylvatica Ecuadorian populations. Two major branches split other populations with moderate support (0.93 posterior probability), grouping Alto Tambo, Lita, Durango and Otokiki in one cluster, and Felfa, Cristóbal Colón, Simón Bolívar, Puerto Quito, Cube, Quingüe, Santo Domingo and La Maná in another cluster.

DISCUSSION
In this study we investigated genetic structure among populations of Oophaga sylvatica in Ecuador. Our results show that the genetic structure of the populations follows their geographical distribution and is split in two main lineages, distributed in the northern and the southern part of their range. The existence of two distinct lineages is supported by mitochondrial haplotype network (Figure 2), Bayesian assignment of individuals using Structure ( Figure 5), DAPC analysis ( Figure 6) and a phylogeny (Figure 7) using ddRAD data. Within each lineage, most populations share mitochondrial haplotypes (with the exception of San Antonio, Felfa and Cristóbal Colón) and have low genetic differentiation (low or non significant F ST values). We hypothesize that gene flow from recent range expansions might be driving the color diversity in these frogs.
The northern lineage (San Antonio, Lita, Alto Tambo, Durango and Otokiki) presents higher mitochondrial diversity with a large amount of unique haplotypes within populations. We observe that four of these populations (San Antonio, Lita, Alto Tambo and Durango) have unique haplotypes, and three (Durango, Lita and Alto Tambo) share haplotypes with the polymorphic population at Otokiki. These latter populations are weakly differentiated based on genetic data and are geographically close (distances from Otokiki are approximately 4km, 5km and 20 km to Alto Tambo, Lita and Durango, respectively). Otokiki individuals present all color patterns found in the surrounding populations in addition to intermediate patterns ( Figures 1A and 1B). The DAPC plot highlights a genetic overlap between some individuals of Otokiki and the three populations from Lita, Alto Tambo and Durango; some individuals were assigned with mixed proportions to these different groups  O. histrionica and O. pumilio (2255 bp). Circles indicate haplotypes, with the area being proportional to the number of individuals sharing that haplotype. Colors refer to the geographical origin of the population and the pie charts represent the percentage of each population sharing the same haplotype. Line thickness between haplotypes is proportional to the inferred mutational steps (or haplotypes). Inferred numbers of mutational steps are indicated inside circles along the line when greater than four steps.
( Figure 6). Taken together, our results show a close genetic interaction between populations presenting low diversity and intermediate phenotypes among populations with overlapping distributions. There at least two scenarios that could explain these patterns. First, these features support the presence of an admixture zone occurring within Otokiki, which in turn promotes the dramatic color diversity of this population. Alternatively, Otokiki could be a large source population where surrounding populations represent recent expansions with less genetic and phenotypic diversity due to founder effects or other selection pressures. Although outside the scope of the current study, more analyses on population structure with more sampling in the northern populations is needed to disentangle these alternative hypotheses.
The southern lineage shows little mitochondrial diversity, in which a few variants radiate from a common haplotype shared by a majority of individuals across different populations (Figure 2). This lineage shows great color diversity among populations and little variation within, suggesting that variation in mitochondrial markers is ineffective to resolve the observed phenotypic diversity. Furthermore, these populations are geographically distant from each other (from 20km between Simón Bolívar and Puerto Quito to as much as 190km between Quingüe and La Maná), have low but significant F ST values from ddRAD data (Suppl. Table 4) and are genetically well differentiated on the DAPC analysis. Since amphibians usually have low dispersal abilities (Vences & Wake 2007;Vences & Köhler 2008;Zeisset & Beebee 2008) and Dendrobatidae are not known to be migratory species, the low levels of differentiation observed at the mitochondrial markers is unlikely to occur through gene flow among most of these populations. In addition, this contrasting pattern of low genetic variation but high phenotypic diversity among geographically distant population suggests that the southern clade radiation has probably occurred relatively recently. This genetic pattern may be accentuated by decades of population bottlenecks and genetic drift, as southern populations have been subject to strong isolation due to human intervention and habitat destruction. Indeed, a similar pattern of recent expansion and isolation with reduced gene flow was proposed as a model to explain the phenotypic diversification within the southern lineage of O. pumilio, among the islands of Panamá (Gehara et al. 2013). Repeated phylogeographical patterns of diversification among amphibians and reptiles across Northwestern Ecuador was recently described by Arteaga and colleagues (2016) and suggests that diversification follows thermal elevation gradients among the Chocoan region and the Andes. There are also several natural barriers to dispersion, including the Mira River in the northern most part of the O. sylvatica range. The Mira River geographically separates the San Antonio population from the rest of the northern lineage, which may explain its genetic differentiation. The Esmeraldas River also separates the distribution range of O. sylvatica in two mains zones: the populations of San Antonio, Lita, Alto Tambo, Durango, Otokiki, Felfa and Cristóbal Colón in the north, and the populations of Simón Bolívar, Puerto Quito, Cube, Quingüe, Santo Domingo and La Maná in the south. Arteaga and colleagues (2016) have identified similar geographical patterns of differentiation for other frog species, such as Pristimantis nietoi and P. walkeri, and populations of the snake species Bothrops punctatus near the Esmeraldas and Guayllabamba Rivers (Arteaga et al. 2016). These large rivers may have acted as barriers precluding or reducing gene flow between the northern and southern lineages. Nevertheless, Felfa and Cristóbal Colón populations, distributed north of the Esmeraldas River, present conflicting genetic data suggesting that the barrier did not isolated these two lineages completely. Both populations have admixed proportions of the two lineages with high proportion of the southern lineage ( Figure 5). Still, mitochondrial data support Felfa's relationship to the northern lineage, but Cristobal Colón likely belong to the southern lineage ( Figure 2). Also, this genetic pattern observed for Cristobal Colón could be supported by the presence of the Toisán Mountains extending in this area from west to east, which might have acted as a barrier, isolating the Cristobal Colón population from the rest of the northern lineage. A critical question is raised from this study: How is high phenotypic diversity in coloration and patterning maintained in the Otokiki population? A wide range of putative predators such as birds, reptiles or arthropods with distinct visual abilities and predatory strategies might act on color selection and diversity among O. sylvatica populations (Maan & Cummings 2012;Crothers & Cummings 2013;Dreher et al. 2015). In addition, sexual selection through non-random courtship within color morphs has been reported in O. pumilio (Reynolds & Fitzpatrick 2007;Maan & Cummings 2008) and Ranitomeya (Twomey et al. 2014), and can possibly promote diversification of color among populations, while decreasing the variation within. Although this phenomenon seems highly dependent on the environmental context and genetic background of the frogs (Richards-Zawacki et al. 2012;Meuche et al. 2013;Medina et al. 2013;Twomey et al. 2014). A recent study reported an example of hybridization promoting new coloration and patterning between two close species O. histrionica and O. lehmanni (Medina et al. 2013). This work also suggests complex roles played by sexual selection as hybrid females present non-random sexual preferences depending on the combination of available males. If similar processes occur in O. sylvatica at Otokiki, the extreme polymorphism within this admixture zone represents a unique opportunity to Figure 6. DAPC scatterplot for ddRAD data. The scatterplot shows the first two principal components of the DAPC of data generated with 3,785 SNPs. Geographic samples are represented in different colors and inertia ellipses according to Figure 1, with individuals shown as dots. Lines between groups represent the minimum spanning tree based on the (squared) distances, and show the actual proximities between populations within the entire space. Right inset shows the inference of the number of clusters using the Bayesian information criterion (BIC). The chosen number of clusters corresponds to the inflexion point of the BIC curve (K=2). Left inset shows the number of PCA eigenvalues retained in black, and how much they accounted for variance. The bottom graph represents the membership probability of each individual to one or more populations. Geographic populations are represented with the same colors as in the DAPC plot.
further test the balance of selective pressure on the evolution of an aposematic trait by selection through predators or by conspecifics through mate choice.

Summary
By evaluating mitochondrial DNA variation and genome-wide SNPs, we have gained four important insights about O. sylvatica: (1) The Ecuadorian populations of O. sylvatica are composed of two clades that reflect their geographic distribution. Further behavioral, ecological and morphological information is required to determine whether the two geographical lineages observed in Oophaga sylvatica represent one polytypic or distinct species. A combination of climatic gradient and structured landscape generating geographic barriers to gene flow could explain the complex patterns of diversification observed in Oophaga and some other Dendrobatidae (Arteaga et al. 2016;Posso-Terranova & Andrés 2016ab). (2) The northern and southern populations show different amounts of structure, which may reflect more recent range expansions in the south. (3) Phenotypic variation in O. sylvatica has evolved faster than mitochondrial mutations can fix in the population. (4) A highly polymorphic population (Otokiki) exists, which provides a unique opportunity for testing hypotheses about the selective pressures shaping aposematic traits. We hypothesize that this polymorphic population arose from either gene flow between phenotypically divergent populations at secondary contact zones or through a range expansion of the polymorphic Otokiki population into surrounding regions. More data is needed to distinguish between these alternative scenarios.