One proposed mechanism of speciation is divergent sexual selection, whereby divergence in female preferences and male signals results in behavioural isolation. Despite the appeal of this hypothesis, evidence for it remains inconclusive. Here, we present several lines of evidence that sexual selection is driving behavioural isolation and speciation among populations of an Amazonian frog (Physalaemus petersi). First, sexual selection has promoted divergence in male mating calls and female preferences for calls between neighbouring populations, resulting in strong behavioural isolation. Second, phylogenetic analysis indicates that populations have become fixed for alternative call types several times throughout the species' range, and coalescent analysis rejects genetic drift as a cause for this pattern, suggesting that this divergence is due to selection. Finally, gene flow estimated with microsatellite loci is an average of 30 times lower between populations with different call types than between populations separated by a similar geographical distance with the same call type, demonstrating genetic divergence and incipient speciation. Taken together, these data provide strong evidence that sexual selection is driving behavioural isolation and speciation, supporting sexual selection as a cause for speciation in the wild.
Sexual selection can cause the rapid divergence of sexually dimorphic traits associated with mate acquisition (Darwin 1871; Fisher 1930; Andersson 1994). Many of these traits are involved in mate recognition and result in assortative mating within and among species (Boake 2000; Panhuis et al. 2001; Coyne & Orr 2004). Therefore, it has been proposed that sexual selection can facilitate divergence of mate recognition systems among populations and thus, incidentally, lead to speciation (Fisher 1930; West-Eberhard 1979; Lande 1981; Andersson 1994; Gray & Cade 2000; Boughman 2001; Panhuis et al. 2001; Masta & Maddison 2002). In fact, some researchers suggest that sexual selection is the primary driver of speciation (Carson 2003).
Several studies provide tantalizing evidence that sexual selection may cause speciation in the wild (Coyne & Orr 1989, 1997, 2004; Barraclough et al. 1995; Seehausen et al. 1997; Gray & Cade 2000; Masta & Maddison 2002), yet support for this hypothesis remains incomplete. To demonstrate conclusively that sexual selection drives speciation, several criteria must be met: male signals must differ significantly among populations or recent sister species; females must prefer local male signals to foreign ones; divergence in male signals must result from selection and not genetic drift; and finally, divergent signals and preferences must be correlated with restricted gene flow at nuclear loci to show that populations are diverging towards distinct species (i.e. incipient species).
Several studies have shown that male signals differ among populations or sister species (e.g. Coyne & Orr 1989, 1997; Seehausen et al. 1997; Masta & Maddison 2002), and some have also demonstrated female preferences for local male signals (e.g. Seehausen et al. 1997; Gray & Cade 2000). Divergence in male signals has proceeded faster than expected by genetic drift in Habronattus jumping spiders (Masta & Maddison 2002), but in this case, females do not prefer local male signals (Hebets & Maddison 2005). Hoskin et al. (2005) also recently demonstrated premating isolation among populations of green-eyed tree frogs, Litoria genimaculata, but in this species, premating isolation and speciation were caused by reinforcing natural selection against hybridization, not by sexual selection. Thus, although the above studies are suggestive, no single study has provided a complete body of evidence for the hypothesis that sexual selection can cause speciation. Here, we show that the above criteria are met in the Amazonian frog Physalaemus petersi.
Recent work has revealed an extraordinary divergence in male advertisement calls used in mate recognition among populations of P. petersi (Boul & Ryan 2004). Males from five sites in Ecuador and Peru only produce simple calls with a single component, the ‘whine’ (figure 1a,b; electronic supplementary material, audio files 1–3). In contrast, males from five other sites in Ecuador, Peru and Brazil facultatively add a second component to their calls, a ‘squawk’, producing a complex call. In the túngara frog (Physalaemus pustulosus), the sister species of P. petersi (Ron et al. 2006), the whine is necessary and sufficient for mate recognition, but the chuck, a call component homologous to the P. petersi squawk, further enhances the call's attractiveness, so that females prefer complex to simple calls (Ryan 1980).
We hypothesized that sexual selection for complex calls in P. petersi has driven call divergence between neighbouring populations, which, in turn, is generating reproductive isolation and initiating speciation. We tested this hypothesis using a combination of: (i) call analysis to characterize differences among populations; (ii) phonotaxis experiments to test whether females prefer complex to simple calls; as well as (iii) whether females prefer simple calls of local males to those of foreign males; (iv) coalescent simulations using a mtDNA phylogeny to test whether among-population call variation was caused by divergent selection or genetic drift; and (v) analysis of microsatellite loci to test whether gene flow is restricted between populations.
2. Material and methods
(a) Call analysis
Advertisement calls of P. petersi males were recorded with a Sennheiser SE66 microphone (frequency response 40–20 000 Hz), a Sony Walkman WM-D6C professional stereo cassette recorder (frequency response 40–15 000 Hz) and metal cassette tapes. Calls were digitized using Signal (Engineering Design, Belmont, MA) at a sampling rate of 25 kHz. Digitized calls, one from each male, were then analysed using batch processing in Signal. Batch processing enforces a degree of standardization that is sometimes lost when calls are analysed individually. The dominant frequency of the entire whine was measured from a fast Fourier transform. The transform length was 8192 points, which, given the sampling rate of 25 kHz, provided a frequency resolution of 3 Hz. The significance of differences in the dominant frequency of whines between populations with simple and complex calls was tested using Kruskal–Wallis tests.
(b) Phonotaxis experiments
We used standard phonotaxis experiments (Ryan 1980) to test female preferences for squawks and local whines. In phonotaxis experiments, female preferences are tested by broadcasting different male calls and allowing females to choose between calls, as indicated by contacting one of the speakers. This is an accurate bioassay for mate choices, because P. petersi females only show phonotaxis to choose a mate. We tested females from Yasuní with simple calls from each population and simple versus complex calls from their own population. La Selva females were tested using the same set of simple calls, and using their simple calls versus an artificial complex call produced by digitally appending a squawk from a Yasuní call to a simple La Selva call. The phonotaxis chamber was 92 cm wide, 188 cm long and 92 cm high. Females were placed in the middle of the chamber and calls were broadcast antiphonally from either end. A choice was made when the female came within 11 cm of either speaker. Detailed phonotaxis methods are provided in the electronic supplementary material.
(c) Phylogenetic analysis
We sequenced samples from 42 individual P. petersi and 4 outgroup taxa (electronic supplementary material, table 1). Tissue samples (liver, muscle and/or toes) were stored in 95% ethanol, tissue buffer, DMSO buffer or frozen. Physalaemus pustulosus is the sister species of P. petersi, which together form a well-supported clade (Edentulus) that is sister to another clade (Duovox) containing all other species in the P. pustulosus species group (Ron et al. 2006). Therefore, for outgroup taxa, we used two P. pustulosus specimens and one each of two species in the Duovox clade, Physalaemus pustulatus and Physalaemus coloradorum, for a total of four individuals. Southern populations of P. petersi in western Brazil, southeastern Peru and Bolivia have been recognized as a different species, Physalaemus freibergi (Cannatella et al. 1998). However, since the geographical boundaries and diagnostic characters of P. freibergi are poorly understood, we refer to all the western Amazonian populations of Physalaemus as P. petersi in the present analysis.
Although Nascimento et al. (2005) resurrected the genus Engystomops for the P. pustulosus group, that action is not consistent with their own analysis of relationships. Ron et al. (2006) uncritically followed the use of Engystomops. However, one of the authors of the latter paper (David C. Cannatella) agrees that the resurrection of Engystomops as a genus was unjustified; therefore, the use of Physalaemus is continued here.
Total genomic DNA was extracted from tissue samples using DNeasy Tissue Kits (Qiagen, Inc., Valencia, CA). Overlapping sets of primers were used to amplify approximately 2.4 kb of the mtDNA 12S and 16S genes and the intervening valine tRNA using the polymerase chain reaction (PCR). We used the primers, PCR conditions and sequencing protocol described by Pauly et al. (2004). Editing and assembly of contigs were completed using Sequencher v. 4.5 (Gene Codes Corp., Ann Arbor, MI).
Initial alignment of DNA sequences was completed in ClustalX (Thompson et al. 1997). Manual adjustments were then made in MacClade v. 4.06 (Maddison & Maddison 2000) to minimize the number of changes required across sites. In the final alignment, there were 17 nucleotide positions in which the majority of individuals had an apparent deletion. Exclusion of these 17 sites from the matrix did not affect the topology of the maximum parsimony tree. Parsimony, likelihood and Bayesian analyses all yielded the same tree topology, so only the likelihood tree, with Bayesian posterior probabilities, is shown (figure 1c). Detailed phylogenetic methods are provided in the electronic supplementary material.
(d) Coalescent analysis
In order to test whether the observed pattern of call variation was caused by divergent selection or genetic drift, we used an indirect method developed by Masta & Maddison (2002) for the five populations from Ecuador with larger sample sizes of sequences (six per population) and calls (mean of 14 males per population). This method assumes that nuclear genes code for the phenotypic character of interest (in our case, complex versus simple calls) and asks whether the degree of fixation of phenotypes observed is likely to occur under the assumption of neutrality. This method involved three steps. First, s of Slatkin & Maddison (1989), a measure of incomplete lineage sorting among populations, was calculated for maximum parsimony (MP) trees inferred from the analysis of mtDNA sequences. Larger s values indicate greater levels of incomplete lineage sorting (suggesting a more recent population divergence) and smaller s values indicate lower levels of incomplete lineage sorting (suggesting an older population divergence). For the two MP mtDNA trees, the value of s was equal to 8.
Second, coalescent simulations (using Mesquite v. 1.06; Maddison & Maddison 2003) were used to estimate the upper 95% confidence limit (CL) for the number of generations since population divergence (scaled by effective population size, Ne) that would be expected to give the observed s value for the mtDNA trees. We used 10 000 simulations and Ne=500 (although the choice of Ne does not affect results, as all calculations are scaled by Ne). Using the upper 95% confidence limit for time since divergence rather than the mean is conservative, because a greater time since divergence increases the probability of fixation of phenotypic traits under neutrality. In our analysis, the upper 95% CL for time since population divergence was 1.15Ne generations under the assumption of simultaneous population divergence and 1.85Ne generations under the assumption of bifurcating divergence.
Finally, estimates of generations since population divergence were used in coalescent simulations of nuclear genes to estimate the probability of complete fixation of call types (s=4 for five populations) by genetic drift, assuming that this trait is encoded by nuclear genes. These simulations used Ne=2000, because the effective population size of nuclear genes is four times the effective population size of mitochondrial genes. Details on the assumptions of this method are provided in the electronic supplementary material.
(e) Gene flow estimation at microsatellite loci
We analysed genetic variation at eight variable microsatellite loci for 30 frogs from each of three populations (total 90) separated by similar geographical distances (21–28 km), but with different call types and female preferences: La Selva with simple calls; and Yasuní and Tiputini with complex calls (figure 1b). Microsatellite primers for P. petersi were developed by Genetic Identification Services (Chatsworth, CA), and primer sequences, GenBank accession numbers and PCR annealing temperatures are shown in table 2 (see electronic supplementary material). DNA was extracted as described in §2c, and PCR and fragment analyses were performed as described by Lampert et al. (2003). Fragment data were scored using GeneMarker v. 1.3 (Soft Genetics, LLC). Two loci could not be amplified in frogs from La Selva. Therefore, we used six loci in all three populations for subsequent analyses.
Standard genetic analyses, including tests for Hardy–Weinberg (HW) proportions, linkage disequilibrium and genic differentiation, and FST calculations, were performed in Genepop (Raymond & Rousset 1995). Genotypic proportions differed significantly from HW expectations in 3 out of 15 possible tests (p<0.05), but no loci consistently had an excess of homozygotes in all populations, as would be expected if there were null alleles. Moreover, only 1 out of 31 possible tests for linkage disequilibrium was significant, less than that expected by chance, indicating that loci were independent.
Gene flow (Nm) was estimated using a maximum-likelihood coalescent-based approach implemented in Migrate (Beerli & Felsenstein 2001). For each run of Migrate, we used 10 short chains of 50 000 sampled and 500 recorded trees, followed by three long chains of 500 000 sampled and 5000 recorded trees. Four-chain heating was used with temperatures set to 1, 1.5, 3 and 6 to improve the sampling of tree space. Initial estimates of theta and 4Nm were generated from FST values. We conducted this analysis four times using different random numbers to ensure that runs converged on similar relative parameter estimates. Migrate and FST gave similar estimates of Nm, hence only estimates from Migrate are shown in the text. Although Migrate estimates Nm in both directions (e.g. population 1 to 2, as well as population 2 to 1), we were only interested in the average Nm between populations; therefore, we show only average Nm estimates between each pair of populations in §3.
(a) Call variation
Analysis of the calls of 101 males from 10 P. petersi populations confirmed that complex calls (whine plus squawk) are restricted to five populations (figure 1b). In addition, whines from populations with squawks had significantly lower dominant frequencies (, n=40) than from those without squawks (644±23 Hz, Kruskal–Wallis test, p<0.001, n=61). The same pattern in whine frequencies held for the two populations in Ecuador that were used for female choice experiments: Yasuní (with squawks) and La Selva (without squawks; figure 1a,b; electronic supplementary material, audio files 1 and 2). Dominant frequencies of whines from Yasuní (436±8 Hz, n=15) were lower than those from La Selva (742±37 Hz, p<0.001, n=24).
(b) Female choice experiments
To determine whether call differences between Yasuní and La Selva influence mating preferences, we conducted phonotaxis experiments. In the first set of tests, females from Yasuní (where males squawk) strongly preferred the simple call plus a squawk to the simple call alone (exact binomial probability, p=0.007, n=15; figure 2a). Thus, there is sexual selection for complex calls in Yasuní. Females from La Selva (where males do not squawk) showed no preference for complex calls over simple calls (p=0.481, n=18), although the preference for squawks was not significantly different between Yasuní and La Selva (Fisher's exact test, p=0.134, n=33). In the second set of experiments, females from both the populations strongly preferred their local simple call to the simple call of foreign males; for the Yasuní females, the preference was unanimous (p=0.001, n=15) and for the La Selva females, it was nearly unanimous (p=0.001, n=19; figure 2c). These results show: (i) the existence of sexual selection for squawks in the population that has squawks, and (ii) strong behavioural isolation between these populations based on the differences in simple calls, despite being separated by only 21 km.
(c) Phylogenetic and coalescent analysis of mtDNA
Phylogenetic analysis of mtDNA sequences indicates that populations throughout the western Amazon basin have become fixed for simple or complex calls multiple times (figure 1c). Coalescent simulations (Masta & Maddison 2002) of call and mtDNA sequence data from five populations with large sample sizes of sequences and calls were used to differentiate the roles of selection and genetic drift in call evolution. If simple and complex call types, assumed to be encoded by nuclear genes, have sorted significantly faster than mtDNA haplotypes, then the hypothesis of divergent selection rather than drift is supported. We found that mtDNA haplotypes, unlike call types, were not completely sorted by population, suggesting relatively recent divergence of populations and implying rapid fixation of call types (figure 1c). The simulations gave very low probabilities (p<0.0001) that nuclear genes would be fixed in different populations, as observed for call types, under neutrality (electronic supplementary material, figure 3). These results suggest that the observed pattern of among-population call variation is better explained by divergent selection.
(d) Gene flow estimation at microsatellite loci
To test whether gene flow is restricted between populations with different call types and female preferences, we estimated gene flow using microsatellite loci among three focal populations separated by similar geographical distances (21–28 km), but with different call types and female preferences: La Selva, at which males only produce simple calls; and Yasuní and Tiputini, at which males produce complex calls (figure 1b). Average gene flow (Nm) was Nm=0.04 (95% CI=0.02–0.06) between La Selva and Yasuní (separated by 21 km); 0.04 (0.02–0.06) between La Selva and Tiputini (separated by 28 km); and 1.19 (1.07–1.32) between Yasuní and Tiputini (separated by 26 km). Thus, gene flow was an average of 30 times lower between populations with different call types than between populations with the same call type, suggesting that populations with different call types and female preferences are incipient species. Moreover, two out of eight microsatellite loci that were successfully amplified with PCR in the Yasuní and Tiputini populations could not be amplified in the La Selva population (even after repeated attempts), indicating high genetic divergence.
Three out of six microsatellite loci (A11, A114, D1r) had non-overlapping allele sizes between La Selva versus Yasuní and Tiputini and were therefore diagnostic of these different populations. Neither La Selva×Yasuní nor La Selva×Tiputini hybrids were detected with these loci.
The above call data, behavioural experiments and genetic analyses provide strong support for the hypothesis that sexual selection has driven the observed pattern of variation in the presence/absence of complex calls in P. petersi, incidentally generating behavioural reproductive isolation among populations and initiating speciation. The first set of female phonotaxis experiments shows selection for squawks at Yasuní, where males produce complex calls, but no preference for squawks at La Selva, where males only produce simple calls. Thus, sexual selection appears responsible for the evolution of squawks. However, preferences for squawks were not significantly different between Yasuní and La Selva with the sample size available, thus additional factors besides differences in squawk preferences could contribute to the lack of squawks at La Selva. The second set of female phonotaxis experiments shows that populations with different call types are strongly behaviourally isolated based on female preferences for the local whine. Furthermore, the coalescent analysis suggests that fixation of call types, assumed to be encoded by nuclear genes, has occurred faster than expected by genetic drift, providing further evidence that sexual selection is driving call evolution. Finally, gene flow is strongly restricted between populations with different call types, suggesting incipient speciation.
Although the behavioural isolation between populations is based on preferences for the whine, and not for squawks, we argue that this is a secondary result of female preference for squawks. We suggest that the evolution of squawks in the Yasuní frogs resulted in lower frequency whines, which, in turn, drove the strong divergence in local call preferences between La Selva and Yasuní (figure 2). Recently, it was shown that ablation of the fibrous mass (a structure associated with the vocal folds) of túngara frogs precludes the production of chucks, supporting a functional relationship between complex calls and morphology (Gridi-Papp et al. 2006). Males from Yasuní, and other populations with complex calls, have larger body sizes and proportionally larger larynges and fibrous masses than males from La Selva and other populations with simple calls (figure 2b; Boul & Ryan 2004). This suggests that selection for squawks has caused an increase in the size of these morphological features in populations with complex calls. The same correlations are exhibited among species in the genus Physalaemus (Drewry et al. 1982) and within species of the P. pustulosus species group (Ryan & Drewes 1990).
The difference in whine dominant frequency is the most probable candidate for generating the local mate preferences between La Selva and Yasuní, as preference based on spectral aspects of the whine is well known in P. pustulosus (Wilczynski et al. 1995; Bosch et al. 2000; Ryan et al. 2003). Other alternatives are that sexual selection for lower frequency whines has resulted in larger larynges, which, in turn, can produce squawks, or that natural selection favours larger body size in some populations. Both of the alternatives seem less probable because fibrous masses, which are necessary for producing squawks (Gridi-Papp et al. 2006), are proportionally larger in populations with squawks after correcting for larynx size (Boul & Ryan 2004). This suggests that selection for larger fibrous masses is primary and the correlated increases in larynx and body size are secondary. However, under all alternatives, our argument that sexual selection is driving behavioural isolation and speciation in P. petersi is supported by the data and the analyses.
Another potential cause of divergence of calls and preferences in P. petersi is reinforcing selection against maladaptive hybridization (Dobzhansky 1937). But hybridization is a necessary prerequisite for the reinforcement hypothesis (Howard 1993), and it was not detected using diagnostic microsatellite markers between La Selva (with simple calls) and Yasuní or Tiputini (with complex calls). Additionally, it is unlikely that frogs from these different populations frequently come into contact owing to low gene flow among P. petersi populations shown in this and another study (Gascon et al. 1998), indicating historically low potential for hybridization. Thus, although we cannot exclude reinforcement, the data suggest that sexual selection is the primary driver of the observed divergence in reproductive characters.
Although the Napo River runs between La Selva (on the north side) and Yasuní and Tiputini (on the south side; figure 1b), two lines of evidence suggest that it does not strongly restrict gene flow in P. petersi. First, if the Napo River were a genetic barrier for P. petersi, then populations on the north side (Cando and La Selva) and south side (Jatun Sacha, Yasuní and Tiputini) of the river should form reciprocally monophyletic clades, but they do not (figure 1b,c). Similarly, if the Napo were responsible for call divergence in P. petersi, then populations on opposite sides of the river should consistently have different call types (e.g. all northern populations with simple calls and all southern populations with complex calls), but this is also not the case (figure 1b). The Napo River also fails to act as a barrier for other frogs (Symula et al. 2003). Second, gene flow was not restricted among populations of P. petersi by another major Amazonian tributary, the Juruá River, in western Brazil (Gascon et al. 1998). The Juruá also fails to act as a barrier for many other taxa (Gascon et al. 2000). The general lack of an effect of rivers in the western Amazon basin on gene flow may be due to lateral channel migration, which causes across-river transfer of large pieces of land (Räsänen et al. 1987). Thus, restricted gene flow in P. petersi between La Selva and Yasuní and Tiputini is probably caused primarily by behavioural isolation, not by the Napo River. Nevertheless, our female choice experiments and coalescent analyses show that sexual selection, not genetic drift, is driving behavioural isolation.
The pattern of geographical variation of call types in P. petersi is in stark contrast to the situation in its sister species, P. pustulosus, which has a Middle American and northern South American distribution. An examination of 30 populations of P. pustulosus along a 5000 km transect throughout its range showed that all males in all populations produce complex calls (Ryan et al. 1996). Thus, the presence of complex calls is variable among populations of P. petersi, but fixed within P. pustulosus. The relationship between time since divergence and mate preferences also differs substantially between P. pustulosus and P. petersi. Average sequence divergence is more than twice as high between the two major P. pustulosus clades (3.7%) than between the Yasuní and La Selva populations of P. petersi (1.5%; Ron et al. 2006). Nevertheless, there is no evidence for local mate preferences in P. pustulosus similar in magnitude to those uncovered here between La Selva and Yasuní (Pröhl et al. 2006).
This study represents one of the most convincing examples of speciation by sexual selection for any species in nature. We have demonstrated: (i) divergence in male calls, (ii) strong female preferences for local calls, (iii) that fixation of call types is not caused by genetic drift, and (iv) severe reduction in gene flow between populations with different call types, showing incipient speciation. We hypothesize that call divergence and behavioural isolation are incidental effects of sexual selection for complex calls.
We thank A. Angulo, F. Ayala, R. Boul, J. Guayasamin, M. Guerra, K. Koe, D. Lombeida, K. Mott, K. Newcomb, S. Padilla, S. Ron and I. Tapia for their field assistance; the community of Cando, Jatun Sacha Biological Station, E. Schwartz (La Selva Lodge), D. Romo and K. Swing (Tiputini Biodiversity Station), the Estación Científica Yasuní, and Kurt Holle (Tambopata Research Centre) for field accommodations; U. Mueller for access to the automated sequencer; S. Boles, B. Caudle and C. Peden for their lab assistance; L. Coloma (Pontificia Universidad Católica del Ecuador) for tissue loans; X. Bernal and B. Dawson for their assistance with call analysis; M. Mahoney, S. Masta, G. Pauly and D. Zwickl for their assistance with phylogenetic and coalescent analyses; S. Arnold, S. Barrett, K. Hoke, M. Kirkpatrick, P. Martin, S. Ron and two anonymous reviewers for providing their helpful comments on our manuscript; and the Ecuadorian Ministerio de Ambiente (nos. 004-IC-FAU-DFP and 004-IC-FAU-DNBAPVS/MA) and the Peruvian Instituto Nacional de Recursos Naturales (no. 013-2004/INRENA-IANP-J-PNBS/RNTAMB) for research permits and P. Barriga (Fundación Yasuní), J. Espinoza (INRENA) and J. Córdova (Museo de Historia Natural, Universidad Nacional Mayor de San Marcos) for their assistance with permits. Sequences have been deposited at GenBank (electronic supplementary material, tables 1 and 2). This research was supported by National Science Foundation grants IBN 98-16564, IBN 99-81631 and IRCEB 00-78150. Animal care procedures were approved by the University of Texas (no. 0302170). This is publication no. 58 of the Yanayacu Natural History Research Group.