Geographic variation in the advertisement calls of Hyla eximia and its possible explanations

Song description

Based on information from public databases we surveyed 51 locations where H. eximia has been previously collected. We found H. eximia in nine of these locations. Surveys were carried out during a single rainy season to avoid bias due to inter-seasonal call variation and to ensure the presence of the species in the breeding pools. Hyla eximia was alone (allopatric) in four of the nine locations that were sampled. In two locations, H. eximia was sympatric and formed breeding choruses with H. arenicolor, and in the remaining three with T. smithii (Table 1, Fig. 1). Here we considered “sympatric” the populations in which two species were found to simultaneously occupy the same microhabitat during the breeding season (form breeding choruses); this preliminary classification that will have to be confirmed with repeated sampling.

During the summer of 2011 we recorded advertisement calls from males of the three species. Calls were recorded in natural breeding aggregations with a directional microphone (Sennheiser™ME66) connected to a digital recorder (Marantz™PMD660). The microphone was held 1 m in front of the calling male, and the recording volume was adjusted to avoid saturation (0–6 dB). Each male was recorded for 1–1.5 min to ensure that adequate call series were obtained. Then we measured the body temperature of the frog by holding it from a hind limb and pointing an infrared thermometer (Extech™42529; ±0.05 °C) at its body, thus preventing heat transfer from the observer. Finally, we measured the frog’s snout-vent length with a digital calliper (±0.05 mm) and collected a toe-clip from the distal phalanx of the fourth right front-leg digit for subsequent DNA extraction (see Gonser & Collura, 1996).

We used only recordings containing at least 10 advertisement calls. In order to maximize the opportunity to detect differences in the songs of H. eximia between populations, we measured a large number of variables from our recordings. Using Avisoft-SASLab Pro™ we quantified the following attributes (definition, abbreviation; units): call duration (length of the call, continuous trace in the sonogram, CD; s), inter-call duration (silent period between consecutive calls, IC; s), call period (sum of CD + IC, CP; s), number of pulses (number of short signals that compose the call, NP; count), pulse duration (length from the beginning to the end of the pulse, PD; s), inter-pulse duration (length of the silent period between the end of one pulse and the beginning of the next, ID; s), pulse amplitude (peak intensity of the pulse, PA; V), pulse peak frequency (highest energy within the pulse, PPF; kHz), pulse period (the sum of PD + IP, PP; s), pulse rate (pulses per second measured from the start of the first pulse to the start of the last pulse, PR; Hz), dominant frequency (frequency with the highest energy within the call, DF; kHz) and fundamental frequency (the lowest, or reference frequency of series of frequencies -called harmonics- that are its integer multiples, FF; kHz). Temporal variables were measured directly on the spectrogram, all the pulse variables were extracted with the Pulse Train Analysis option, and the frequency variables were obtained with Power spectrum logarithmic function (Fig. S1). Individual frogs were represented in the analyses by the average value of each call attribute.

To detect the call variation in H. eximia we defined three chorus compositions: (1) allopatric, composed of H. eximia males only, (2) H. eximia in mixed choruses with H. arenicolor, and (3) H. eximia in mixed choruses with T. smithii. We also analyzed the calls of (4) H. arenicolor males in mixed choruses with H. eximia, and (5) T. smithii males in mixed choruses with H. eximia.

Call variation

To reduce the number of variables included in multivariate analyses in an objective way, we calculated for each chorus composition the correlation matrix of the twelve variables and discarded one attribute of any pair that had a correlation coefficient of ≥0.7. Since temperature is known to affect several amphibian call attributes (Blair, 1958; Zweifel, 1959; Zweifel, 1968; Gerhardt & Mudry, 1980; Gayou, 1984; Gerhardt & Huber, 2002) we performed an analysis of covariance (ANCOVA) for each species and variable to assess the potential effect of temperature on call attributes from different chorus compositions. Since ANCOVAs showed a significant effect of temperature on call variables, we used ANCOVA residuals for all subsequent analyses.

Once the number of variables had been reduced and the effect of temperature removed, we compared the attributes of the calls of H. eximia from the three chorus compositions using a multivariate analysis of variance (MANOVA) to determine whether and to what extent call variation was determined by locality and/or chorus composition. The absence of a significant effect of locality on call attributes might suggest that song variation occurs at a larger geographical scale, such as by basin, thus when there was no locality effect we repeated the analysis replacing locality with major river basins (Panuco, Lerma, Magdalena and Balsas). All analyses were performed using R (R Development Core Team, 2011). The above analyses were designed to evaluate whether the vocalizations of H. eximia produced at different chorus compositions were different in some -or in a combination of- attributes. To visualize such differences we performed a canonical discriminant analysis on the calls of H. eximia, and used the resulting discriminant function to plot the distribution of the calls of H. arenicolor and T. smithii in the same canonical space. This allowed the visualization of any significant differences revealed by the MANOVA as displacements away the calls of sympatric species in a canonical multivariate space. Finally we ran an analysis of variance to compare the first two canonical scores of H. eximia amongst basins—a major geographic source of variation. Analyses in this section were conducted with NCSS™ (Kaysville, UT).

Genetic variation

Using mitochondrial and nuclear genes we analyzed phylogenetic patterns to explore the extent to which geographic variation of songs across populations is related to genetic variation. Total genomic DNA of 30 individuals of Hyla eximia from 12 localities was extracted from EtOH-preserved toe clippings using the DNeasy Tissue Kit (Qiagen™Inc., Valencia, CA, USA). Standard PCRs were carried out in 25 µl reaction mixes with final concentrations of 1 µM of each primer, 1.5 mM of MgCl2, 0.2 mM of each dNTP, 1xNH4 reaction buffer (50 mM Tris–HCl pH 8.8; 16 mM [NH4]2 SO4), 1.25 units of Taq DNA polymerase, and 1–4 µl of DNA. Successful amplifications were performed using the following protocol: initial denaturing for 5 min at 95 °C, denaturing for 45 s at 94 °C, annealing for 45 s at 48 °C (for Cyt b) and 57 °C (for ATP), extension for 1 min at 72 °C, final extension min for 7 min at 72 °C; denaturing, annealing and the first extension stage were cycled 35–40 times.

We amplified two mitochondrial gene regions; cytochrome-b (754 bp) with the primers MVZ16-H AAATAGGAARTATCAYTCTGGTTTRAT and MVZ15-L GAACTAATGGCCCACACWWTACGNAA (Moritz et al., 1992; Stöck et al., 2008; Stöck et al., 2011), and ATPase subunits 8 and 6 (778 bp) with the primers LysAF CAACCACCCTTGATGAATGCC and C3FR GGGCTGGGGTTKACTATGTG (Stöck et al., 2008; Bryson et al., 2010; Table S2). These mitochondrial gene regions have been found to be informative at different levels of variation within hylid frogs (references given above). We also amplified two nuclear gene regions using the same protocol as above, with a slight modification (see below). The Pro-opiomelanocortin gene region (POMC; 447 pb) with the primers POMC_DRV_F1 ATATGTCATGASCCAYTTYCGCTGGAA and POMC_DRV_R1 GGCRTTYTTGAAWAGAGTCATTAGWGG (annealing temperature = 58 °C; modified from Wiens et al., 2005; Bryson et al., 2010), and the Rhodopsin (Rho) gen region (291 pb) with the primers Rho1A ACCATGAACGGAACAGAAGGYCC and Rho1C CCAAGGGTAGCG AAGAARCCTTC (Bossuyt & Milinkovitch, 2000; Faivovich et al., 2005; Klymus et al., 2010; Table S2).

PCR products were sent for sequencing to the University of Washington High-Throughput Genomics Unit, Department of Genome Sciences, to the Service of DNA Sanger Sequencing. Sequences were aligned with Clustal X (Jeanmougin et al., 1998). Some regions were difficult to align, and were therefore adjusted manually using BioEdit (Hall, 1999). DNAsp (Librado & Rozas, 2009) was used to determine the number of unique haplotypes and invariable and polymorphic sites.

Phylogenetic analyses were performed for each gene using the maximum likelihood algorithm (ML; Felsenstein, 1981). We used the Hasegawa–Kishino–Yano model (HKY ) − G + I (gamma distribution with invariable sites) for mitochondrial genes and the GTR + I + G for nuclear genes. In both analyses the nearest neighbor interchange (NNI) was used. The best model was identified with the program Modeltest 3.7 and MrModeltest modified version of Modeltest version 3.6 (Posada, 2008) using the Akaike information criterion (AIC).

We also constructed a haplotype network using TCS 2.1.1 for preliminary exploration of the phylogeographic relationships of the obtained haplotypes (Clement, Posada & Crandall, 2000). Genetic diversity was calculated as nucleotide (π) and haplotype (h) diversity for each of the sub-regions of the basins (northern, central and southern) using Arlequin 3.5.1.2 (Excoffier & Lischer, 2010), for both pairs of concatenated genes (mitochondrial and nuclear). We also calculated the population differentiation index F_ST and its significance levels, between pairs of basins, using Arlequin 3.5.1.2 (Excoffier & Lischer, 2010). To visualize the differences in nuclear genes between basins, the genetic structure was investigated using Structure 2.3 (Pritchard, Stephens & Donnelly, 2000) with the Admixture model and correlated frequencies. A 500,000 burnin was employed and 5,000,000 generations of K 1 to 10 were performed with 10 iterations for each gene. We used the Evanno method to determine the true value of K (Evanno, Regnaut & Goudet, 2005) with the program Structure Harvester (Earl & vonHoldt, 2012). This analysis was employed only for nuclear genes, as it is not designed to find structure in linked or haploid genes.

Phenotypic, genetic and geographical distance

We performed Mantel tests and partial Mantel tests with 1000 resamplings using R v. 2.15.1 (R Development Core Team, 2011) to evaluate the correlation between acoustic (phenotypic), genetic and geographic distances. To calculate the acoustic matrix we employed the distances between canonical variables for each basin obtained from a discriminant analysis. The geographic distance was constructed with the geo- reference points taken at each locality; as the basin could include several localities, an average geo-reference was employed. The genetic distance matrix was built using the D_xy values.

The project was reviewed and approved by the Consejo Técnico de la Investigación Científica (CTIC) of UNAM.

Song description

We obtained 169 good quality recordings of male advertisement calls. These came from nine Hyla arenicolor, 20 Tlalocohyla smithii and 140 Hyla eximia from both sympatric and allopatric localities. A detailed description of the calls of H. eximia based on our complete sample is available in the Supplemental Information (see song description S1) along with descriptions of those of H. arenicolor and T. smithii which are also summarized in Tables 2 and 3. Briefly, the calls of H. eximia consist of a single note that lasts less than a quarter of a second, has a dominant frequency of about 2.6 kHz, and is composed of ca. 18 pulses. The call of Hyla arenicolor (Table 3) is roughly three times longer, (0.88 s), although it is composed of only a few more (about 20 pulses) but these pulses are longer and louder than the calls of H. eximia. Calls of H. arenicolor are lower-pitched than those of H. eximia (dominant frequency = 1.26 kHz). The structure of the call of male T. smithii is very different from those of H. eximia and H. arenicolor. Often the males emit calls with two elements; a long note followed by an extremely short one. Most individuals, however, produce more long than short notes, and in several recordings short notes are absent. Some males produced a third note with different structure and duration than the previous two, but it is not clear whether this is also part of the advertisement call. For our analyses we considered just the first note, as it is the most constant element in our recordings, and it is the one more likely to interfere with the advertisement calls of H. eximia calls. Spectrograms of typical advertisement calls of the three species are shown in Fig. 2.

Call variation

The following variables were highly correlated in the calls of Hyla eximia in allopatry: pulse amplitude and interval between calls (0.997), and call period and inter-call duration (0.993). In sympatry with H. arenicolor, the correlated variables of calls of H. eximia were: call period and inter-call duration (0.991), pulse rate and inter-pulse duration (−0.984), pulse period and inter-pulse duration (0.954), and pulse rate and pulse period (−0.953). Correlated variables in the contact zone between H. eximia and T. smithii were: call period and inter-call duration (0.999), pulse period and pulse duration (0.996), pulse rate and inter-pulse duration (−0.966), dominant frequency and pulse peak frequency (0.943), and pulse rate and number of pulses (0.920).

The following attributes of the calls of H. arenicolor were intercorrelated: pulse rate and inter-pulse duration (−0.985), pulse period and inter-pulse duration (0.965), pulse rate and pulse period (−0.941), and call period and inter-call duration (0.917). In the case of T. smithii the highest correlations amongst variables were: call period and inter-call duration (0.998), pulse period and inter-pulse duration (0.990), pulse rate and inter-pulse duration (−0.979) and pulse rate and pulse period (−0.967).

Considering that correlations ≥0.7 indicate redundant variables, we dropped from the analyses the following variables: call period (CP), pulse period (PP), pulse rate (PR), inter-pulse duration (ID) and pulse peak frequency (PPF). Consequently, analyses were performed using seven variables: call duration (CD), inter-call duration (IC), number of pulses (NP), pulse duration (PD), pulse amplitude (PA), dominant frequency (DF) and fundamental frequency (FF).

Covariance analyses revealed that temperature, chorus composition and the interaction between them affect call attributes of each species in different ways (Table S1). To control for the effect of temperature, in cases where the interaction between temperature and chorus composition was significant we used the residuals from the ANCOVA, rather than the original variable, in all subsequent analyses.

MANOVA results show that call variability in H. eximia is due to the locality (F_(6,42) = 3.91, P = 1.66e⁻¹⁴), whereas the fact that individuals call in allopatry or sympatry does not contribute significantly to call variation. However, when populations are grouped according to watershed (Panuco, Lerma, Magdalena y Balsas), it is amongst basins, rather than localities, that calls diverge significantly (F_(3,21) = 7.7241, P < 2e⁻¹⁶). The seven call variables entered in the MANOVA analyses allow for the discrimination between the four basins with a 34% reduction in classification error, classifying correctly 51% of the calls according to the basin of origin. The first two canonical functions explain 99.8% of the variance and define a multivariate acoustic space where the calls of H. eximia from Balsas and Magdalena are distinct from those from the other basins (Figs. 2 and 3). Subsequent ANOVAs and Tukey tests confirm that these differences are significant (function 1, F_3,134 = 87.67, P < 0.0001; function 2, F_3,134 = 6.84, P = 0.0002; Fig. 4). Discriminant function 1 distinguishes amongst calls from different basins based on the interval between calls (r = −0.809), call duration (r = −0.274) and number of pulses (r = 0.254), while the discriminant function 2 discriminates by dominant frequency (r = 0.974). Applying the functions that discriminate between calls of H. eximia from different basins to the call attributes (or the ANCOVA residuals) of H. arenicolor and T. smithii reveals no overlap of canonical acoustic space between the calls of H. eximia from the Balsas basin and those of co-occurring T. smithii, and a partial overlap between the calls of H. eximia at Magdalena those of the local T. smithii (Fig. 3). There is complete overlap between the canonical acoustic space of calls of H. eximia from the several populations in the Panuco basin, whether they are sympatric with H. arenicolor or not, and between those and the calls of H. eximia in the Lerma basin (Fig. 3).

Genetic variation

We found a total of 31 different haplotypes of concatenated mitochondrial genes (Table 4). Nucleotide diversity, number of analyzed samples and total haplotypes per basin for each concatenated gene pair are shown in Table 4. Mitochondrial nucleotide diversity was generally low, and the high values from Magdalena (π = 0.462) and Balsas (π = 0.269) are probably not reliable since standard deviations are high due to small sample size (Table 4). Mitochondrial haplotype diversity was very high in all regions (Table 4). Since in most cases the number of haplotypes is identical to the sample size, it does not represent the populations’ haplotype diversity.

The phylogenetic reconstruction based on mitochondrial genes recovers four major lineages of H. eximia loosely grouped by basin, with some admixture (Fig. 5). The basal clade (A) includes frogs from the southern reaches of the Panuco drainage, at the edge of the ancient closed basin of Apan, as well as frogs from the Balsas and one from Magdalena. This is the best -supported clade with bootstrap values of 56–96%. The next mitochondrial clade to diverge (B) includes frogs from central Lerma, Magdalena and Balsas basins. The most recent, sister clades include only Lerma (C) frogs, or frogs from north and central Panuco and central Lerma (D). Bootstrap values are higher in mitochondrial genes than in nuclear genes, and for the basal groups more than for the recent groups. As expected from the above reconstruction, the populations comprising clade A (southern Panuco, Balsas and Magdalena) are the most differentiated; they differ mostly (and significantly) from populations in clade D; central and northern Panuco, as well as from central Lerma as shown by the F_ST values (see Table 5).

The mitochondrial haplotype network seems to reflect the distribution of individuals into the major basins (Fig. 6). Again, haplotypes from southern Panuco are not connected to any other basin. Two of the three Balsas haplotypes are also isolated from those from other basins, except for one haplotype from Magdalena (the other Balsas haplotype is linked to haplotypes from central Lerma). By and large, the haplotypes from central and northern Panuco are closely linked, and relatively closely linked also to those from the three Lerma sub basins (Fig. 6).

We found 25 different haplotypes of concatenated nuclear genes (Table 4). The greatest nucleotide diversity of nuclear genes was found in central Panuco (π = 0.200) central (π = 0.190) and southern Lerma (0.157), but the latter has a high standard deviation (Table 4). Just as with mitochondrial genes, nuclear haplotype diversities were often ≃1 and are certainly overestimations.

The reconstruction based on nuclear genes is much less clear, with very shallow branches, although it also suggests a link between Lerma and both Balsas and Panuco (Fig. 5). The degree of nuclear genetic differentiation of H. eximia between regions is roughly consistent with that shown by mitochondrial genes: populations from southern Lerma (mitochondrial clade A) are significantly differentiated from those in northern and central Panuco, as well as central Lerma (mitochondrial clade D), yet this latter clade is formed by frogs from regions that are significantly differentiated by nuclear genes. Frogs from Magdalena (mitochondrial clade B) are significantly differentiated from those from northern Panuco (mitochondrial clade D; Table 5).

The Evanno method obtained a true value of K = 3 which is, however, not too sharply delimited. Both genes show different patterns among the populations analyzed and little correspondence to the patterns inferred from mitochondrial genes (Fig. 7). This could be due to recombination or to the fact that Structure fails to find differentiation when the population structure is too subtle (Latch et al., 2006), which could be our case as the reduced sample size failed to produce significant values of population differentiation.

Phenotypic, genetic and geographical distance

There was no significant correlation between call variation and geographic distance (Mantel test r = 0.125, P = 0.306) even when controlling for mitochondrial (Mantel partial test r = 0.126, P = 0.218) or nuclear (r = 0.185, P = 0.234), D_xy. Nuclear D_xy was not correlated with geographic distance (r = 0.097, P = 0.311). This pattern remains when controlling for call variation (r = 0.167, P = 0.293). We did find a significant correlation between geographic distance and mitochondrial D_xy (r = 0.372, P = 0.035), and this association remained when controlling for call variation (r = 0.372, P = 0.042). Call variation was not correlated with either mitochondrial (r = 0.022, P = 0.409) or nuclear (r = 0.426, P = 0.108), D_xy, even when controlling for geographic distance (mitochondrial D_xy, r = 0.026, P = 0.419; nuclear D_xy, r = 0.443, P = 0.128).

Song description

Intraspecific variation in the attributes of the advertisement calls of H. arenicolor and of the H. eximia species group is common (e.g., Blair, 1960; Klymus et al., 2010), but its origin and biological significance are not fully understood. A comprehensive study of call variation across the whole geographic range of the Canyon treefrog revealed that such differences are not large enough to promote speciation, although a degree of assortative mating preferences was evident over large geographic scales (Klymus & Gerhardt, 2012), and at least some of that variation may be linked to introgression with species in the H. eximia group (see Klymus et al., 2010). Three populations of Hyla wrightorum, a member of the latter group, were found to produce calls of different dominant frequency (and two differed in the duration of their calls; Gergus, Reeder & Sullivan, 2004). Pulse rate, an attribute often involved in mate recognition in these treefrogs (e.g., Klymus & Gerhardt, 2012), was not different between populations, thus the authors suggested that call variation in H. wrightorum is unlikely to promote mating isolation (although we note that differences in dominant frequency may be due to local adaptation to facilitate transmission (Littlejohn, 1970) and could thus lead to assortative mating).

Variation in the calls of H. eximia has been less comprehensively studied. An early report by Blair (1960) classified some populations as either producing “slow” (PR ∼50 Hz) or “fast” (PR ∼100 Hz) calls. That study was not concerned with exploring the possible origin of such variation but we now know that two populations then categorized as producing slow calls correspond to chorus where H. eximia is sympatric with T. smithii; in these we recorded songs with low PR and relatively high PPF and DF. In fact, the calls of frogs recorded in sympatry with T. smithii were more variable (mean CV = 29%) than those recorded elsewhere (see Table 2 and Fig. 2). This is unlikely to be an artifact of a small sample size, as this estimate is based on a larger sample than that used to assess variation in the attributes of calls produced in sympatry with H. arenicolor. In mixed choruses, calls of H. eximia may be masked by the long and/or the short notes of T. smithii calls. This may mean that males of H. eximia confront unpredictable note periodicity, since calls of T. smithii are made of different notes, both of variable length, and each with a different and also variable period (Table 3). This would promote high variation of inter-call interval in H. eximia, which may explain the high coefficient of variation of IC in those sympatric localities (Table 2).

We could not discern a link between type of variable (temporal, frequency, structural) and the amount of variation it harbors (Table 2; Fig. 2), but pulse peak frequency, dominant frequency, call duration and pulse duration were the least variable elements of H. eximia calls, and may thus be useful elements for species recognition.

Call variation

In his comprehensive work, Duellman (1970) and Duellman (2001) suggested that H. eximia lacks a typical call. This was puzzling since anuran advertisement calls are often stereotypical, which is why they often constitute pre-mating barriers between closely-related species. We did find a substantial degree of variation in calls, but mostly restricted to two geographic regions: Balsas and Magdalena (Fig. 3). Calls of H. eximia from these two basins were significantly different from calls elsewhere, mostly in structure/length (IC, CD, NP; discriminant function 1). Number of pulses (NP) was highly correlated with pulse rate (PR; e.g., H. eximia in contact with T. smithii: 0.9204), which together with dominant frequency plays an important role in species recognition (Littlejohn, 1965; Gerhardt & Davis, 1988; Gerhardt, 1994) and mate discrimination in other hylids (e.g., Dendropsophus ebraccatus (Hyla ebraccata), Wollerman, 1998). There are no published studies of mate choice in H. eximia, but in the related H. arenicolor female mate choice is influenced by call PR (Klymus & Gerhardt, 2012), which is a key element for maintaining isolation between lineages of H. arenicolor and might also lead, through its correlation with NP, to reproductive isolation between populations of H. eximia. Call (CD) and inter-call duration (IC) reflect calling effort, and the energy required to increase them depends on body size (Gerhardt, 1994). We found no differences between populations or basins in body size, thus it is possible that variation of call structure in the Balsas localities is adaptive, enabling H. eximia to avoid interference from the calls of T. smithii (see below).

Calls of H. eximia in Balsas and Magdalena were also different in frequency (DF; function 2) from calls elsewhere. This is a static trait (Gerhardt, 1991) that, like other such attributes of treefrog calls, is determined by morphology (several body structures, muscle size, larynx structure, among others; Dubley & Rand, 1991; McClelland, Wilczynski & Ryan, 1996; Wells, 2001; Gerhardt & Huber, 2002), resulting in structural similarity of notes and seemingly simple calls (compared to the more complex bird song). A change in this attribute would be energetically costly, and it is therefore likely that the observed variation in DF is the result of selection to improve communication efficiency, which may include avoiding interference from heterospecific calls.

In the Balsas and Magdalena basins, H. eximia shares choruses with T. smithii. Calls of H. eximia from Balsas are very variable, and their canonical acoustic space (as defined by discriminant functions 1 and 2; Fig. 3) is different from that of the local T. smithii, which seems to indicate reproductive character displacement, a possibility that requires corroboration from more intense sampling and phonotaxis studies. Interestingly, the distribution in the same canonical space of the calls of H. eximia from Magdalena shows that they are shifted in the opposite canonical direction than in Balsas. This can be caused by random differences being selected in different populations, as long as they reduce overlap with T. smithii. However, we note that in Magdalena H. eximia also faces interference with another frog species; Gastrophryne sp. Probably because this is a larger frog, its calls are much longer and lower-pitched (Duellman, 2001) than to those of H. eximia. This results in a more limited acoustic space available for H. eximia to accommodate its call to avoid masking, being simultaneously limited by the temporal and spectral attributes of the calls of T. smithii and of Gastrophryne sp.

While at Balsas and Magdalena the calls of H. eximia differ significantly from those of their conspecifics elsewhere (and are completely different from those of T. smithii in Balsas), calls from localities sympatric with H. arenicolor (mostly in the Panuco basin) were similar to those from allopatric populations. This was unexpected since, although both species are morphologically very different (H. arenicolor has a rough, granulated skin and is somewhat larger than smooth-skinned H. eximia, which is more similar to the smaller T. smithii), they belong to the same phylogenetic group, share the same mating ponds, and in the Balsas (where we did not sample mixed choruses of these two species) there is evidence suggesting ancient mitochondrial introgression (Bryson et al., 2010; Klymus & Gerhardt, 2012), thus mating interference between them was expected. It may be that, since the calls of H. eximia and H. arenicolor are substantially different, there is no, or little, interference between their calls to drive character displacement, although it may also be that the interaction of the two species is too recent in northern and central Panuco (see below).

Genetic variation

While call variation between some regions covered in our sample is supported by mitochondrial phylogeny and population structure, it cannot be properly interpreted as a consequence of the phylogeny of H. eximia. The mitochondrial reconstruction is generally well supported, although it draws a puzzling phylogenetic history of H. eximia. It suggests, very preliminarily, that the populations in our sample originated in the south-east of the current species distribution, and expanded clockwise in their colonization of what are now the basins of rivers Balsas, Magdalena (Ameca), and central (lower) Lerma, from where they may have reached the northern (Santiago) and southern (upper) Lerma, as well as the central and upper Panuco basins (Fig. 1). Published phylogeographic hypotheses for H. arenicolor often include a few sequences of H. eximia and cover a limited portion of its geographic range, so that no phylogeographic hypotheses for this species can be drawn from them (see Bryson et al., 2010; Klymus & Gerhardt, 2012). Thus we provisionally propose the above phylogenetic scenario of H. eximia based on our best -resolved (mitochondrial) reconstruction (Fig. 1).

Phenotypic, genetic and geographical distance

Call variation between the regions covered in our sample may not be random (i.e., due to drift); it was greater in particular basins (Balsas and Magdalena) and chorus compositions (H. eximia–T. smithii), but did not correlate with genetic or geographic distance. Absence of a correlation between acoustic and genetic distance has been reported in other studies. For instance, in Oophaga (Dendrobates) pumilio (Pröhl et al., 2007), the correlation is lost when the analysis controls for geographic distance, and a similar effect occurred in some populations of the Túngara frog Engystomops (Physalaemus) pustulosus (Ryan, Rand & Weigt, 1996). These findings are perhaps not surprising, since one of the forces that most effectively shape anuran call attributes is female mate choice (Boul et al., 2007), and this is most effective when mistakes are penalized as in a hybrid zone. Thus we should in fact expect greater call variation between populations when reproductive character displacement is favored in at least one of them (Pfennig and Pfennig, 2010), than when populations are distant geographically but do not face mating call interference. Indeed, we rule out the possibility that a correlation between genetic and acoustic distance exists but was undetected, since we have adequate call samples, and, as expected, mitochondrial genetic distance was correlated with geographic distance, making it likely that our genetic screening is unbiased.

We have argued that call variation in the regions where H. eximia and T. smithii coexist is consistent with the expected effects of acoustic interference, but it might also be the result of morphological differentiation in attributes that influence call production (e.g., the body size in Hyla leucophyllata, Lougheed et al., 2006; see also Emerson, 2001; Gerhardt & Huber, 2002; Castellano et al., 2002; Hoskin et al., 2005). However, an exploratory analysis failed to detect significant differences in body size among populations of H. eximia, thus it is unlikely that the variation in calls found between basins is due to differences in the body size. Other causes of call differentiation include local differences in call assemblage (Wollerman & Wiley, 2002), in signal transmission efficiency (related to the type of habitat), or in sexual selection (acting through reproductive isolation mechanisms). The first is compatible with our explanation of call differentiation to avoid masking by the calls of T. smithii, whereas the reported lack of differentiation amongst the populations not sympatric with T. smithii, which cover a much larger and presumably environmentally variable geographic area, argues against the other two. We thus provisionally propose that the geographic variation of advertisement calls of H. eximia in the Balsas and Magdalena basins populations is due to reproductive character displacement produced by sharing mating choruses with the related and morphologically similar T. smithii and with the distantly related Gastrophryne sp.