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

Populations of species occupying large geographic ranges are often phenotypically diverse as a consequence of variation in selective pressures and drift. This applies to attributes involved in mate choice, particularly when both geographic range and breeding biology overlap between related species. This condition may lead to interference of mating signals, which would in turn promote reproductive character displacement (RCD). We investigated whether variation in the advertisement call of the mountain treefrog (Hyla eximia) is linked to geographic distribution with respect to major Mexican river basins (Panuco, Lerma, Balsas and Magdalena), or to coexistence with its sister (the canyon treefrog, Hyla arenicolor) or another related species (the dwarf treefrog, Tlalocohyla smithii). We also evaluated whether call divergence across the main river basins could be linked to genetic structure. We found that the multidimensional acoustic space of calls from two basins where H. eximia currently interacts with T. smithii, was different from the acoustic space of calls from H. eximia elsewhere. Individuals from these two basins were also distinguishable from the rest by both the phylogeny inferred from mitochondrial sequences, and the genetic structure inferred from nuclear markers. The discordant divergence of H. eximia advertisement calls in the two separate basins where its geographic range overlaps that of T. smithii can be interpreted as the result of two independent events of RCD, presumably as a consequence of acoustic interference in the breeding choruses, although more data are required to evaluate this possibility.


INTRODUCTION
Populations of species occupying large geographic ranges are likely to experience different selective pressures (West-Eberhard, 1983;Panhuis et al., 2001;Coyne & Orr, 2004) which, together with drift (Wiens, 2004), may result in phenotypic and genotypic differences between populations (e.g., Avise, 2000;Laugen et al., 2003;Amézquita et al., 2009). Local differences in ecology, such as prey availability (e.g., Arnold, 1980;Bonansea & Vaira, 2007) or habitat structure (e.g., Relyea, 2002;Skelly, 2004) can lead to differential adaptation between populations (Newman, 1992), but in species with generalized diets and habitat requirements such variation would have a limited effect on differentiation (e.g., Virgós, Llorente & Cortés, 1999). On the other hand, environmental variation in the attributes that determine the transmission and reception of signals used in social (Wells, 1977;Sullivan & Wagner, 1988;Wagner, 1989) or sexual contexts (Searcy & Andersson, 1986) may lead to rapid population differentiation (Jennions & Petrie, 1997;Seehausen, van Alphen & Witte, 1997;Lougheed et al., 2006;Boul et al., 2007). This is because signals used by animals in a breeding context may convey information about the species, sex, breeding status, and even the condition of the sender (Gerhardt, 1992;Wilczynski & Chu, 2001), all of which may be relevant to conspecifics searching for mating partners (Butlin & Ritchie, 1994;Emerson, 2001;Forsman & Hagman, 2006). Geographic variation in mating signals has been widely reported in studies of character displacement, where specific traits (e.g., morphological, behavioral, ecological or physiological traits) differ among sympatric and allopatric populations due to the risk of maladaptive hybridization with related species (Brown & Wilson, 1956;Grant, 1972). Reproductive character displacement (RCD) is often studied in systems where sister species meet at typically narrow hybrid or tension zones (Butlin, 1987;Howard, 1993;Butlin & Ritchie, 1994). However an approach involving the study of signal variation across wide geographic areas, coupled with genetic and geographic data, could be helpful in identifying evolutionary patterns of signal evolution and tracing the links between micro and macro evolution of mating signals (Avise et al., 1987).
Amphibians are good models for studying signal evolution and reproductive isolation since their dispersal is restricted by habitat, and they thus are highly philopatric, facilitating the study of genetic divergence amongst populations (c.f. Gamble, MacGarigal & Compton, 2007). Amphibian philopatry results in populations being genetically structured over short geographic distances, which facilitates inferences about the historic events that caused their current distribution (Zeisset & Beebee, 2008). Although other modalities (e.g., visual;Hartmann et al., 2005;Reynolds & Fitzpatrick, 2007;Taylor et al., 2008) are sometimes involved, it is common in anurans that species recognition and female preference depend on a single acoustic signal; the advertisement call Wells, 1977;Cocroft & Ryan, 1995;Wells & Schwartz, 2007. This is produced by males, typically during the mating season. Several attributes of the calls vary within a species as a function of male morphology and condition. For instance, the frequency or pitch (Hertz) of the call is inversely correlated with body size (Zug, 1993;Gerhardt & Huber, 2002). Calls are also affected by temperature, which influences the length of the call as well as the number of pulses that compose the call, and the pulse rate (Gerhardt & Davis, 1988;Wilczynski & Chu, 2001;Gerhardt & Huber, 2002;Yamaguchi et al., 2008). Additionally, males at a chorus often adjust the timing of their calls to avoid masking by other vocalizations, and when choruses are formed by more than one species, acoustic interference can select for the use of different acoustic space (Brush & Narins, 1989;Schwartz, Buchanan & Gerhardt, 2002;Chek, Bogart & Lougheed, 2003;Hoskin & Higgie, 2010). Acoustic space is a multidimensional space defined by acoustic properties of the signals, such as frequency, duration, amplitude, temporal structure, etc.; see, for instance, Ryan & Rand (2003). In spite of variation, anuran songs are sufficiently stereotyped to permit females to identify and assess males (e.g., the Hylidae; Gerhardt, 1991;Gerhardt, 1992), so much so that Hyla spp. have often been used as models to study pre-mating reproductive isolation amongst sister species (e.g., Littlejohn, 1965;Ball & Jameson, 1966;Littlejohn, 1970;Duellman, 1973;Gerhardt, 1994;Höbel & Gerhardt, 2003;Gerhardt, 2005).
Hyla eximia is a relatively small (ca. 3 cm) tree-frog endemic to the geologically complex Trans-Mexican Volcanic Belt (TMVB;de Cserna, 1989;Ferrari et al., 2012) and adjacent Mexican High Plateau (Duellman, 1970;Duellman, 2001). It gives its name to the H. eximia species group (Smith et al., 2007), which includes eleven other species (Faivovich et al., 2005). The very similar H. wrightorum replaces H. eximia to the north (through Arizona), and H. euphorbiacea is found to the south of the TMVB (allopatric with H. eximia). H. plicata is endemic to the TMVB and is partially sympatric with H. eximia (Smith et al., 2007), and H. arenicolor is sympatric with H. eximia in parts of the Mexican Plateau. Like H. eximia, H. arenicolor has a wide geographic range; in the Mexican Plateau it is mainly found in the mountain zones of Balsas Basin, at elevations of 300-3000 m (Duellman, 2001). Tlalocohyla smithii (formerly Hyla smithii, Faivovich et al., 2005) is distributed in the Pacific lowlands of Mexico up to elevations of about 1000 m, from central Sinaloa to southern Oaxaca, and inland within the Balsas Basin (Duellman, 2001). Tlalocohyla smithii, is not part of the H. eximia group, but occupies the same matings pools.
Previous work has shown evidence of geographic variation in H. eximia calls. Indeed, variation in pulse rate, call duration and the dominant frequency among recordings from several populations (Blair, 1960;Duellman, 1970;Sullivan, 1986;Duellman, 2001) was so great that Duellman (1970) andDuellman (2001) concluded that the species lacks a typical call, though sample sizes per population were small. It has also been suggested that some variation in advertisement calls of H. eximia may be linked to RCD. Based on phonetic data Cortés-Soto (2003) suggested that H. eximia and H. plicata have evolved different advertisement calls in the 500 m altitudinal band (2400-2900 m asl) where their ranges overlap. There, the calls of H. eximia are shorter and contain fewer pulses than when the species are found in allopatry, suggesting that RCD has occurred. Here we describe the variation of the advertisement calls of H. eximia across a substantial part of its geographic range, and explore whether this variation is linked to genetic structure (based on mitochondrial and nuclear DNA sequences), geography (i.e., hydrographic basins) and/or range overlaps with its sister species, the canyon treefrog (H. arenicolor), and a related species, the Mexican dwarf treefrog (Tlalocohyla smithii).

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  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 TM ME66) connected to a digital recorder (Marantz TM 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 TM 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 TM 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 The distribution of frogs belonging to clades A-D from a phylogenetic analysis of two concatenated mitochondrial genes (see text) is indicated with globes of different color (clade A, turquoise blue; clade B, red; clade C, yellow; D, green). Dotted arrows indicate the possible direction of the colonization, and are not intended to show the exact paths followed by the frogs. The occurrence of one individual from Magdalena in clade A is indicated with a dotted turquoise-blue line. Thin grey dotted lines represent boundaries between the major basins (Panuco, Lerma, Balsas and the currently closed Valley of Mexico), whereas rivers are indicated with continuous grey lines. Lerma River flows northwards until it meets the Santiago, which drains in the Pacific Ocean to the west. Balsas drains southwards, also in the Pacific, whereas the Panuco drains to the east, in the Gulf of Mexico. See Table 1 for population codes. 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 harmonicsthat 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

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 TM (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 TM 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.
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 −14 ), 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 −16 ). 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).

Basin Northern panuco Central panuco Southern panuco Central lerma Southern lerma Magdalena Balsas
Northern

Notes.
Mitochondrial F ST below, and nuclear F ST above the diagonal.
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) andDuellman (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 closelyrelated 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.