Rapid differentiation of sexual signals in invasive toads: call variation among populations

Advertisement calls tend to differ among populations, based on morphological and environmental factors, or simply geographic distance, in many taxa. Invasive cane toads (Rhinella marina) were introduced to Australia in 1935 and their distribution has expanded at increasing rates over time. Rapid evolution occurred in morphological and behavioural characters that accelerate dispersal, but the effects of rapid expansion on sexual signals have not been examined. We collected advertisement calls from four populations of different ages since invasion, and analysed the geographic differentiation of seven call parameters. Our comparisons indicate that the calls of R. marina differ among Australian populations. The signal variation was not simply clinal with respect to population age, climate, or morphological differentiation. We suggest that selection on signalling among populations has been idiosyncratic and may reflect local female preferences or adaptation to environmental factors that are not clinal such as energy availability.

We followed the multivariate comparisons with univariate comparisons of each call parameter among populations. Most parameters differed significantly among populations ( Table 1). The NEQ population had higher dominant frequency, longer call duration with longer pulse length and longer inter-pulse interval than the others ( Fig. 2; Table 1). The three remaining populations were all relatively similar in these parameters, but with  The first (highly significant) axis accounted for 61.0% and the second axis accounted for 32.1% of the total variation among sites. Call parameters with correlations greater than 0.45 of standardised coefficient for each canonical axis are indicated in the text boxes along each axis. The first axis was correlated strongly with inter-pulse interval length and dominant frequency and the second axis was most strongly correlated with dominant frequency, call duration and pulse length. relatively lower dominant frequencies and shorter inter-pulse intervals than the NEQ population ( Fig. 2; Table 1). The SQ population was near the mean for the remaining parameters, pulse length and call duration, however the CQ and WA populations were significantly differentiated along this axis; CQ toads had relatively greater pulse lengths and shorter call durations, while WA toads diverged in the opposite direction, with shorter pulse lengths and longer call durations ( Fig. 2; Table 1).
Calls vary among individuals as well as among populations, so we examined the among-versus within-population variation in individual call parameters. Dominant frequency and pulse rate varied little among individuals within populations, as a static property, typically less than about 5% 27 (Table 2). Inter-call interval length showed the highest variability among individuals within populations ( Table 2). The variation in dominant frequency, call duration, inter-call interval, and inter-pulse interval, were higher among than within populations (i.e., values of CVa/CVp, were higher) ( Table 2), whereas the variation in proportion of time spent calling and pulse rate were about equal within and among populations (as shown in Table 1). Pulse length showed the largest ratio of overall to within-population variation ( Table 2).
Differences in time since colonisation among sites, or climate gradients among sites could have caused clinal variation in call parameters. Given the locations of the sites, our results were not consistent with the existence of a gradient related to time since colonisation (r 2 = − 0.248, F 1,4 = 0.006, P = 0.944). In addition, none of the distances between sites in terms of eight different climatic variables (mean annual, maximal and minimal temperatures, mean temperature warmest and coldest quarters, mean annual rainfall, temperature seasonality and rainfall seasonality) were significantly correlated with the distances in call parameters between locations; the most highly correlated climate factor was mean annual temperature, explaining 25% of the variation in call parameters between populations, which was not significant (r 2 = − 0.250, F 1,4 < 0.001, P = 0.992).

Discussion
The acoustic signals of invasive cane toads differed significantly among populations colonized at different times. This geographic variation indicates that there has been rapid differentiation of call properties since the initial introduction of toads to Australia, and its pattern indicates that differentiation has not been not simply clinal, driven by environmental factors, or related to the differentiation in other characters that has resulted from spatial sorting.
Most parameters, except for pulse rate and proportion of time spent calling, differed significantly among populations after we removed size and temperature effects (Table 1). Our examination of relative levels of variability in parameters indicated that dominant frequency and pulse rate varied less within populations. Gerhardt 27 suggested that call properties with low levels of within-population variation may be those that are strongly

Table 2. Coefficient of variation (%) in within-individual (CVi), within-population (CVp) and overall mean values (CVa) for seven call parameters from 100 males recorded from four localities.
maintained by female choice. Our comparisons of the ratios between among-and within-population coefficients of variation, however, indicated that dominant frequency had relatively high levels of variation among populations. These results suggest that dominant frequency is strongly selected within populations, in directions that differ among populations. Pulse rate, on the other hand, tended not to differ among populations and also varied little within individuals. This call trait often differs among species, driven by evolutionary differentiation for species recognition 1 . Given the absence of closely related frog species within toad distribution in Australia, or even species that sound like them 28 , there may be little selective pressure for species-level differentiation. The degree of variability of this parameter (a range of 10 to 20 pulses/sec) may be limited by morphology. Morphological constraints and/or very strong and consistent sexual selection could account for the relatively low levels of variation of this trait at every level. Proportion of time spent calling did not quite differ significantly among populations, and was relatively less variable among populations than within populations than were call duration and inter-call interval. Thus, proportion of time spent calling has some degree of flexibility within populations, but this flexibility is similar across populations. The highest ratio of among-population to within-population CVs occurred in pulse length, suggesting that variation in this call trait may be strongly reduced by selection within populations, while diverging among populations. Another invasive anuran, the Coqui frog (Eleutherodactylus coqui) showed differences in call properties among populations at different elevations in its introduced range in Hawaii after about 34.5 generations, but these differences were caused by differences in body size and temperature 29 . Differences among populations in terms of dominant frequency and pulse rate, for instance, can generally be explained by morphological differences (body size or mass) and environmental temperatures, respectively 1 . Our data, however, do not support the hypothesis that cane toad sexual signals varied along simple clinal gradients of temperature or rainfall, as reported in many studies of other anurans [13][14][15] . The call differentiation we observed among old and new populations with similar colonisation times (Fig. 1) does not suggest spatial or temporal clines caused by toad arrival times (or geographic distance). Similarly, pure drift based on geographic or temporal distance does not account for the similarity between SQ and the newer populations. Thus, the call variation we observed does not seem to be driven by broad clinal patterns of environmental characteristics or by genetic drift related to time since colonisation.
Divergence in female choice among sites could have driven the differences in male signals we observed among populations 30 . Alternatively, this divergence may be a response to selection caused by non-clinal habitat heterogeneity in environmental conditions, including factors affecting signal propagation in different habitats. In addition, metabolic or energetic resources allocated to signalling may influence the character of the signal 13 , and the allocation of resources to signalling may vary among populations as a function of resource availability, or population density, or both 31,32 .
In summary, our study revealed significant divergence in calls among populations of toads within their invasive range. Our study suggests that, like morphology and dispersal behaviour, sexual signals have also experienced rapid evolutionary changes. Those changes, however, do not appear to be the product of spatial sorting but may instead reflect adaptation to localised sets of behavioural, ecological or environmental conditions. Local modification of acoustic signals in invasive animals may accelerate population growth and facilitate invasion by improving mating success.

Methods
Ethics statement. All procedures in this study were approved by the Animal Ethics Committee at James Cook University (permit no. A1838). All procedures undertaken were in accordance with approved guidelines.
Recording of toad advertisement calls. Advertisement calls of 32 male cane toads in the Townsville region were recorded in March and April 2014 (NEQ; Fig. 1). We sampled calls in the Brisbane region of south-eastern Queensland (SQ), in central western Queensland near Longreach (CQ), and near Kununurra in Western Australia (WA; Fig. 1). Toad advertisement calls from 12 males in SQ, 28 males in CQ and 28 males in WA were collected in February and March 2014.
Toad calls were recorded using a Marantz PMD 661 compact digital audio recorder (D&M Professional, Itasca, USA), equipped with a NTG3 shotgun microphone (RØDE, Australia). We recorded consecutive advertisement calls from each male in .wav sound format with 96 kHz sample rate and 24 bit-resolution with manual level adjustment. Immediately following each recording, we captured the focal animal and measured its body temperature to an accuracy of 0.1 °C using a digital non-contact infrared thermometer (QM-7221, DIGITECH, Australia), then placed it in a plastic bag for transport to the laboratory. After we finished all recordings at a site, we measured microhabitat conditions, including water and air temperatures (°C) and relative humidity (%, using a whirling hygrometer P2520, G H Zeal Ltd, UK). In the lab, collected toads were euthanized and frozen. We quantified body condition of each toad by measuring body size as snout-urostyle-length (SUL, mm) using plastic callipers (Pittsburgh, Taiwan), and as body mass (g) using an electronic balance (PM400, Mettler, Switzerland). Acoustic analysis. Recordings were analysed using Raven Pro 1.4 (Cornell Lab of Ornithology, Ithaca, USA) and Avisoft SASLab Pro 5.2.02 (Avisoft Bioacoustics, Berlin, Germany) to measure or calculate seven different call properties from each of 286 recorded advertisement calls from NEQ (range 2 to 16 calls per male, mean = 8.9 ± 2.7 SD), 95 calls from SQ (3 to 10 calls/male, 7.9 ± 2.2 SD), 202 calls from CQ (3 to 10 calls/male, 7.2 ± 2.2 SD), and 181 calls from WA (4 to 10 calls/male, 6.5 ± 1.6 SD). First we downsampled the sampling rate of each call 96 kHz to 44.1 kHz using r8brain v1.9, in order to reduce frequency grid spacing from 93.8 Hz to 43.1 Hz in Raven Pro, and thus obtain more detailed values of dominant frequency. We measured the following temporal variables of each recorded call from the waveform and the spectrogram (Fig. 3a): (1) call duration (section); (2) inter-call interval (section); (3) pulse duration (ms); and (4) inter-pulse interval (ms). We calculated (5) proportion of time spent calling by dividing the sum of call duration by the sum of call duration plus all inter-call Scientific RepoRts | 6:28158 | DOI: 10.1038/srep28158 intervals. We also counted the number of pulses per call (Fig. 3b), and calculated (6) pulse rate (number of pulses per sec) by dividing pulse number by call duration. We measured (7) the dominant frequency (Hz) using Raven Pro's spectrogram function (1024 points fast-Fourier transform [FFT], overlap 75%, Hamming's sampling window with a frequency resolution of 56 Hz; Fig. 3c). We obtained average values (mean ± SD) of each of the seven call parameters for each individual male.
Geographic and climatic data. Distance in km between each location were determined using Google maps, and climate data were taken from the Australian Bureau of Meteorology databases (www.bom.gov.au/), using long term climate averages from weather stations close to the field sites. Climate variables used were; mean annual temperature, mean maximal and minimal temperatures, mean temperature warmest and coldest quarters, temperature seasonality (standard deviation of monthly means), mean annual rainfall, and rainfall seasonality (coefficient of variation in monthly means).
Statistical analysis. All statistical analyses were performed using the software package R (version 3.1.2., R core development team 2014). Since body size of calling males and temperature often affect advertisement call traits 1 , we performed multiple regressions using each call variable as the dependent variable and body size (SUL) and body temperature as independent variables, across all male toads from all populations. We found that The dimensions correspond to amplitude (milliunits, vertical) and time (sec, horizontal). This figure includes call duration (sec) and inter-call interval length (sec) in two calls and pulse numbers, pulse length (ms) and inter-pulse interval length (ms) in a single call. (b) Spectrogram of the same single call. The vertical dimension shows frequency in kilohertz (kHz) and the horizontal shows time (sec). Colour of the spectrogram indicates the relative intensity of the sound and frequency. (c) Power spectrum of that call (same setting as above) with amplitude (dB) and frequency (kHz) axes. Dominant frequency (kHz) is determined as the frequency at the highest amplitude. dominant frequency and pulse length were significantly related to body size (N = 100, r 2 = 0.226, F 2,98 = 29.91, P < 0.001 for dominant frequency, N = 100, r 2 = 0.165, F 2,98 = 20.56, P < 0.001 for pulse length; Supplemental  Fig. S1), and that pulse rate and inter-pulse interval were significantly related to temperature (N = 100, r 2 = 0.221, F 2,98 = 29.06, P < 0.001 for pulse rate, N = 100, r 2 = 0.117, F 2,98 = 14.12, P < 0.001 for inter-pulse interval; Supplemental Fig. S1). To remove those effects, and thus remove the possibly confounding effects of body size and thermal environment from the effects of population, analyses to examine population effects were carried out on the residuals of each of those four variables from the applicable regression. Before those analyses, all variables were standardized to z-scores. Those z-scores were used for the following analyses.
To investigate differences of call characters among populations we initially performed a multivariate analysis of variance (MANOVA) including all seven call properties and residuals, with hypothesis testing via Wilks' Lambda. We used the first two canonical discriminant axes defined by that MANOVA to create a biplot illustrating the degree and manner in which our four populations differed. Following the significant MANOVA, we used univariate analyses of variance (ANOVAs) to test for significant differences among populations for each of the call parameters; pairwise Tukey's HSD tests were used to determine which populations differed significantly for each call parameter.
We calculated coefficients of variation (CV = 100% × (SD/mean)) for each call parameter within individuals (CVi), within populations (CVp), and across all populations (CVa) using the individual, population or grand mean and standard deviations 33 . Based on these CV values, we determined the ratio of overall (or among) populations and within-population variation as CVa/CVp.
To examine possible effects of geographic and climate differences on how toad calls differed among populations, we calculated Euclidean distance matrices among populations for geographic distance and for each of the eight climatic variables using the first two axes of principle component analysis (PCA) (the "vegdist" function in the "vegan" R package). We summarised the overall variation in toad call parameters using a similar technique by determining the Euclidean distances between the centroids (the first two axes, accounting for 93.1% of the variability) from the canonical discriminant space. This created three distance matrices (toad calls, climate, geographic distance), including all six between site comparisons. We then performed linear regressions to test for linear relationships between pairwise site distances in call variable space and pairwise distances in the geographic and climate variable spaces.