Advertisement calls play an important role in anuran reproduction, carrying information about mate location, mate quality, and species identity (Wells 2007). Because advertisement calls are often species-specific signals, call differences between populations or species might prevent recognition of potential mates and, therefore, result in reproductive isolation (Gerhardt 1974; Oldham and Gerhardt 1975; Telford and Passmore 1981; Backwell and Jennions 1993; Ryan and Rand 1993; Lemmon 2009). Call differences might also be relevant to speciation, as there are several examples where speciation in frogs appears to have been driven by call divergence (Platz 1989; Boul et al. 2007; Guerra and Ron 2008; Funk et al. 2009). By examining call variation among populations of a single species, we might understand how potential reproductive barriers evolve during early stages of population divergence and speciation. Although intraspecific call variation in frogs is relatively common (Wilczynski and Ryan 1999), it is often unclear whether or not such variation produces reproductive isolation.

Ranitomeya imitator is a species of dart-poison frog that has undergone a Müllerian mimetic radiation; that is, different populations have evolved distinct color-patterns that mimic several different congeneric species (Symula et al. 2001; Yeager 2009; Twomey et al. 2013). There are four mimetic morphs of R. imitator, all of which occur in different geographic regions in north-central Peru. Where different morphs come into contact, “transition zones” occur that are characterized by the presence of phenotypic intermediates (Twomey 2014; Twomey et al. 2014). This system is ideal for studying how divergence in an ecologically relevant trait (mimetic color-pattern) can drive speciation in that there are three different mimetic transition zones; thus, we can study each transition zone separately for the presence of reproductive isolating barriers. Recently, Twomey et al. (2014) found that mimetic divergence might be driving incipient speciation between two mimetic morphs (striped and varadero), and furthermore, that these morphs have distinct calls despite occurring <4 km apart. However, whether or not these call differences contribute to reproductive isolation between the two morphs remains unknown. Furthermore, call variation across the other two mimetic transition zones has not yet been investigated.

Although mimetic divergence itself should not drive call differentiation (most divergence should occur only in traits related to mimicry), there are a few processes by which call differentiation between mimetic morphs might arise. First, if divergent selection is sufficient to reduce gene flow between morphs, divergence in traits unrelated to mimicry might occur via genetic drift (Feder and Nosil 2010). In these situations, divergence in most traits occurs as a consequence of reduced gene flow between morphs. A similar pattern could arise via divergence in allopatry followed by secondary contact. Second, call divergence could occur if traits subject to divergent ecological selection have pleiotropic effects on calls. For example, selection for different body sizes can, in turn, influence call frequencies and pulse rates (Duellman and Trueb 1986; Wilczynski and Ryan 1999). Third, because aposematic traits are subject to positive frequency-dependent selection (Endler and Greenwood 1988; Greenwood et al. 1989), different morphs should become fixed in different geographic areas, creating parapatry between morphs. Thus, because different morphs often occupy different geographic areas, they might also occupy distinct habitats, which could drive call divergence directly (e.g., if different calls were better suited to different local conditions) or indirectly (e.g., if different habitats caused divergence in body size, indirectly driving call divergence). Finally, it is possible that reproductive character displacement could cause call divergence in areas where different lineages come into contact (Höbel et al. 2003; Lemmon 2009). Because cross-morph matings can produce nonmimetic offspring, individuals that mate with their own mimetic morph could have higher fitness than individuals that mate at random (Noor 1999). Therefore, call differences between morphs could arise as an adaptation to avoid cross-morph matings. A key prediction of the last hypothesis is that call differences between morphs should be greater between sympatric (or parapatric) populations versus allopatric populations (Lemmon 2009).

In this study, we investigate call variation in R. imitator to determine whether or not any differences in call parameters exist and, if so, to what degree they coincide with changes in mimetic color-pattern. We also investigate the influence of elevation and body size on call variation in R. imitator in order to understand potential mechanisms explaining call variation across the mimetic transition zones.

Materials and Methods

We recorded 75 individuals of R. imitator from 13 localities in the departments of San Martin and Loreto, Peru. We included recordings from an additional 58 R. imitator from 7 localities published in Twomey et al. (2014) corresponding to the striped-varadero transect. Our sampling localities span the three mimetic transition zones in R. imitator (banded-striped, 9 localities, 61 recordings; spotted-striped, 6 localities, 36 recordings; striped-varadero, 7 localities, 58 recordings; Table 1), allowing us to investigate whether aspects of the advertisement call vary across mimetic transition zones. Transect position was calculated as the distance from a fixed point representing the estimated transition zone center. See Twomey et al. (2014) and Twomey (2014) for details on transect position calculations.

Table 1.

Sampling localities, transect positions, and call sample sizes of the transition zones for each of three mimetic transition zones of Ranitomeya imitator. For all coordinates, GPS datum  =  WGS84.

Call recordings were obtained by searching for calling males in the field. Ranitomeya imitator has two call types: a longer advertisement call (Fig. 1), and a shorter call used for close-range courtship (Mayer et al. 2014). The latter call is used after courtship has already been initiated, and might have a lesser importance in attracting a potential mate. Therefore, we focused on the longer advertisement calls that are recognized by having longer note lengths and internote intervals than courtship calls. Calls were recorded with a Marantz PMD660 digital recorder with a Sennheiser ME66-K6 microphone. We extracted three bioacoustic variables commonly used in dendrobatid studies (Brown et al.2011): dominant frequency (frequency at which peak amplitude is reached; Hz), note length (length of a note measured from the start of the first pulse to the end of the last pulse, where a note is described as a discrete bundle of pulses; s), and pulse rate (pulses/s, calculated by dividing the number of pulses contained within a note by the note length). Call measurements were taken in Raven Pro v1.3 (Charif et al. 2008). For each male, a recording generally consisted of several notes. Measurements were always taken on ≥3 notes and averaged for each male. To account for any influence of temperature on calls, we took temperature measurements alongside each call recording in the microhabitat of the calling male. We then standardized each of the three call variables (dominant frequency, note length, and pulse rate) by calculating regression residuals. Any additional statistical analyses were performed on these temperature-standardized variables rather than on raw measurements. Residuals were calculated in SPSS (v20; SPSS Inc., Chicago, IL) on the entire data set simultaneously (i.e., not separated per transect). Comparisons of call parameters between the varadero and striped-transition populations were performed with a one-way analysis of variance in SPSS. This comparison was made because these two populations exhibit strong call differences across a short geographic distance.

Fig. 1.

Advertisement calls of Ranitomeya imitator, showing waveforms (above) and spectrograms (below). (a) varadero morph, recording locality Varadero Forest 1, recorded at 27°C (note length  =  0.983 s, pulse rate  =  37.6 pulses/s, dominant frequency  =  5059 Hz); (b) striped morph, recording locality Varadero South Bank, recorded at 26°C (note length  =  0.692 s, pulse rate  =  34.7 pulses/s, dominant frequency  =  5668 Hz); (c) banded morph, recording locality Sauce, recorded at 26°C (note length  =  1.044 s, pulse rate  =  28.7 pulses/s, dominant frequency  =  5003 Hz); (d) spotted morph, recording locality San Jose, recorded at 22°C (note length  =  0.777 s, pulse rate  =  33.5 pulses/s, dominant frequency  =  5194 Hz). The nearly ubiquitous noise occurring at 7 kHz in all recordings is background noise caused by calling insects.

To describe spatial variation (i.e., variation along each transect) in call parameters, we fit three candidate models—a flat, linear, and sigmoid model—to the data and evaluated each using an information theoretic approach (ΔAIC_c). A flat model describes a scenario where a call parameter is constant across the transect and consists of a single parameter, the population mean (defined as the grand mean of all individuals). A linear model describes a constant, smooth change across the transect, and has two parameters—slope and y-intercept. The linear model was fit with linear regression. Finally, a sigmoid model describes a scenario where a call parameter changes abruptly along the transect. This is expected, for example, if two mimetic morphs sharing a transition zone have distinct calls. It is possible that an abrupt change could occur anywhere along the transect, however, not just at a mimetic transition zone. Therefore, we constrained the center parameter to equal that of the center of the color-pattern clines found in previous studies (Twomey 2014; Twomey et al. 2014). For the banded-striped transect, the center was constrained to 0.99 km; for the spotted-striped transect, the center was constrained to 1.44 km; and for the striped-varadero transect, the center was constrained to 0.52 km. By constraining call-centers to color-pattern centers, we tested the hypothesis that there is an abrupt change in call parameters that occurs at the same geographical location as the change in color-pattern. For the sigmoidal model, we performed nonlinear regression using a 4-parameter tanh function

where c is the center of the cline, w is the width of the cline, and y_max and y_min are the maximum and minimum trait values. Parameter searches were done using the solver function in EXCEL using a least-squares optimality criterion. Solver was run using the GRG nonlinear algorithm with the following settings: convergence  =  0.0001; central derivatives; multistart on; population size  =  100. To evaluate which of the three models was a better fit to the data (1-, 2-, or 4-parameter), we calculated ΔAIC_c and Akaike weights (w_i) for each model using the residual sum of squares divided by the sample size as the likelihood criterion (Burnham and Anderson 2002).

To describe relationships between call parameters and male body mass (mass data from (Twomey 2014; Twomey et al. 2014), we regressed the population average of the mass of male frogs with the population average of each call parameter. We used population-level data (i.e., averaging across individuals) because we did not have corresponding mass data for each call recording. We also regressed average male mass on elevation to test for a relationship between mass and elevation. Linear regressions were calculated in SPSS.

Results

For the banded-striped transect, note length showed no consistent pattern of variation along the transect, whereas a gradual, linear increase in dominant frequency occurred from the banded to striped morph (Fig. 2a). Pulse-rate variation was best explained by a sigmoidal model, even with the center of the cline constrained a priori to 0.99 km. Thus, there was evidence that a stepped increase in pulse rate was geographically coincident with the color-pattern cline, although the width of the pulse rate cline was quite wide (11.63 km; Table 2).

Fig. 2.

Clines in advertisement call parameters across three mimetic transition zones of Ranitomeya imitator. (a) Banded-striped transition zone, (b) spotted-striped transition zone, (c) striped-varadero transition zone. In all panels, call parameter values for individual R. imitator (represented by dots) are plotted along the sampling transect (x-axis). For all panels, the fit line represents the best-supported model according to the AIC_c (see Table 2).

Table 2.

Model fit results for call clines of three mimetic transition zones of Ranitomeya imitator. In all cases, the best supported model is shown in bold. For cases where the sigmoid model was best supported, point estimates for cline parameters are given.

For the spotted-striped transition, we found that neither pulse rate nor note length showed a consistent pattern of variation along the transect (Fig. 2b). We observed a weak, linear increase in dominant frequency from the spotted to striped morph, although the relationship between transect position and dominant frequency was marginal (R²  =  0.09; Fig. 2b).

For the striped-varadero morph, variation in both note length and dominant frequency was best explained by a sigmoidal model, even when the center of the cline was constrained a priori to 0.52 km (Fig. 2c). Thus, along the transect in the direction of the varadero morph, there was evidence of a stepped increase in note length and decrease in dominant frequency being geographically coincident with the color-pattern cline. The call-parameter cline widths were relatively narrow: 4.89 km for note length and 3.52 km for dominant frequency (Table 2). Variation in pulse rate was best described by a linear model (Fig. 2c), showing a slight decrease in pulse rate along the transect (from the striped to varadero morph).

Across all populations, we found that mean male mass was not correlated to elevation (R²  =  0.11, P  =  0.22; Fig. 3a). This lack of a mass–elevation relationship is exemplified by lowland populations (elevation <300 m), showing a range of mass variation exceeding variation shown across all other populations. Neither note length nor pulse rate were correlated to average male mass across all populations (note length, R²  =  0.04, P  =  0.44; pulse rate, R²  =  0.03, P  =  0.54; Fig. 3b,c). However, dominant frequency was negatively correlated to average male mass across all populations (R²  =  0.31, P  =  0.03; Fig. 3d).

Fig. 3.

Regressions between (a) elevation and average male mass, (b) average male mass and note length residuals, (c) average male mass and pulse rate residuals, and (d) average male mass and dominant frequency residuals of Ranitomeya imitator. Individual dots represent population averages.

When comparing call parameters between the varadero and striped-transition populations, we found that these populations differed with respect to note length (F_1,38  =  8.91, P  =  0.005) and dominant frequency (F_1,37  =  17.55, P <0.001), but not pulse rate (F_1,35  =  0.26, P  =  0.62; Fig. 4).

Fig. 4.

Boxplots comparing call parameters between the striped morph (open boxes) and varadero morph (shaded boxes) at the striped-varadero transition zone. Asterisks indicate differences in call parameters between mimetic morphs of Ranitomeya imitator.

Discussion

We documented intraspecific variation in a number of call traits among populations of R. imitator. Dominant frequency varied across all three transects, although variation in pulse rate and note length also occurred (Fig. 2). Our sampling revealed clinal variation in call traits across two mimetic transition zones. Clinal variation in advertisement calls has been found in a number of frogs (Narins and Smith 1986; Ryan et al. 1996; Hasegawa et al. 1999; Bernal et al. 2005) and often corresponds with an underlying environmental gradient such as elevation or latitude. In R. imitator, an underlying elevation gradient might explain some of the variation in dominant frequency across the banded-striped and spotted-striped transition zones. In these transition zones, the lowland striped morph generally has a higher dominant frequency than either the banded or spotted morphs, which are both from highland habitats. Thus, there appears to be a general transition toward lower dominant frequency moving up the elevation gradient, even when controlling for recording temperature. One possible explanation for this pattern might involve a change in body mass (e.g., larger highland frogs producing lower frequency calls). We found no correlation between body size and elevation, indicating that the transition toward a lower frequency in the highland populations is likely not driven by an elevation effect on body size (Fig. 3). In addition, we found that clinal variation in advertisement calls tracks clines in mimetic color-pattern in two of the three transition zones. In the banded-striped transition zone, we found a stepped change in pulse rate corresponding to the cline in mimetic color-pattern. The change occurs over a relatively wide area (11.63 km) but involves an increased pulse rate in the striped morph. Thus, these two mimetic morphs appear to be somewhat distinct with respect to pulse rate. There is some evidence that the striped and banded morphs assortatively mate, and that color-pattern specifically is one of the cues (Twomey 2014). However, whether call differences also contribute to assortative mating should be investigated.

In the striped-varadero transition zone, we documented a change in both note length and dominant frequency, both coincident with the relatively narrow (<5-km) mimetic transition zone. Thus, at the striped-varadero transition zone, the two mimetic morphs are well-differentiated in two aspects of their advertisement calls (Fig. 4). In addition to changes in call parameters across this transition zone, there is also a strong change in body size, with striped males being much smaller than varadero males (Twomey et al. 2014). Thus, the change in dominant frequency is likely caused by the change in mass, reflecting a well-established negative relationship between body size and call frequency (Duellman and Trueb 1986). This relationship was also confirmed for R. imitator in the current study (Fig. 3d). We found no such relationship between male mass and note length (Fig. 3b), however, which indicates that the change in note length cannot be explained by differences in body size between the two morphs. In a previous study, Twomey et al. (2014) found that the striped and varadero morphs are somewhat genetically isolated, and that morph-based assortative mating occurs at their transition zone. Thus, one possibility, related to the observed call differences, is that this isolation is mediated by differences in advertisement calls.

Advertisement calls in the vanzolinii clade of Ranitomeya (a species group containing R. cyanovittata, R. flavovittata, R. imitator, R. sirensis, R. vanzolinii, and R. yavaricola) are remarkably conserved (Brown et al. 2011): all species have a short, musical trill. For example, the breadth of variation in calls within R. imitator meets or exceeds the variation seen across the entire species group. Across all recorded R. imitator individuals, note length ranged from 0.27 to 1.38 s, pulse rate ranged from 29 to 47 pulses/s, and dominant frequency ranged from 4703 to 5977 Hz. Regarding call parameters for other species in the clade, note length ranges from 0.57 s (in R. vanzolinii) to 2.2 s (R. sirensis), pulse rate from 24 pulses/s (R. sirensis) to 32 pulses/s (R. yavaricola), and dominant frequency from 5010 Hz (R. sirensis) to 6000 Hz (R. yavaricola; Brown et al. 2011). Thus, in many cases, there is less variation in call parameters across the entire clade than observed among individual R. imitator. The exception is R. sirensis, which can produce a longer call (up to 2.2 s in length), and a lower pulse rate (down to 24 pulses/s), than has been observed in R. imitator. Overall, this indicates that speciation in the vanzolinii clade has taken place without substantial interspecific call divergence. Color-pattern, however, varies considerably among many of the species in this clade. Color-pattern divergence is frequently implicated as a driver of population divergence in dendrobatids (Wang and Summers 2010; Rudh et al. 2011; Twomey et al. 2014), but research in this group has generally focused on how color-pattern divergence among populations of a single species can lead to within-species reproductive isolation. Thus far, no studies in dart-poison frogs have attempted to address whether a known radiation of species (i.e., a scenario where speciation is already complete) has been driven by color-pattern divergence. Given that most species in the vanzolinii group show similar advertisement calls, other factors—such as color-pattern divergence or biogeographic history—might have been important in producing this pattern of diversification.