Molecular Ecology (2005) doi: 10.1111/j.1365-294X.2005.02707.x ? 2005 Blackwell Publishing Ltd mBlackwell Publishing, Ltd. Biogeography of the t?ngara frog, Physalaemus pustulosus : a molecular perspective LEE A. WEIGT, * ANDREW J . CRAWFORD, ? A . STANLEY RAND ? and MICHAEL J . RYAN ?? * Smithsonian Institution, National Museum of Natural History, Washington, D.C. 20250, ? Smithsonian Tropical Research Institute, Apartado 2072, Balboa, Anc?n, Republic of Panama, ? Section of Integrative Biology, 1 University Place C0930, University of Texas, Austin, TX 78712 Abstract Physalaemus pustulosus , a small leptodactylid frog with South American affinities, ranges across northern South America through Middle America to southern Mexico. To investi- gate its geographic variation and evolutionary origins, we analysed the presumptive gene products of 14 allozyme loci and sequenced a portion of the mitochondrial COI gene from individuals sampled throughout the distribution. Generally, allozyme dissimilarities and sequence divergences are correlated with each other and with geographic proximity. The greatest discontinuity in genetic variation was found between populations in Middle America vs. South America + Panama. Based on two Bayesian MCMC (Markov chain Monte Carlo) divergence time estimates involving two independent temporal constraints, the timing of the separation of northern and southern t?ngara frog lineages is significantly older than the time since completion of the current Panama land bridge. P . pustulosus first invaded Middle America from South America about 6?10 million years ago giving rise to the northern lineage. The southern lineage then invaded Panama independently after land bridge comple- tion. Despite millions of years of independent evolution, the multilocus allozyme data revealed that western Panama populations represent a contact zone containing individuals with alleles from both groups present. Keywords : allozymes, biogeography, intergradation, mtDNA Received 26 May 2005; revision accepted 25 July 2005 Introduction The t?ngara frog, Physalaemus pustulosus , ranges from Veracruz and Yucatan in southern Mexico down the Pacific coast of Middle America and across the Caribbean coast of northern South America to Trinidad and Guyana. Many dry- forest amphibians and reptiles show a similar distribution in Middle America and northern South America (Dunn 1940). The family Leptodactylidae, subfamily Leptodacty- linae, and genus Physalaemus are largely South American (Savage 2002). The genus Physalaemus contains at least 45 species (Frost 2004). All but P . pustulosus are confined to South America and most species in this genus occur east of the Andes and south of Amazonia. Early work on the species focused on mate recognition signals, female mating preferences, sexual selection, and life history (see Ryan 1985; Ryan & Rand 2003 for reviews). Recent investigations include genetic studies of relationships of several members of the genus (Cannatella et al . 1998), geographic variation in mating calls (Ryan et al . 1996), and microsatellite primer development (Pr?hl et al . 2002). All of these studies assume P . pustulosus is a single taxonomic entity. The last decade has seen a great increase in our under- standing of the geological events culminating in the rise of the Isthmus of Panama and its effect on the biota. Recon- structed from diverse lines of evidence, the emergence took millions of years, with the seaway closing with the surfacing of the Panama land bridge in late Pliocene (Coates et al . 1992, 2004; Farrell et al . 1995; Coates & Obando 1996; D?az de Gamero 1996 ? see Table 5, p. 744 for chronology of events). Movements of animals between North and South America in the past have been complex and their interpretation con- troversial (Stehli & Webb 1985). Vanzolini & Heyer (1985) suggest that the emergence of the Panama land bridge did not produce the sort of massive interchange between South Correspondence: Lee A. Weigt, Fax: 301 238 3059; E-mail: weigtl@si.edu 2 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology , 10.1111/j.1365-294X.2005.02707.x and North American herpetofaunas that is described for mammals (Webb 1985). They assert that most of the her- petofaunal interchange, including the leptodactylid frogs, took place prior to the emergence of the land bridge in Panama. Duellman (1966) and Savage (1966, 1982) express a contrasting view. They include P . pustulosus among the species that they think invaded Middle America from South America after the present Panama land bridge emerged. Modern-day western Panama populations of P . pustulosus are separated from the rest of the Central American popu- lations by a gap in distribution in southwestern Costa Rica (Savage 2002). In the context of the above geological history, we examined nuclear allozyme data and mitochondrial DNA (mtDNA) sequences for one species of frog that occurs in both South and Middle America. Molecular data provide additional information, along with data on present-day distributions, the fossil record, and morphology for reconstructing migra- tion patterns and vicariant events (e.g. Bermingham & Avise 1986). Such reconstructions not only provide information about the history of the taxa under study, but also constrain interpretations of the geological history of the region. Reconstructing the geographic origins of a given species requires knowledge of the phylogenetic relationship of that species with its nearest relatives. In their analysis of relationships among species within the pustulosus group, Cannatella et al . (1998) found that the inferred phylo- genetic position of P . pustulosus differed among data sets. The COI data suggested that P . pustulosus is sister to the other six species of the pustulosus species group, while the morphology, 12S mitochondrial DNA sequences, allozyme data and the combined data all favoured P . pustulosus as the sister to the pair, Physalaemus petersi + Physalaemus freibergi . Among these data sets, only the morphology gave strong support to a particular sister relationship for P . pustulosus . In the present study, we revisit the relationships among species, and among populations of P . pustulosus . Materials and methods Collection sites and sampling Figure 1 and Table 1 show the collecting localities sampled throughout the range of Physalaemus pustulosus . Ryan et al . Fig. 1 Map showing collecting localities for Physalaemus pustulosus (circled dots) and the mountain ranges whose origins are used to constrain divergence time estimates. Each locality is labelled with the same four-letter code used in Table 1. A close-up map of the Panama Canal Zone localities (inset) is shown on the lower left. In one divergence time analysis, the Ecuadorian Andes are assumed to have split the four Pacific coast Physalameus from the two Amazonian species plus P. pustulosus. In a second, independent analysis the northern Andes, comprised of the Serran?a de Perij? and M?rida Andes, are assumed to have divided the CALA, LMAR, and MARI populations at least 2 million years ago. Note, Physalaemus is absent from the lowland rainforest of the Choc?. Elevation is indicated by five different shades, with the darkest shade representing 1?500 m and the lightest shade representing elevations above 2000 m. B IO G E O G R A P H Y O F T ? N G A R A F R O G S 3 ? 2005 B lackw ell P u blishing L td , M olecular E cology , 10.1111/ j.1365-294X .2005.02707.x Table 1 Populations, locality information and sample sizes of 30 populations of Physalaemus pustulosus analysed in this study. Right-hand columns show biosys -1 output statistics on genetic variability within populations (standard errors in parentheses) Population Locality information Sample size Mean sample size/locus Mean no. of alleles/locus Percentage of polymorphic loci Mean heterozygosity (left, direct count; right, Hardy?Weinberg expected)(Latitude) (Longitude) VERC Laguna Verde, Veracruz, Mexico 11 10.9 (0.1) 1.2 (0.10) 21.40 (0.06) 0.10 (0.06) 0.10 19.73 ? 96.43 TEHU Tehuantepec, Chiapas, Mexico 11 11.0 (0.0) 1.4 (0.20) 35.70 (0.05) 0.12 (0.05) 0.13 16.35 ? 95.28 TAPA Tapachula, Chiapas, Mexico 16 15.0 (0.1) 1.5 (0.20) 35.70 (0.05) 0.11 (0.05) 0.11 14.86 ? 92.22 GUAT Taxisco, Santa Rosa, Guatemala 10 10.0 (0.0) 1.4 (0.20) 35.70 (0.06) 0.11 (0.04) 0.10 14.03 ? 90.18 ESAL San Miguel, El Salvador 12 12.0 (0.0) 1.3 (0.10) 28.60 (0.05) 0.10 (0.04) 0.08 13.49 ? 88.18 NICA Tipitapa, Managua, Nicaragua 12 12.0 (0.0) 1.4 (0.20) 28.60 (0.03) 0.07 (0.03) 0.06 12.20 ? 86.07 CRIC Liberia, Guanacaste, Costa Rica 11 11.0 (0.0) 1.4 (0.10) 35.70 (0.06) 0.11 (0.05) 0.09 10.61 ? 85.45 PARM Puerto Armuelles, Chiriqui, Panama 15 14.9 (0.1) 1.5 (0.10) 50.00 (0.07) 0.14 (0.04) 0.12 8.27 ? 82.86 GUAL Gualaca, Chiriqui, Panama 14 14.0 (0.0) 2.0 (0.30) 64.30 (0.07) 0.23 (0.06) 0.20 8.53 ? 82.29 SANT Santiago, Veraguas, Panama 16 15.4 (0.6) 1.8 (0.30) 42.90 (0.06) 0.19 (0.06) 0.17 8.13 ? 80.98 ANTO Anton, Cocle, Panama 15 15.0 (0.0) 1.6 (0.20) 42.90 (0.05) 0.13 (0.05) 0.12 8.40 ? 80.24 GATW Gatun West Bank, Col?n, Panama 14 14.6 (0.2) 1.4 (0.10) 42.90 (0.07) 0.18 (0.06) 0.17 9.25 ? 79.95 GATE Gatun East Bank, Col?n, Panama 15 14.9 (0.1) 1.5 (0.20) 42.90 (0.06) 0.18 (0.06) 0.19 9.28 ? 79.92 BCIC Barro Colorado Island, Col?n, Panama 14 13.9 (0.1) 1.4 (0.20) 28.60 (0.06) 0.11 (0.05) 0.11 9.17 ? 79.85 PRPH Pipeline Road, Col?n, Panama 15 14.9 (0.1) 1.5 (0.10) 50.00 (0.05) 0.13 (0.05) 0.12 9.16 ? 79.73 GAMB Gamboa, Col?n, Panama 30 29.6 (0.3) 1.6 (0.10) 57.10 (0.06) 0.13 (0.05) 0.12 9.12 ? 79.70 4 L . A . W E IG T E T A L . ? 2005 B lackw ell P u blishing L td , M olecular E cology , 10.1111/ j.1365-294X .2005.02707.x GBRG Gamboa Bridge, Col?n, Panama 16 16.0 (0.0) 1.5 (0.20) 35.70 (0.05) 0.12 (0.06) 0.14 9.11 ? 79.69 SUMM Summit Gardens, Panama, Panama 13 12.9 (0.1) 1.5 (0.10) 50.00 (0.06) 0.16 (0.06) 0.15 9.07 ? 79.65 CHIV Chiva Chiva Road, Panama, Panama 15 14.8 (0.2) 1.8 (0.20) 64.30 (0.06) 0.17 (0.05) 0.17 9.02 ? 79.59 COCO Coco Solo, Panama, Panama 14 15.0 (0.0) 1.6 (0.20) 50.00 (0.06) 0.15 (0.05) 0.15 8.97 ? 79.59 KOBE Ft. Kobbe, Panama, Panama 14 13.8 (0.2) 1.4 (0.20) 35.70 (0.06) 0.13 (0.05) 0.12 8.90 ? 79.59 ITAB Isla Taboga, Panama, Panama 14 14.9 (0.1) 1.3 (0.10) 28.60 (0.04) 0.07 (0.04) 0.07 8.80 ? 79.55 IREY Isla del Rey, Las Perlas, Panama 15 15.0 (0.0) 1.4 (0.10) 35.70 (0.03) 0.05 (0.04) 0.07 8.45 ? 78.85 METE Metet?, Darien, Panama 15 14.9 (0.1) 1.5 (0.20) 42.90 (0.08) 0.17 (0.06) 0.13 8.50 ? 77.97 ELRE El Real, Darien, Panama 15 14.4 (0.6) 1.5 (0.10) 50.00 (0.08) 0.18 (0.06) 0.15 8.13 ? 77.73 MARI Mariquita, Tolima, Colombia 15 14.1 (0.9) 1.6 (0.30) 35.70 (0.09) 0.22 (0.07) 0.18 5.18 ? 74.90 LMAR L. Maricaibo, Venezuela 2 2.0 (0.0) 1.2 (0.10) 21.40 (0.06) 0.11 (0.06) 0.11 8.56 ? 71.63 CALA Calobozo, Guarico, Venezuela 12 11.9 (0.0) 1.5 (0.20) 42.90 (0.07) 0.13 (0.04) 0.10 8.98 ? 67.35 CARU Car?pano, Sucre, Venezuela 14 13.6 (0.2) 1.4 (0.20) 28.60 (0.06) 0.12 (0.06) 0.14 10.64 ? 63.22 TRIN Trinidad 5 3.7 (0.5) 1.5 (0.20) 35.70 (0.09) 0.23 (0.07) 0.18 10.63 ? 61.28 Population Locality information Sample size Mean sample size/locus Mean no. of alleles/locus Percentage of polymorphic loci Mean heterozygosity (left, direct count; right, Hardy?Weinberg expected)(Latitude) (Longitude) Table 1 Continued B I O G E O G R A P H Y O F T ? N G A R A F R O G S 5 ? 2005 Blackwell Publishing Ltd, Molecular Ecology , 10.1111/j.1365-294X.2005.02707.x (1996) reported and compared geographic call dissimilarity and Nei unbiased genetic distances ( D N ) for these 30 popu- lations. The four-letter population abbreviations used here follow Ryan et al . (1996). Within a locality, male frogs, and occasionally females, were sampled from sites within an area radius of less than 300 m. Material from Trinidad came from Smithsonian Museum of Natural History collections. For this study one individual, of the 15 animals per site used for allozymes, was dissected in the field and flash frozen in liquid nitrogen, or prepared in the laboratory, and stored frozen at ? 80 ? C for mtDNA analyses. Molecular methods Allozyme methods and results presented in Ryan et al . (1996) are not repeated here ? only new analytical results and their methods are presented. Buffers and products resulting from standard starch gel electrophoresis (Murphy et al . 1990) are reported in Table 2. We consider the allozyme data to reflect genome-wide (nuclear) gene differences (Garcia-Paris et al . 2003). Allozyme frequency data were analysed with the biosys -1 (Swofford & Selander 1981), and multidimensional scaling was performed using ntsys - pc version 1.8 (Rohlf 1993). Maximum-parsimony phylo- genetic analysis of allozyme frequencies was calculated with freqpars (Swofford & Berlocher 1987) and paup * 4.1b (Swofford 1998). Mantel tests were used to determine correlations between all matrices ( mantel version 2.0, Mantel nonparametric calculator, Copyright 1999, Adam Liedoff). The probability of rejecting the null hypothesis was based on 1000 randomization simulations. The new comparisons not presented in Ryan et al . (1996) were those involving the mtDNA divergence matrix. For mtDNA analysis, thigh muscle was minced in DNA storage buffer (Seutin et al . 1991). Genomic DNA was extracted from a 100 mg wet-weight portion of tissue with the Gnome DNA extraction kit (BIO 101)) with all kit volumes reduced to one-third volume. The final resuspen- sion was in 100 ? L of sddH 2 O. Polymerase chain reaction (PCR) thermal cycler profiles, methods, and primers are listed in Appendix II. Double-strand sequences were unambiguously aligned with sequencher (Gene Codes). In addition to the 30 P . pustulosus samples obtained here, the final data set also included nine (one sample from each) of the following members of the P . pustulosus species group: Physalaemus colouradorum , Physalaemus pustulatus (Ecuador) P . sp. B (?caicai?), and P . sp. C, Physalaemus petersi and Physalaemus freibergi (Cannatella et al . 1998). These latter two species comprise the sister clade to P . pustulosus in all data sets of Cannatella et al . (1998) except COI. We also included as outgroups one COI sequence each from Physalaemus enesefae, Physalaemus ephippifer, and P. sp. A (Roriama, Brazil). Choice of outgroup was based on Cannatella & Duellman (1984) and Cannatella et al. (1998). Name usage follows Ryan et al. (1996) and Cannatella et al. (1998), but names of P. pustulatus from west of the Andes are in the process of being changed by Ron et al. (2004; 2005). Analytical methods The 39 aligned COI sequences were subjected to the following phylogenetic analyses. Fitch (1971) parsimony (FP) trees were inferred from heuristic searches using paup* 4.0b10 for Unix (Swofford 1998). Starting trees were obtained by stepwise addition of taxa in random order in 50 000 replicate searches, with a maximum of 5000 trees used for tree-bisection?reconnection (TBR) branch swap- ping in each replicate. Clade support in parsimony trees was evaluated by two methods. First, we employed the bootstrap procedure (Felsenstein 1985) using 2000 replicates with five random addition sequences and TBR branch swapping in each replicate, with a 1000-tree limit placed on each replicate. Second, we calculated decay indices (Bremer 1988) for each node using treerot version 2c (Sorenson 1999) and paup*. Maximum-likelihood (ML) inference began with a test of significant heterogeneity in base frequencies in the COI data across n taxa using a ?2[3(n ? 1)] test. Fifty-six models of sequence evolution were then evaluated using three criteria: by hierarchical likelihood-ratio test (hLRT; Felsenstein 1981) and Akaike information criteria (AIC; Akaike 1974), both implemented in modeltest version 3.5 (Posada & Crandall 1998), and by the Bayesian information criterion as implemented in the program dt-modsel (version 13-Aug-02) by Minin et al. (2003). This latter tech- nique tends to suggest simpler models than the former two criteria. Because branch lengths are an integral part of phylogeographic analyses, we would prefer that model suggested by the DT method over the hLRT or AIC model should they differ. In all cases, ML models were evaluated and parameters calculated from neighbour-joining (NJ) trees (Saitou & Nei 1987). In modeltest, this was based on JC69 distances (Jukes & Cantor 1969), whereas dt-modsel used an NJ tree based on LogDet distances (e.g. Lockhart et al. 1994). We also re-ran modeltest starting from a con- sensus topology of the most parsimonious trees. ML inference also used heuristic searches and TBR swapping, but started from NJ trees. ML analyses used the model/s suggested by the above criteria and assumed the given parameter values. These analyses were repeated using full searches in which both the topology and the parameter values for the given model/s were estimated. We tested for evolutionary rate heterogeneity across the phylogeny by a LRT test comparing the difference in support between trees with and without a molecular clock enforced. Rather than attempt a bootstrap analysis of clade support on ML trees, we conducted Bayesian MCMC 6 L . A . W E IG T E T A L . ? 2005 B lackw ell P u blishing L td , M olecular E cology, 10.1111/ j.1365-294X .2005.02707.x Table 2 Allele frequencies of polymorphic loci in 30 populations of Physalaemus pustulosus. Three gene products (ME, SOD1, SOD2) were invariant for all populations. Horizontal shading indicates (non-rare) alleles showing intergradation and vertical double lines show the western Panama populations where the intergradation occurs Gene product Sample size Populations Allele VERC TEHU TAPA GUAT ESAL NICA CRIC PARM GUAL SANT ANTO GATW GATE BCIC PRPH GAMB GBRG SUMM CHIV COCO KOBE ITAB IREY METE ELRE MARI LMAR CALA CARV TRIN (N) 11 11 16 10 12 12 11 15 14 16 15 14 15 14 15 30 16 13 15 15 14 14 15 15 15 15 2 12 14 5 Phosphogluconate dehyrogenase 6PGDH A 0.03 B 0.14 0.07 C 0.81 0.63 0.67 0.46 0.73 0.80 0.83 0.62 0.58 0.46 0.63 0.58 0.03 1.00 0.54 0.50 0.60 D 0.04 0.04 0.30 E 0.18 0.04 0.06 0.33 0.33 0.54 0.27 0.20 0.17 0.38 0.42 0.54 0.37 0.42 0.97 0.32 0.50 0.75 0.04 0.50 F 0.07 0.25 0.96 1.00 0.50 G 1.00 0.82 0.97 0.90 0.79 0.75 0.36 0.93 0.96 0.13 0.03 H 0.03 0.10 0.21 0.17 0.64 Aspartate transaminase ATA A 1.00 1.00 1.00 1.00 0.88 1.00 1.00 0.82 0.82 0.34 0.37 0.30 0.23 0.03 0.12 0.16 0.08 0.23 0.17 0.18 0.03 0.03 0.10 0.05 B 0.03 0.07 0.20 C 0.12 0.18 0.18 0.59 0.63 0.70 0.77 1.00 0.97 0.88 0.81 0.92 0.77 0.83 0.75 1.00 1.00 0.97 0.97 0.90 0.75 0.86 0.46 D 0.06 0.25 0.09 0.54 0.80 Triose phosphate isomerase TPI A 0.50 0.03 B 0.27 0.50 0.17 C 0.22 0.03 D 0.50 0.86 0.84 0.95 1.00 1.00 1.00 0.17 0.39 0.07 0.07 E 0.20 F 0.03 0.30 0.77 0.32 0.25 0.30 0.56 0.54 0.53 0.57 0.36 0.17 0.20 G 0.05 H 0.14 0.13 0.83 0.61 0.72 0.93 0.70 0.23 0.68 0.75 0.70 0.44 0.46 0.40 0.37 0.64 0.83 0.80 0.73 0.50 1.00 1.00 0.83 1.00 0.80 I 0.06 Lactate dehydrogenase LDH1 A 0.07 B 0.23 0.12 0.55 0.04 0.90 0.86 0.70 0.90 1.00 1.00 1.00 0.97 0.95 1.00 1.00 1.00 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 C 0.41 0.73 0.88 0.45 1.00 0.96 1.00 0.10 0.14 0.30 0.03 0.03 0.05 0.03 D 0.59 0.04 LDH2 A 0.14 0.03 B 0.04 0.03 0.03 C 1.00 0.86 0.66 1.00 1.00 1.00 0.96 1.00 1.00 1.00 1.00 0.97 1.00 1.00 1.00 1.00 1.00 1.00 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 D 0.31 B IO G E O G R A P H Y O F T ? N G A R A F R O G S 7 ? 2005 B lackw ell P u blishing L td , M olecular E cology, 10.1111/ j.1365-294X .2005.02707.x Phosphoglucomutase PGM A 0.08 B 0.05 0.08 0.04 0.03 C 0.04 D 1.00 1.00 1.00 0.90 0.92 0.92 0.96 0.90 0.96 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.96 0.97 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.25 0.04 0.36 E 0.05 0.04 0.10 0.75 0.96 0.64 1.00 Glucose phosphate isomerase GPI A 0.03 0.10 B 0.57 0.12 0.20 C 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.07 0.50 1.00 1.00 1.00 1.00 1.00 1.00 0.98 1.00 1.00 0.90 1.00 1.00 1.00 0.97 1.00 0.80 0.40 1.00 0.88 1.00 0.60 D 0.04 0.10 E 0.04 0.02 0.03 0.20 0.10 F 0.93 0.39 G 0.04 Malate dehyrogenase MDH1 A 0.04 B 0.03 0.25 C 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.54 0.53 0.93 0.90 0.98 0.97 0.85 0.87 0.97 0.93 0.96 0.97 0.97 0.93 1.00 1.00 1.00 0.71 1.00 D 0.46 0.47 0.07 0.10 0.02 0.03 0.15 0.10 0.03 0.07 0.03 0.03 0.07 MDH2 A 0.23 0.06 B 0.36 0.17 0.17 0.03 0.03 0.40 C 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.64 0.88 0.77 1.00 0.77 1.00 0.93 0.97 1.00 0.92 0.97 0.83 1.00 1.00 0.97 0.93 0.97 1.00 1.00 1.00 1.00 D 0.03 E 0.06 0.07 0.03 0.08 0.03 0.07 0.60 F 0.03 Isocitrate dehydrogenase ICD1 A 0.03 B 0.25 0.50 0.25 0.20 0.29 0.13 0.36 0.47 0.14 0.78 0.97 0.29 0.60 0.68 0.57 0.55 0.53 0.62 0.87 0.83 0.96 0.37 0.27 0.50 0.57 0.50 0.46 0.36 0.50 C 0.07 0.03 D 0.75 0.50 0.75 0.80 0.71 0.87 0.64 0.53 0.71 0.03 0.71 0.40 0.29 0.43 0.45 0.44 0.38 0.13 0.17 0.04 0.63 0.73 0.50 0.43 0.50 1.00 0.54 0.64 0.50 E 0.07 0.22 ICD2 A 0.04 0.04 B 1.00 1.00 1.00 1.00 1.00 1.00 0.96 1.00 0.96 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 Gene product Sample size Populations Allele VERC TEHU TAPA GUAT ESAL NICA CRIC PARM GUAL SANT ANTO GATW GATE BCIC PRPH GAMB GBRG SUMM CHIV COCO KOBE ITAB IREY METE ELRE MARI LMAR CALA CARV TRIN (N) 11 11 16 10 12 12 11 15 14 16 15 14 15 14 15 30 16 13 15 15 14 14 15 15 15 15 2 12 14 5 Table 2 Continued 8 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x (Markov chain Monte Carlo) analyses of phylogenetic relationships (Rannala & Yang 1996; Yang & Rannala 1997) using mrbayes version 3.0b4 for Unix (Huelsenbeck & Ronquist 2001). Each MCMC analysis involved six chains with default heating. After a conservative burn-in period of 1001 generations, 9000 trees were sampled at a frequency of one per 200 generations. To estimate divergence times, we employed the highly parametric Bayesian MCMC method of Thorne et al. (1998; Kishino et al. 2001) as implemented in the software package multidistribute version 05-Aug-03 (Thorne & Kishino 2002). In applying this method to the pustulosus group, we compared and contrasted two different sets of assump- tions based on two different geological events that may reasonably be considered independent. First, following the preferred phylogeny of Cannatella et al. (1998), we assumed that the rise of Ecuadorian Andes separated the ancestor of the Pacific coast species (P. coloradorum, P. pustulatus, P. sp. B ?caicai?, and P. sp. C) from the ancestor of the Amazonian species (P. freibergi, P. petersi) + P. pustulosus. We further assumed that this vicariant event happened at some point in the mid-Miocene, 16.4?11.2 million years ago (Ma) (Hoorn 1993; Hoorn et al. 1995; Gregory-Wodzicki 2000). We repeated the analysis removing the 16.4 Ma constraint on maximal age of divergence. Second, we ran an independent analysis assuming that the northernmost Andes split ancestral P. pustulosus popu- lations prior to 2.0 Ma. The eastern chain of the northern Andes terminate in a Y-shaped pattern, with the Sierran?a de Perij? along the Colombian border forming the western fork of the Y, and the M?rida Andes of Venezuela forming the eastern fork. Among our samples, Mariquita (MARI) lies to the west of the Andes, Lago Maricaibo (LMAR) in the centre of the ?Y?, and Calabozo (CALA), Carupano (CARU), and Trinidad (TRIN) lie to the east. The northern Andes are thought to have reached their current elevation by around 2.7 Ma (Gregory-Wodzicki 2000), but the Sierra de Perij? and M?rida Andes are probably much older than the Eastern Cordillera (Hoorn et al. 1995). Therefore, our constraint of ? 2.0 Ma is conservative. The COI sequence data in Cannatella et al. (1998) and the COI results in this study (see below) supported P. pustulosus as sister to the other six species of the P. pustulosus group. However, most data sets as well as the combined data set presented in Cannatella et al. (1998) suggested that P. pustulosus is the sister to the P. freibergi + P. petersi clade. Given these two alternative topologies, we chose to assume the latter hypothesis for our divergence time estimation because this topology would again bias our results towards obtaining a younger age for t?ngara frogs. We also assumed that our two island samples, Trinidad (TRIN) and Isla del Rey (IREY), did not share a common ancestor with any mainland sample more recently than 10 000 years ago. Priors were chosen following the guidelines in the multidistribute documentation. Mean of the prior distribution of evolutionary rates was based on a rate calibration of 0.69% uncorrected sequence divergence per lineage per million years (Myr) from a study of Mongolian toads by Macey et al. (1998). To compare our non-clock divergence time estimation with a traditional molecular clock calibration approach, we also applied this rate of 0.69% to our data by averaging all pairwise uncorrected genetic distances across a particular node. Thus we report three wholly independent divergence time estimates: (i) Ecuadorian Andes splitting Pacific and Amazonian species in the mid-Miocene, (ii) northern Andes splitting Colombian and Venezuelan populations 2.0 Ma or more, and (iii) standard application of a ?frog clock?. Results Allozymes Variation in the presumptive gene products of 14 loci is shown in Table 1. There are no significant differences among populations in any of the statistics. Levels of heterozygo- sity are close to those reported by Nevo & Beiles (1991) for seven tropical anurans (mean = 0.069 ? 0.055). Subsequent years sampling in Gamboa (GAMB) resulted in nearly identical allozyme frequencies (pooled), which suggests that the differences we observe among populations do not simply reflect year-to-year variation. The 11 polymorphic loci are presented in Table 2 as allele frequencies and in Appendix I as genotypes. Intergradations at four putative gene loci (6PGD, LDH, AAT, and TPI) is highlighted in Table 2 and is illustrated in Fig. 2 as an ordination of the allozyme frequency data Fig. 2 Multidimensional scaling (MDS) plot of three principal ordinations for allozyme data on Physalaemus pustulosus populations. B I O G E O G R A P H Y O F T ? N G A R A F R O G S 9 ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x by multidimensional scaling (MDS) (Lessa 1990). Multi- dimensional scaling allows populations to be shown as intermediate rather than with one group or another group (Lessa 1990). This is an advantage over a phenetic clustering, which has other shortcomings as well (de Queiroz & Good 1997). There is a mean Nei?s unbiased genetic distance (DN) of 0.294 (range 0.171?0.434) between Middle and South American + Panama populations (Table 3). This distance is 3- to 10-fold greater than the mean DN of 0.026 and 0.083 within these two groups of populations. Previously, a matrix of geographic distances was calcu- lated as the shortest distance between pairs of populations (appendix A in Ryan et al. 1996). Here we present a Mantel test comparing the matrices of pairwise geographic and both pairwise genetic distances (only mtDNA data are new; other data from Appendix A, B, and C in Ryan et al. 1996), along with call dissimilarities matrix between all populations showed a strong correlation for the two measure of genetic distances, allozymes and mtDNA, with each other and with geography. Table 4 shows correlations among all matrix comparisons. Mitochondrial DNA The 564 bp of COI data presented here (GenBank Acces- sion nos DQ120012?DQ120050) correspond to bases 8114?8677 of the Xenopus laevis mitochondrial genome (GenBank Accession no. NC_001573; Roe et al. 1985). All phylogenetic analyses included additional members of the P. pustulosus species group plus three outgroup taxa. In this alignment, 369 characters were constant, 19 sites contained singletons, leaving 176 parsimony-informative characters. Base frequencies across taxa were not heterogeneous ( = 24.82, P = 1.0). On both the NJ and the MP starting trees, the hLRT suggested the most complex, 10-parameter model, GTR + ? + I, which assumes unequal base frequencies, an independent reversible rate for each of the six possible types of nucleotide substitutions, heterogeneity among nucle- otide sites described by a discrete ? distribution with shape parameter, ? (Yang 1994), plus a proportion of invariable sites, I (Hasegawa et al. 1987). The AIC suggested a simpler nine- parameter model (TVM + ? + I) that assumes a single rate for both types of transitions. Values of ML model parameters were similar when estimated from NJ versus parsimony starting trees. The decision theory (DT) technique of model selection, which seeks to minimize error branch length estimation, also recommended the ? and I parameters, but found that only two substitution rates were sufficient in the model, one each for transitions and transversions. This cor- responds to a HKY + ? + I model (Hasegawa et al. 1985). Based on the above results, we conducted six likelihood- based phylogenetic analyses. An ML search with fixed para- meter values was conducted under the TVM + ? + I model Table 3 The upper portion is a group averaging of Nei?s (1978) unbiased genetic distances (DN) from allozyme data. The lower portion is a group averaging of the HKY-IG sequence divergences (*100). ***** indicate no within-area comparisons, as only one population represents the area Mexico to Costa Rica Western Panama Central Panama Eastern Panama Colombia Western Venezuela Eastern Venezuela and Trinidad Mexico to Costa Rica 0.03 0.15 0.27 0.28 0.41 0.33 0.37 3.9 Western Panama 10.8 0.04 0.19 0.19 0.21 0.22 0.25 1.1 Central Panama 13.4 4.3 0.03 0.03 0.12 0.13 0.18 0.6 Eastern Panama 11.1 6.2 7.8 0.02 0.10 0.10 0.17 0.1 Colombia 10.7 5.5 6.4 5.7 ***** 0.20 0.25 Western Venezuela 10.7 5.9 6.6 5.8 3.5 ***** 0.05 Eastern Venezuela and Trinidad 13.4 6.4 6.7 6.6 5.3 5.8 0.05 0.9 ?[ ]114 2 Table 4 Matrix correlation table, Mantel tests using 1000 itera- tions. The first three rows are taken from Ryan et al. (1996). G, standard normal variate; Z, Mantel coefficient; R, correlation coefficient between matrices; P, probability Comparison matrices G Z R P Advertisement call vs. geography 4.87 571 540.52 0.502 0.001 Allozymes vs. advertisement call 5.32 211.75 0.430 0.001 Allozymes vs. geography 7.37 247 214.34 0.732 0.001 mtDNA vs. advertisement call 4.29 93.83 0.353 0.001 mtDNA vs. allozymes 9.78 14.40 0.781 0.001 mtDNA vs. geography 6.34 104 084.00 0.643 0.001 10 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x suggested by the AIC. Full ML and Bayesian analyses were conducted using the GTR + ? + I model. The fixed parameter, full ML, and Bayesian analyses were repeated under the simpler HKY + ? + I model suggested by the DT criterion. The full ML heuristic searches yielded a tree with support ?Ln = 3417.31974 under the GTR + ? + I model, and ?Ln = 3422.30150 under the HKY + ? + I model. Figure 3 presents the ML tree obtained under the HKY + ? + I model, with levels of statistical support obtained by Bayesian analysis under the HKY + ? + I model, parsimony bootstrap analysis, and decay indices. Parsimony analyses recover six most-parsimonious trees of 621 steps. Only one conflicting difference is observed among ML and consensus parsimony topologies, and this involved the root node of the northwestern clade [Mexico?Costa Rica; (ML = GUAT or TEHU + TAPA; parsimony = VERC; Bayesian = trichotomy)]. Our phylogenetic results corroborate the break between northern and southern lineages of P. pustulosus (Ryan et al. 1996) and reveal significant differences among populations in western, central, and eastern Panama. The southern samples of t?ngara frogs can be divided into at least five or six clades. Western and central Panama form a well- supported group, whereas phylogenetic relationships among eastern Panama, Colombia, western Venezuela, and eastern Venezuela + Trinidad are unresolved (Fig. 4). The molecular clock hypothesis was not rejected under either model of sequence evolution tested. Under HKY + ? + I with fixed parameter values, ML tree support, ?Ln = 3422.3187, clock enforced tree ?Ln = 3443.334, for a dif- ference of 21.015 ? 2 = 42.030. Under GTR + ? + I with fixed parameter values, ML tree ?Ln = 3418.028, clock tree ?Ln = 3439.333, for a difference of 21.305 ? 2 = 42.610. With (n ? 2) = 37 degrees of freedom, = 48.363. All optimal trees placed the t?ngara frogs, Physalaemus pustulosus, as sister to the rest of the pustulosus species group (not shown), in agreement with the COI gene tree in Cannatella et al. (1998), although their other data sets placed the t?ngara frogs as sister to the Physalaemus freibergi + Physalaemus petersi sister pair. To test whether our COI topology was significantly better than one constraining the monophyly of (P. freibergi, P. petersi, all P. pustulosus samples), we performed a paired-sites test of Shimodaira & Hasegawa (1999) as implemented in paup*. The SH test used RELL optimization on ML trees obtained with fixed parameter values under both HKY + ? + I and GTR + ? + I models. Despite the ?100% certainty? of the Bayesian anal- ysis, the difference in ?ln L support between optimal and constrained trees was nonsignificant under either model, and this also held true when the 12S data of Cannatella et al. (1998) was also included and re-analysed using dt- modsel or mrbayes. Fig. 3 Phylogeography of Physalaemus pustulosus population samples (indicated by four-letter abbreviations) inferred from ML analysis of COI gene sequences under the HKY + ? + I model. Additional species of the P. pustulosus group and outgroups used in the analyses are not shown. Statis- tical support for branches is indicated by three numbers, separated by slashes. The first number shows the Bayesian marginal posterior probabilities (MPP) multiplied by 100, and assuming an HKY + ? + I model. The middle number reflects the decay index, a.k.a. Bremer support value. The third number shows the percent parsimony bootstrap support (BSS). A dash indicates in the first position an MPP of less than one-half, in the middle position a node not found in the most-parsimonious tree, and in the second position a level of BSS less than 50%. This tree corresponds to node (2) in Fig. 4. ?( . )P=0 1 2 B IO G E O G R A P H Y O F T ? N G A R A F R O G S 11 ? 2005 B lackw ell P u blishing L td , M olecular E cology, 10.1111/ j.1365-294X .2005.02707.x Table 5 Pairwise HKY-IG sequence divergences (? 100) VERC VERC TEHU 6.7 TEHU TAPA 5.9 2.2 TAPA GUAT 4.4 4.5 3.8 GUAT ESAL 4.3 5.4 4.4 3.1 ESAL NICA 3.5 4.9 4.0 3.1 2.2 NICA CRIC 4.1 5.6 4.7 2.9 2.4 0.9 CRIC PARM 11.5 12.4 12.4 9.8 11.6 10.2 10.1 PARM GUAL 10.9 11.7 11.1 9.1 11.0 9.7 9.5 1.1 GUAL SANT 14.1 14.2 13.7 10.9 14.1 12.7 13.2 4.3 3.9 SANT ANTO 13.9 14.1 13.5 10.8 14.7 12.5 13.1 4.4 3.9 1.8 ANTO GATW 13.6 13.8 13.2 10.4 14.3 12.2 12.7 4.2 3.7 1.6 0.2 GATW GATE 13.9 14.1 12.8 10.8 14.7 12.5 13.1 4.4 3.9 1.8 0.4 0.2 GATE BCIC 13.9 14.1 13.5 10.8 14.7 12.5 13.1 4.4 3.9 1.8 0.0 0.2 0.4 BCIC PRPH 14.3 14.5 13.8 11.1 15.0 12.8 13.4 4.6 4.2 2.0 0.2 0.4 0.5 0.2 PRPH GAMB 13.9 14.1 13.5 10.8 14.7 12.5 13.1 4.6 4.2 2.0 0.2 0.4 0.5 0.2 0.4 GAMB GBRG 14.6 14.8 14.2 11.4 15.3 13.2 13.7 4.8 4.4 2.2 0.3 0.5 0.7 0.3 0.5 0.5 GBRG SUMM 13.9 14.8 14.2 10.7 14.6 12.5 13.0 4.8 4.4 2.2 0.3 0.5 0.7 0.3 0.5 0.5 0.7 SUMM CHIV 14.5 14.4 13.7 11.3 15.3 13.1 13.6 4.7 4.2 2.1 0.2 0.4 0.5 0.2 0.4 0.4 0.5 0.5 CHIV COCO 14.3 14.5 13.8 11.1 15.0 12.8 13.4 4.6 4.2 2.0 0.2 0.4 0.5 0.2 0.3 0.4 0.5 0.5 0.4 COCO KOBE 14.5 14.4 13.7 11.3 15.3 13.1 13.6 4.7 4.2 2.1 0.2 0.4 0.5 0.2 0.4 0.4 0.5 0.5 0.0 0.4 KOBE ITAB 13.9 14.1 13.5 10.8 14.7 12.5 13.1 4.4 3.9 1.8 0.0 0.2 0.4 0.0 0.2 0.2 0.3 0.3 0.2 0.2 0.2 ITAB IREY 12.2 11.1 10.5 10.4 11.1 9.7 10.2 6.3 5.8 6.7 7.6 7.3 7.1 7.6 7.9 7.9 8.1 7.6 7.7 7.9 7.7 7.6 IREY METE 12.7 11.6 11.0 10.9 11.5 10.2 10.6 6.6 6.1 7.1 8.0 7.7 7.4 8.0 8.2 8.3 8.5 8.0 8.1 8.2 8.1 8.0 0.2 METE ELRE 12.7 11.6 11.0 10.9 11.5 10.2 10.6 6.6 6.1 7.1 8.0 7.7 7.4 8.0 8.2 8.3 8.5 8.0 8.1 8.2 8.1 8.0 0.2 0.0 ELRE MARI 11.9 12.1 10.1 9.2 12.0 10.0 9.9 6.0 5.1 6.5 6.3 6.0 5.8 6.3 6.5 6.6 6.8 6.8 6.6 6.0 6.6 6.3 5.5 5.8 5.8 MARI LMAR 12.0 11.6 11.0 9.0 11.5 10.1 10.0 6.1 5.6 6.6 6.4 6.2 6.4 6.4 6.7 6.7 6.9 6.9 6.8 6.7 6.8 6.4 5.6 5.9 5.9 3.5 LMAR CALA 15.2 14.7 12.7 11.4 14.0 13.8 13.7 6.0 6.0 6.0 6.0 5.8 5.5 6.0 6.3 6.3 6.0 6.5 6.3 6.3 6.3 6.0 6.5 6.8 6.8 4.8 5.4 CALA CARU 14.2 15.1 13.7 11.0 13.6 13.1 12.6 6.5 7.0 7.0 7.1 6.8 6.5 7.1 7.3 7.3 7.6 7.0 7.4 7.3 7.4 7.1 6.5 6.8 6.8 5.7 6.3 1.4 CARU TRIN 14.9 14.4 13.0 11.0 13.6 13.4 13.3 6.2 6.7 6.7 6.8 6.5 6.3 6.8 7.0 7.0 6.8 7.3 7.1 7.0 7.1 6.8 6.2 6.5 6.5 5.5 5.6 0.5 0.9 TRIN 12 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x To avoid biasing our results towards obtaining a very old t?ngara lineage, we conservatively ran the Bayesian MCMC divergence time analyses on a topology corresponding to the COI ML tree (HKY + ? + I model), but with P. freibergi + P. petersi placed as the sister lineage to P. pustulosus. Among all three Bayesian analyses, the youngest mean divergence times were obtained with the Ecuadorian Andes scenario with both upper and lower temporal constraints applied. These results are shown in Fig. 4. Yet even this analysis placed the 95% confidence interval (CI) of the divergence time of northern vs. southern clades of P. pustulosus nearly 1 Ma prior to the rise of the isthmus (3.1?2.8 Ma, Coates & Obando 1996). Under the northern Andes time constraint, the 95% CI around the basal split within the t?ngara clade preceded the rise of the isthmus by almost 2 Ma (Table 6). Applying a traditional molecular clock also resulted in a P. pustulosus split happening almost 2 Ma prior to the rise of the isthmus. Repeating the analysis using the optimal COI topology with P. pustulosus as sister group to all other pustulosus group species, t?ngara populations would be older still. Discussion Our analyses of COI divergence in Physalaemus pustulosus provide the newest and perhaps strongest molecular evid- ence to date that organisms seemingly intolerant of salt water were able to move between Middle and South America in the late Tertiary well prior to the date of the final comple- tion of the land bridge 3.1 Ma. Previous studies arguing for intercontinental migration prior to the completion of the Panamanian isthmus have relied on external calibrations and fixed-rate molecular clocks (e.g. Zamudio & Greene 1997; Bermingham & Martin 1998; Zeh et al. 2003). In contrast, our study uses not one but two independent temporal Fig. 4 Estimation of absolute divergence time for each node using the Bayesian MCMC method of Kishino et al. (2001) given the ML HKY + ? + I topology and assuming the rise of the Ecuadorian Andes (Fig. 1) between 16.4 and 11.2 million years ago precipitated the split at the ingroup root node (1). Heavy black bars on either side of node (1) illustrate this temporal constraint. Nodes are placed according to their mean divergence time relative to the geological timescale indicated above (Quat., Quaternary). Black nodes received higher statistical support than open nodes (Fig. 4). Error bars on each node denote the central 95% of the estimated posterior probability distribution of divergence time. The long thick vertical grey bar indicates the rise of the Isthmus of Panama between 3.1 and 2.8 million years ago. Node (2) represents the split between northeast and southwest lineages of Physalaemus pustulosus. Note: the estimated divergence time of node (2) is roughly 1 million years older (Table 6) when the Ecuadorian Andes constraint is removed and instead the northern Andes are constrained to be at least 2 million years old (corresponding to the two nodes connecting MARI, LMAR, and CALA). The two heavy vertical black lines next to IREY and TRIN indicate the constraint that these two continental island populations diverged from their mainland relatives at least 10 000 years ago. B I O G E O G R A P H Y O F T ? N G A R A F R O G S 13 ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x constraints internal to the data at hand, and does not rely on the assumption of a strict molecular clock (Table 6). Within P. pustulosus, both nuclear and mitochondrial genetic results define two distinct groups of populations: the northern group that includes populations from Mexico to northern Costa Rica, and the southern group that includes those from western Panama to Trinidad. With the exception of the western Panama populations at this junction, localities that are close geographically are generally similar in allozyme frequencies and have a small DN. Exceptions include Mariquita (MARI), the Colombian population and Lago Maricaibo (LMAR), the western Venezuela population, which are more distinct. Related frogs in Argentina (Pleurodema spp.) showed much less allozyme or call differentiation than observed here (McLister et al. 1991). The DN between P. pustulosus populations in the Middle American group and the southern group (mean = 0.294) are comparable to the value obtained by Cannatella et al. (1998) for Physalaemus petersi vs. Physalaemus freibergi (DN = 0.323). The level of allozyme differentiation between the northern and southern lineages (mean DN = 0.294) is nearly equi- valent to the level of differentiation observed between most closely related of the recognized species in the P. pustulosus species group (lowest DN = 0.323; Cannatella et al. 1998) although these results do not share the same suite of loci, so interpretations should be limited. The mtDNA divergences are more striking. The average sequence divergence between the northern and southern groups is 12.6% and within each group is only 3.9% and 4.5%, respectively. Nuclear allozyme results indicate populations in western Panama are intermediate with bidirectional intergradation (Table 2). The wide geographic distribution of this repre- sentative of an otherwise South American genus without any conspicuous morphological differentiation within Middle America (Freeman 1967) and only partial clinal variation in some advertisement call parameters, would suggest that either P. pustulosus is a recent invader from South America or that there has been substantial gene flow among its populations. The genetic data do not support either of these inferences. How we interpret the informa- tion for these lineages depends on the species concept employed (see de Queiroz 1998). Timing of events Estimates of dates of divergence and rates of evolution can be highly variable, yet we found a remarkable correspond- ence in dates and rates among our three analyses (Table 6). Our two independent rate calibrations offer reciprocal support for one another, and are further corroborated by the molecular clock from Mongolian toads (Macey et al. 1998). In fact, the toads from Tibet appear to have faster rates of mtDNA sequence evolution than three lineages of Neotropical anurans: Physalaemus (Table 6), Eleutherodactylus (Crawford 2003), and Bufo (Mulcahy & Mendelson 2000). Of primary interest is whether P. pustulosus crossed from South to Middle America prior or subsequent to the emergence of the isthmus 2.8?3.1 Ma (Coates & Obando 1996). Our data show convincingly that the northern and southern groups had diverged prior to the oldest hypoth- esized isthmus formation (Fig. 4, Table 6). Both the geolog- ical and palaeontological data suggest that the final closure of the isthmus took place in eastern Panama about 3.1? 2.8 Ma (Whitmore & Stewart 1965; Coates & Obando 1996; Coates et al. 2004). The contact zone between northern and southern t?ngara frogs lies in western Panama. Therefore, we posit that P. pustulosus entered lower Middle America in two or more distinct events. If we average our two COI divergence time estimates (Table 6), the first event involved the basal split of P. pustulosus around 9 Ma (4.5?14.6 Ma) and the arrival of the northern lineage of t?ngara frogs in Middle America. The Panama microplate completed its Table 6 Comparing three independent estimates of rates and dates. The first four columns compare estimates of dates of divergence for four lineage-splitting events. The Sierra de Perij? separates sample localities Mariquita (MARI) and Lago Maracaibo (LMAR), whereas the M?rida Andes separate these two localities from Calabozo (CALA), Carupano (CARU), and Trinidad (TRIN) to the east. The Ecuadorian Andes separate the Pacific coast species from those in the Amazon. Each row corresponds to one analysis, and within each row assumptions and results are shown in bold and plain text, respectively. Rows 1 and 2 show results (with 95% confidence intervals) from ?multidivtime? analyses. Row 1 shows the same results as Fig. 4, plus estimates of the rates of evolution. Row 2 shows the resulting date and rate estimates assuming only that the northern Andes separated Colombian and Venezuelan lineages not less than 2 Myr. Row 3 shows the divergence date estimates (? two standard deviations from the mean) assuming a molecular clock of 0.69%/lineage/Myr Dates in millions of years (Myr) Rates as percentage divergence per lineage per Myr Northern Andes Ecuadorian Andes North vs. South P. pustulosusCalibration Sierra de Perij? M?rida Andes Ecuadorian Andes (F84 + ?) 3.53 (1.03?7.53) 4.83 (1.74?9.31) 11.2?16.4 8.61 (4.07?13.3) 0.49% (0.10?1.3%) Northern Andes (F84 + ?) ? 2.00 ? 2.00 15.5 (10.4?21.8) 9.73 (4.95?15.9) 0.43% (0.085?1.1%) Toad clock (uncorrected percentage) 2.31 (2.31?2.31) 3.47 (2.96?3.98) 12.5 (10.6?14.5) 6.35 (5.04?7.66) 0.69% 14 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x collision with South America by 7.1 Ma (Coates et al. 2004) and the mammalian interchange began at this time with raccoons migrating south and ground sloths moving north 9.3?8.0 Ma (Marshall et al. 1979; Marshall 1985; Webb 1985). Although amphibians are intolerant of salt water, they can on rare occasion disperse across oceans (Hedges 1996; Vences et al. 2003, 2004). One dispersal mechanism for t?ngara frogs could be on a raft launched from the R?o Magdalena in Colombia yet situated north of Panama. The ocean currents along the Caribbean coast of South America then, as now, flow westwards so discharge from the Magdalena flows towards Panama (D?az de Gamero 1996). Our results support the earlier findings of Bermingham & Martin (1998) who were forced to rely on an externally calibrated molecular clock. In their analysis of three genera of freshwater fishes, they conclude that a wave of coloniza- tion from South America into Panama took place 7?4 Ma, with one or two additional invasions following the rise of the isthmus at ?3 Ma. We used two independent geological events to obtain rate and date estimates without assuming a molecular clock and obtained similar results but with even an older dispersal date. The mean mtDNA sequence divergence rate among branches we estimated at 0.46% per lineage per Myr, a rate 33% slower than an often-cited toad rate (0.69% ? Table 6). Because of their smaller size, faster generation time, and perhaps their tropical environment, one might have expected t?ngara frogs to show the faster rate (Adachi et al. 1993; Martin & Palumbi 1993; Rand 1994). This observation suggests that frogs may have similar rates of mtDNA evolution, independent of size. The cane toad, Bufo marinus, and the t?ngara frog are geographically codis- tributed across the t?ngara frog?s range. The largest genetic break among the cane toad populations corresponds to the northern Andes. Slade & Moritz (1998) proposed that the rise of the northern Andes 2.7 Ma (around the area of Bogot?, Colombia) divided east and west populations. Under this model two conclusions are reached. First, the overall rate of DNA sequence evolution at the ND3 gene and tRNA genes is 1% per lineage per Myr in cane toads, more than double the mean rate of COI evolution in P. pustulosus. Second, unlike t?ngara frogs, there is no evidence to sug- gest that cane toads entered Middle America prior to the emergence of the land bridge. A discontinuity between Middle and South American animals in western Panama is not unique to P. pustulosus. Birds (Wetmore 1965; McDonald et al. 2001) and arthropods (Wilcox et al. 1997; Zeh et al. 2003) show similar patterns. A pseudoscorpion has a range similar to P. pustulosus, and exhibits similar COI genetic distances (Wilcox et al. 1997). The pseudoscorpion populations from western and central Panama, separated by only 400 km, are referred to as undergoing ?incipient speciation? owing to reproductive incompatibility and 8.2% sequence divergence and the authors date the separation at 6 Ma. Limited sampling of the bushmaster, genus Lachesis, shows similar patterns (Zamudio & Greene 1997). They recognize the smallest of their Central American diver- gences (5.3%) as distinct evolutionary species. They use a size-specific reptilian clock calculation range of 0.47?1.32% per Myr to indicate that this predates the land bridge emergence with the two groups divergence time of 6.4? 17.9 Ma and the Costa Rican ?species? at 4?11 Ma. What do these divergences and intergradations mean for the frog populations in Panama? Wake & Jockusch (2000) describe problematic salamander taxonomy in the face of multiple data sets. They promote Ghiselin?s (1997; p. 99) species concept: ?Biological species are populations within which there is, but between which there is not, sufficient cohesive capacity to preclude indefinite diver- gence?, and state the critical nature of knowledge of popu- lation interaction and its role in inferring relationships. Under these criteria, the pseudoscorpion data on pre- and postzygotic isolating mechanisms would support their specific designation, whereas the bushmasters? elevation to specific status, lacking such study, would be premature. For P. pustulosus, the intergradation in western Panama populations in the nuclear allozyme data is pivotal. The fact that t?ngara frogs show this cohesion along with little to no ecological, morphological or call difference leads us to conclude that these interactions will contribute to the process of ?de-differentiation? (Wake & Jockusch 2000). Many geological scenarios (Van Andel et al. 1971; Iturralde & MacPhee 1999) indicate connections of various blocks involved in the formation of the rising isthmus around 10 Ma. By 16 Ma, palaeontological data show that North American mammals were well established in central Panama (Whitmore & Stewart 1965). Savage (2002) includes Physalae- mus with a South American element and characterizes the genus as a ?recent contributor to the faunal diversity?. We feel, however, that an earlier movement into Middle America is more likely, with the southern group moving again more recently into Middle America to re-establish contact with the earlier group. These two groups do not presently occupy portions of Costa Rica where suitable habitat seems to exist (Savage 2002), but frogs must have crossed this gap to produce the admixture seen in western Panama. Conclusions The pronounced genetic differences between populations of Physalaemus pustulosus in Middle and South America suggest that they have been separated from one another for 4?13 Myr. A comparison with differences among species in the pustulosus species group separated by the Andes suggests that the isolation between Middle and South American populations of P. pustulosus predates the com- pletion of the Panama land bridge 2.8?3.1 Ma. In spite of the large molecular differences, the allozyme data alone B I O G E O G R A P H Y O F T ? N G A R A F R O G S 15 ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x support the existence of some reproductive compatibility demonstrated by admixture in western Panama. Combined with the absence of documented pre- or post-mating barriers to gene flow, or mate choice preferences, recognition of two species is premature at this time. References Adachi J, Cao Y, Hasegawa M (1993) Tempo and mode of mito- chondrial DNA evolution in vertebrates at the amino acid sequence level: rapid evolution in warm-blooded vertebrates. Journal of Molecular Evolution, 36, 270?281. Akaike H (1974) A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19, 716?723. Bermingham E, Avise JC (1986) Molecular zoogeography of fresh- water fishes in the southeastern United States. Genetics, 113, 939?965. Bermingham E, Martin AP (1998) Comparative mtDNA phylo- geography of Neotropical freshwater fishes: testing shared history to infer the evolutionary landscape of lower Central America. Molecular Ecology, 7, 499?517. Bremer K (1988) The limits of amino acid sequence data in angiosperm phylogenetic reconstruction. Evolution, 42, 795?803. Cannatella DC, Duellman WE (1984) Leptodactylid frogs of the Physalaemus pustulosus group. Copeia, 1984, 902?921. Cannatella DC, Hillis DM, Chippindale et al. (1998) Phylogeny of frogs of the Physalaemus pustulosus species group, with an exam- ination of data incongruence. Systematic Biology, 47, 311?335. Clayton JW, Tretiak DN (1972) Amine citrate buffers in starch gel electrophoresis. Journal of the Fisheries Research Board of Canada, 29, 1169?1172. Coates AG, Collins LS, Aubry M-P, Berggren WA (2004) The geology of the Darien, Panama, and the late Miocene?Pliocene collision of the Panama arc with northwestern South America. Geological Society of America Bulletin, 116, 1327?1344. Coates AG, Jackson JBC, Collins LS et al. (1992) Closure of the Isthmus of Panama, the near-shore marine record of Costa Rica and western Panama. Bulletin of the Geological Society of America, 104, 814?828. Coates AG, Obando JA (1996) The geologic evolution of the Cen- tral American isthmus. In: Evolution and Environment in Tropical America (eds Jackson JBC, Budd AF, Coates AG), pp. 21?56. University of Chicago Press, Chicago, Illinois. Crawford AJ (2003) Relative rates of nucleotide substitution in frogs. Journal of Molecular Evolution, 57, 636?641. D?az de Gamero ML (1996) The changing course of the Orinoco River during the Neogene: a review. Paleogeography, Paleoclima- tology, Paleoecology, 123, 385?402. Duellman WE (1966) The Central American herpetofauna: an ecological perspective. Copeia, 1966, 700?719. Dunn ER (1940) Some aspects of herpetology in lower Central America. Transactions of the New York Academy of Science, 2, 156? 158. Farrell JW, Raffi I, Janecek TR et al. (1995) Late neogene sedimen- tation patterns in the eastern equatorial Pacific Ocean. Proceedings of the Ocean Drilling Program, Scientific Results, 138, 717?756. Felsenstein J (1981) Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution, 17, 368?376. Felsenstein J (1985) Confidence limits on phylogenies: an approach using the bootstrap. Evolution, 39, 783?791. Fitch WM (1971) Toward defining the course of evolution: minimal change for a specific tree topology. Systematic Zoology, 20, 406?416. Freeman HL (1967) Geographic variation in Engystomops pustu- losus (Amphibia: Leptodactylidae) in Middle America. MA Thesis, University of Kansas, Lawrence, Kansas. Frost Darrel R (2004) Amphibian Species of the World: An Online Reference. Version 3.0 (22 August 2004). Electronic Database accessible at http://research.amnh.org/herpetology/amphibia/ index.html. American Museum of Natural History, New York. Garcia-Paris M, Alcobendas M, Buckley D, Wake DB (2003) Dis- persal of viviparity across contact zones in Iberian populations of fire salamanders (Salamandra) inferred from discordance of genetic and morphological traits. Evolution, 57, 129?143. Ghiselin MT (1997) Metaphysics and the Origin of Species. State Uni- versity of New York Press, Albany, New York. Gregory-Wodzicki KM (2000) Uplift history of the Central and Northern Andes: a review. Geological Society of American Bulletin, 112, 1091?1105. Hasegawa M, Kishino H, Yano T (1985) Dating of the human?ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution, 22, 160?174. Hasegawa M, Kishino H, Yano T (1987) Man?s place in Hominoidea as inferred from molecular clocks of DNA. Journal of Molecular Evolution, 26, 132?147. Hedges SB (1996) The origin of West Indian amphibians and reptiles. In: Contributions to West Indian Herpetology: A Tribute to Albert Schwarz (eds Powell R, Henderson RW), Contributions to Herpetology, volume 12, pp. 95?128. Society for the Study of Amphibians and Reptiles, Ithaca, New York. Hoorn C (1993) Marine incursions and the influence of Andean tectonics on the Miocene depositional history of northwestern Amazonia: Results of a palynostratigraphic study. Palaeoclima- tology, Palaeogeography, Palaeoecology, 105, 267?309. Hoorn C, Guerrero J, Sarmiento GA, Lorente MA (1995) Andean tectonics as a cause for changing drainage patterns in Miocene northern South America. Geology, 23, 237?240. Huelsenbeck JP, Ronquist F (2001) mrbayes: Bayesian inference of phylogenetic trees. Bioinformatics, 8, 754?755. International Union of Biochemistry, Nomenclature Committee (1992) Enzyme Nomenclature. Academic Press, San Diego, California. Iturralde MA, MacPhee RDE (1999) Paleogeography of the Carib- bean region: implications for Cenozoic biogeography. Bulletin of the American Museum of Natural History, 238, 1?96. Jukes TH, Cantor CR (1969) Evolution of protein molecules. In: Mammalian Protein Metabolism (ed. Munro HN), pp. 21?123. Academic Press, New York. Kessing B, Croom H, Martin A et al. (1989) The Simple Fool?s Guide to PCR, Version 1.0. University of Hawaii, Honolulu. Kishino H, Thorne JL, Bruno WJ (2001) Performance of a diver- gence time estimation method under a probabilistic model of rate evolution. Molecular Biology and Evolution, 18, 352?361. Lessa E (1990) Multidimensional analysis of geographic genetic structure. Systematic Zoology, 39, 242?252. Lockhart PJ, Steel MA, Hendy MD, Penny D (1994) Recovering evolutionary trees under a more realistic model of sequence evolution. Molecular Biology and Evolution, 11, 605?612. Macey JR, Schulte JAII, Larson A et al. (1998) Phylogenetic rela- tionships of toads in the Bufo bufo species group from the eastern escarpment of the Tibetan Plateau: a case of vicariance and dis- persal. Molecular Phylogenetics and Evolution, 9, 80?87. Marshall LG (1985) Geochronology and land-mammal biochronol- ogy of the transamerican faunal interchange. In: The Great 16 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x American Biotic Interchange (eds Stehli FG, Webb SD), Plenum Press, New York. Marshall LG, Buttler RF, Drake RE, Curtis GH, Telford RH (1979) Calibration of the Great American Interchange. Science, 204, 272?279. Martin AP, Palumbi SR (1993) Body size, metabolic rate, genera- tion time, and the molecular clock. Proceedings of the National Academy of Sciences, USA, 90, 4087?4091. McDonald DB, Clay RP, Brumfield RT, Braun MJ (2001) Sexual selection on plumage and behavior in an avian hybrid zone: experimental tests of male?male interactions. Evolution, 55, 1443?1451. McLister JD, Lougheed SC, Bogart JP (1991) Electrophoretic and vocalization comparisons among three leptodactylid frogs (Pleurodema spp.) from northwestern Argentina. Canadian Journal of Zoology, 69, 2397?2403. Minin V, Abdo Z, Joyce P, Sullivan J (2003) Performance-based selection of likelihood models for phylogeny estimation. Sys- tematic Biology, 52, 674?683. Mulcahy DG, Mendelson JR III (2000) Phylogeography and speci- ation of the morphologically variable, widespread species Bufo valliceps, based on molecular evidence from mtDNA. Molecular Phylogenetics and Evolution, 17, 173?189. Murphy RW, Sites JW Jr, Buth DG, Haufler CH (1990) Proteins I: isozyme electrophoresis. In: Molecular Systematics (eds Hillis DM, Moritz C), pp. 45?126. Sinauer Associates, Sunderland, Massachusetts. Nei M (1978) Estimation of average heterozygosity from a small number of individuals. Genetics, 89, 583?590. Nevo E, Beiles A (1991) Genetic diversity and ecological hetero- geneity in amphibian evolution. Copeia, 1991, 564?592. Posada D, Crandall KA (1998) modeltest: testing the model of DNA substitution. Bioinformatics, 14, 817?818. Pr?hl H, Adams RMM, Mueller U, Rand AS, Ryan MJ (2002) Polymerase chain reaction primers for polymorphic microsatel- lite loci from the t?ngara frog Physalaemus pustulosus. Molecular Ecology Notes, 2, 341?343. de Querioz K (1998) The general lineage concept of species, species criteria, and the process of speciation: a conceptual unification and terminological recommendations. In: Endless Forms: Species and Speciation (eds Howard DJ, Berlocher SH), pp. 57?75. Oxford University Press, New York. de Querioz K, Good DA (1997) Phenetic clustering in biology: a critique. Quarterly Review of Biology, 72, 3?30. Rand DM (1994) Thermal habit, metabolic rate and the evolution of mitochondrial DNA. Trends in Ecology & Evolution, 9, 125?131. Rannala B, Yang Z (1996) Probability distribution of molecular evolutionary trees: a new method of phylogenetic inference. Journal of Molecular Evolution, 43, 304?311. Ridgway GJ, Sherburne SW, Lewis RD (1970) Polymorphism in the esterases of Atlantic herring. Transactions of the American Fisheries Society, 99, 147?151. Roe BA, Ma D-P, Wilson RK, Wong JF-H (1985) The complete nucleotide sequence of the Xenopus laevis mitochondrial DNA genome. Journal of Biological Chemistry, 260, 9759?9774. Rohlf FJ (1993) NTSYS-PC, V. 1.8. Numerical Taxonomy and Multivariate Analysis. Exeter Software, Setauket, New York. Ron SR, Cannatella DC, Coloma LA (2004) Two new species of Physalaemus (Anura: Leptodactylidae) from western Ecuador. Herpetologica, 60, 261?275. Ron RR, Coloma LA, Cannatella DC (2005) A new, cryptic species of Physalaemus (Anura: Leptodactylidae) from Western Ecuador with comments on the call structure of the P. pusulosus group. Herpetologica, 61, 178?198. Ryan MJ (1985) The T?ngara Frog. A Study in Sekual Selection and Communication. University of Chicago Press, Chicago. Ryan MJ, Rand AS (2003) Mate recognition, in t?ngara frogs: a review of some studies of brain, behavior, and evolution. Acta Zoologica Sinica, 49, 721?726. Ryan MJ, Rand AS, Weigt LA (1996) Allozyme and advertisement call variation in the t?ngara frog, Physalaemus pustulosus. Evolu- tion, 50, 2435?2453. Saitou N, Nei M (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution, 4, 406?425. Savage JM (1966) The origins and history of the Central American herpetofauna. Copeia, 1966, 719?766. Savage JM (1982) The enigma of the Central American herpeto- fauna: dispersals or vicariance? Annals of the Missouri Botanical Gardens, 69, 464?547. Savage JM (2002) The Amphibians and Reptiles of Costa Rica. Univer- sity of Chicago Press, Chicago, Illinois. Selander RK, Smith MH, Yang SY, Johnson WE, Gentry JB (1971) Biochemical polymorphism and systematics in the genus Peromyscus. I. Variation in the old-field mouse (Peromyscus polionotus). Studies in Genetics VI University of Texas Publication No. 7103, 49?90. Seutin G, White BN, Boag PT (1991) Preservation of avian blood and tissue samples for DNA analyses. Canadian Journal of Zoology, 69, 82?90. Shimodaira H, Hasegawa M (1999) Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Molecular Biology and Evolution, 16, 1114?1116. Slade RW, Moritz C (1998) Phylogeography of Bufo marinus from its natural and introduced ranges. Proceedings of the Royal Society of London. Series B, Biological Sciences, 265, 769?777. Sorenson MD (1999) TREEROT, Version 2. Boston University, Boston, Massachusetts. Stehli FG, Webb SD, eds. (1985) The Great American Biotic Inter- change. Plenum Press, New York. Swofford DL (1998) PAUP*4 0b Phylogenetic Analysis Using Parsimony, Version 4. Sinauer Associates, Sunderland, Massachusetts. Swofford DL, Berlocher SH (1987) Inferring evolutionary trees from gene frequency data under the principle of maximum par- simony. Systematic Zoology, 36, 293?325. Swofford DL, Selander RB (1981) biosys-i: a FORTRAN program for the comprehensive analysis of electrophoretic data in population genetics and systematics. Journal of Heredity, 72, 281? 283. Thorne JL, Kishino H (2002) Divergence time and evolutionary rate estimation with multilocus data. Systematic Biology, 51, 689? 702. Thorne JL, Kishino H, Painter IS (1998) Estimating the rate of evolution of the rate of molecular evolution. Molecular Biology and Evolution, 15, 1647?1657. Van Andel TH, Heath GR, Malfait BT, Hendricks DF, Ewing JI (1971) Tectonics of the Panama Basin, eastern Pacific. Geological Society of America Bulletin, 82, 1489?1508. Vanzolini PE, Heyer WR (1985) The American hepetofauna and the interchange. In: The Great American Biotic Interchange (eds Stehli FG, Webb SD), pp. 475?487. Plenum Press, New York. Vences M, Vieites DR, Glaw F et al. (2003) Multiple overseas dis- persal in amphibians. Proceedings of the Royal Society of London. Series B, Biological Sciences, 270, 2435?2442. B I O G E O G R A P H Y O F T ? N G A R A F R O G S 17 ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x Vences M, Kosuch J, R?del M-O et al. (2004) Phylogeography of Ptychadena mascareniensis suggests transoceanic dispersal in a widespread African?Malagasy frog lineage. Journal of Biogeo- graphy, 31, 593?601. Wake DB, Jockusch EL (2000) Detecting species borders using diverse data sets: examples from plethodontid salamanders in California. In: The Biology of Plethodontid Salamanders (eds Bruce RC, Jaeger RG, Houck LD), pp. 95?119. Kluwer Academic/ Plenum Publishers, New York. Webb SD (1985) Main pathways of mammalian diversification in North America. In: The Great American Biotic Interchange (eds Stehli FG, Webb SD), Plenum Press, New York. Wetmore A (1965) The birds of the Republic of Panama. Part 1. Tinamidae (Tinamous) to Rynchopidae (Skimmers). Smithsonian Miscellaneous Collections 150, 1?105. Whitmore FC Jr, Stewart RH (1965) Miocene mammals and Cen- tral American seaways. Science, 148, 180?185. Wilcox TP, Hugg L, Zeh JA, Zeh DW (1997) Mitochondrial DNA sequencing reveals extreme genetic differentiation in a cryptic species complex of Neotropical pseudoscorpions. Molecular Phylogenetics and Evolution, 7, 208?216. Yang Z (1994) Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. Journal of Molecular Evolution, 39, 306?324. Yang Z, Rannala B (1997) Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method. Molecu- lar Biology and Evolution, 14, 717?724. Zamudio KR, Greene HW (1997) Phylogeography of the bush- master (Lachesis muta: Viperidae): imiplications for neotropical biogeography, systematics, and conservation. Biological Journal of the Linnean Society, 62, 421?442. Zeh DW, Zeh JA, Bonilla MM (2003) Phylogeography of the giant harlequin beetle (Acrocinus longimanus). Journal of Biogeography, 30, 747?754. Lee A. Weigt pursues molecular systematics, cryptic species and biodiversity while running molecular facilities. Andrew J. Craw- ford studies evolutionary genetics and biogeography of Neotrop- ical frogs, while A. Stanley Rand has long studied tropical biology, herpetology and behavioural ecology. Michael J. Ryan?s research interests are the evolution of behaviour, animal communication and sexual selection. 18 L . A . W E IG T E T A L . ? 2005 B lackw ell P u blishing L td , M olecular E cology, 10.1111/ j.1365-294X .2005.02707.x Appendix I Genotypes by population with electrophoretic conditions noted for each gene product (Buffers: RW, Ridgway et al. 1970; CT, Clayton & Tretiak 1972; TC , Selander et al. 1971; Enzyme numbers from IUCBN 1992) Gene product Genotype Populations VERC TEHU TAPA GUAT ESAL NICA CRIC PARM GUAL SANT ANTO GATW GATE BCIC PRPH GAMB GBRG SUMM CHIV COCO KOBE ITAB IREY METE ELRE MARI LMAR CALA CARV TRIN Phosphogluconate dehyrogenase 6PGDH AB 1 1.1.1.44 BC 4 1 buffer = CT CC 5 7 5 4 7 9 18 8 3 2 5 4 15 3 1 5 CD 7 CE 1 4 6 5 5 6 7 4 9 7 9 7 1 5 5 CG 2 1 DD 1 DG 1 DH 1 EE 3 1 5 1 1 4 1 3 1 2 14 2 1 1 1 EF 1 1 3 EG 4 1 FF 11 14 1 FG 2 GG 11 7 15 8 7 7 1 13 13 GH 1 2 5 3 6 HH 4 Aspartate transaminase ATA AA 11 11 16 10 9 12 11 10 8 1 1 2 1 2.6.1.1 AB 1 buffer = RW AC 3 3 5 8 9 5 7 1 7 4 2 7 3 5 1 1 3 1 AD 1 BC 2 BD 2 CC 1 5 5 8 8 14 14 23 11 11 8 11 7 15 15 14 14 12 1 8 5 CD 1 1 2 1 DD 6 3 Triose phosphate isomerase TPI AA 4 5.3.1.1 AD 3 1 buffer = CT BB 4 1 BH 8 7 2 CH 7 1 DD 4 9 11 9 12 12 11 2 DF 1 DG 1 DH 1 4 5 7 1 2 EH 2 FF 2 9 3 5 4 4 6 1 1 1 FH 1 5 5 9 7 12 8 6 7 5 8 3 4 HH 1 10 5 7 13 8 1 5 7 15 3 3 2 2 5 11 10 7 4 15 2 9 13 3 HI 2 Lactate dehydrogenase LDH1 AB 2 1.1.1.27 BB 2 12 10 8 12 15 15 14 14 27 16 13 15 14 14 15 15 15 15 15 2 12 14 5 buffer = TC BC 5 4 7 1 3 4 5 1 1 3 1 CC 1 5 12 1 12 11 11 2 CD 7 1 DD 3 LDH2 AC 3 1 BC 1 1 1 CC 11 8 7 10 12 12 10 15 14 16 15 14 15 14 15 30 16 13 14 15 14 15 15 15 15 15 2 12 14 5 CD 6 DD 2 B IO G E O G R A P H Y O F T ? N G A R A F R O G S 19 ? 2005 B lackw ell P u blishing L td , M olecular E cology, 10.1111/ j.1365-294X .2005.02707.x Phosphoglucomutase PGM AD 2 2.7.5.1 BD 1 2 1 1 buffer = TC CD 1 DD 11 11 16 8 10 10 10 12 13 16 15 15 15 14 15 30 16 12 14 15 14 15 15 15 15 15 2 DE 1 1 3 1 1 6 EE 1 11 6 5 Glucose phosphate isomerase GPI AB 1 5.3.1.9 AC 1 buffer = CT BB 3 BC 10 3 2 CC 11 11 16 10 12 12 11 1 3 16 15 15 15 14 15 29 16 13 12 15 14 15 14 15 10 1 2 10 14 1 CD 3 CE 1 1 4 1 CF 8 DE 1 EE 1 FF 14 1 FG 1 Malate dehyrogenase MDH1 AB 1 1.1.1.42 BC 1 6 buffer = RW BD 1 CC 11 11 16 10 12 12 11 15 14 16 15 2 5 12 12 29 15 9 12 14 12 13 14 14 13 15 2 12 7 5 CD 11 6 2 3 1 1 4 2 1 2 1 1 2 DD 1 4 MDH2 AC 7 2 BB 2 BC 10 5 5 1 BE 8 CC 11 11 16 10 12 12 11 15 4 12 8 15 8 14 13 28 16 11 14 10 14 14 14 13 15 2 12 13 5 CD 1 CE 2 2 2 2 1 2 CF 1 EE 5 Isocitrate dehydrogenase ICD1 AB 1 1.1.1.42 BB 3 9 14 1 5 7 5 6 4 5 12 10 11 1 2 2 buffer = TC BC 1 BD 5 5 1 4 7 3 8 14 4 1 6 8 4 7 20 8 5 2 5 1 9 4 15 13 11 10 BE 7 CE 2 DD 5 3 1 6 5 9 3 1 8 7 2 2 3 3 3 2 1 5 9 2 1 4 ICD2 AB 1 1 BB 11 11 16 10 12 12 10 15 13 16 15 15 15 14 15 29 16 13 15 15 14 15 15 15 15 15 2 12 14 2 ME, SOD1, AA 11 11 16 10 12 12 11 15 14 16 15 15 15 14 15 30 16 13 15 15 14 15 15 15 15 15 2 12 13 SOD2 ME, Malate dehydrogenase (+ NADP), 1.1.1.40, buffer = RW. SOD, Superoxide dismutase, 1.15.1.1, buffer = TC. Gene product Genotype Populations VERC TEHU TAPA GUAT ESAL NICA CRIC PARM GUAL SANT ANTO GATW GATE BCIC PRPH GAMB GBRG SUMM CHIV COCO KOBE ITAB IREY METE ELRE MARI LMAR CALA CARV TRIN Appendix I Continued 20 L . A . W E I G T E T A L . ? 2005 Blackwell Publishing Ltd, Molecular Ecology, 10.1111/j.1365-294X.2005.02707.x Appendix II PCR thermal cycler profiles and COI primers used in this study PCR amplifications were 50 ?L total volume with conserved or ?universal? primers for cytochrome oxidase I (COI) COIf (5?-CCTGCAGGAGGAGGAGAYCC-3?) and COIa (5?- AGTATAAGCGTCTGGGTAGTC-3?) (Kessing et al. 1989) or COIa2 (5?-CCTGCYARYCCTARRAARTGTTGAGG-3?), and Taq or Tfl at 1.25 units/50 ?L volume. The thermal cycler profile was: 1.5 min at 94 ?C; 5 cycles: 30 s at 94 ?C, 30 s at 45 ?C, 45 s at 72 ?C; 25 cycles: 30 s at 94 ?C, 30 s at 52 ?C, 45 s at 72 ?C. After thermal cycling, 5 ?L was run on a 1.5% agarose gel and bands were excised and resuspended in 400 ?L, heated to 65 ?C and 1 ?L was used as template for a second round of PCR: 1.5 min at 94 ?C; 25 cycles: 30 s at 94 ?C, 30 s at 52 ?C, 45 s at 72 ?C. PCR products were purified via GeneClean (Bio 101) and cycle sequenced (50 ?C annealing temperature) with ABI dye terminator chemistry following the manufacturer?s instructions. Additional internal sequencing primers included: PP6 (5?-TCTGCAACAATAATYATYCGCAATTCCAAC-3?), PP7 (5?-GTTGGA ATTGCRATGATTATTGTTGCAGA-3?), PP8 (5?-TCTCTAGAYATTGTATTACATGA-3?), and PP9 (5?-TCATGTAATACAATGTCTAGAGA-3?)