Historical biogeography of lowland species of toads (Bufo) - UTA [PDF]

Journal of Biogeography (J. Biogeogr.) (2006) 33, 1889–1904

ORIGINAL ARTICLE

Historical biogeography of lowland species of toads (Bufo) across the TransMexican Neovolcanic Belt and the Isthmus of Tehuantepec Daniel G. Mulcahy*, Benson H. Morrill and Joseph R. Mendelson III 

Department of Biology, Utah State University, Logan, UT, USA

ABSTRACT

Aim In this study, we investigate phylogeographic structure in two different species groups of lowland toads. First, we further investigate strict parapatry of the Pliocene-vicariant Bufo valliceps/B. nebulifer species pair. Secondly, we test for similar phylogeographic structure in the distantly related toad B. marinus, a species we hypothesize will show a Pleistocene dispersal across the same area. Location The eastern extension of the Trans-Mexican Neovolcanic Belt (TMNB) contacts the Atlantic Coast in central Veracruz, Mexico. Although it is not a massive structure at this eastern terminus, the TMNB has nonetheless effected vicariance and subsequent speciation in several groups of animals. The Isthmus of Tehuantepec unites the North American continent with Nuclear Central America and is also known to be a biogeographic barrier for many taxa. Methods We use sequence data from two mitochondrial DNA genes (c. 550 base-pairs (bp) of 16S and c. 420 bp of cyt b) from 58 individuals of the B. valliceps/nebulifer complex, collected from 24 localities. We also present homologous sequence data from 23 individuals of B. marinus, collected from 12 localities. We conduct maximum-parsimony, maximum-likelihood and Bayesian analyses to investigate phylogeographic structure. We then use parsimony- and likelihood-based topology tests to assess alternative phylogenetic hypotheses and use a previously calibrated molecular rate of evolution to estimate dates of divergence. Results Our results further define the parapatric contact zone across the TMNB between the Pliocene-vicariant sister species B. valliceps and B. nebulifer. In contrast, phylogenetic structure among populations of B. marinus across the TMNB is much shallower, suggesting a more recent Pleistocene dispersal in this species. In addition, we found phylogeographic structure associated with the Isthmus of Tehuantepec in both species groups.

*Correspondence: Daniel G. Mulcahy, Department of Herpetology, California Academy of Sciences, 875 Howard Street, San Francisco, CA 94103-3009, USA. E-mail: [email protected]

Main conclusions The existence of a Pliocene–Pleistocene seaway across the Isthmus of Tehuantepec has been controversial. Our data depict clades on either side of the isthmus within two distinct species (B. valliceps and B. marinus), although none of the clades associated with the isthmus, for either species, are reciprocally monophyletic. In the B. valliceps/B. nebulifer complex, the TMNB separation appears to predate the isthmian break, whereas in B. marinus dispersal across the TMNB has occurred subsequent to the presence of a barrier at the Isthmus of Tehuantepec.

 Present address: Department of Herpetology, Zoo Atlanta, 800 Cherokee Ave SE, Atlanta, GA 30315-1440, USA.

Keywords Bufo nebulifer, Bufo marinus, Bufo valliceps, Bufonids, Central America, Mexico, mitochondrial DNA, phylogeography.

ª 2006 The Authors Journal compilation ª 2006 Blackwell Publishing Ltd

www.blackwellpublishing.com/jbi doi:10.1111/j.1365-2699.2006.01546.x

1889

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III INTRODUCTION The Trans-Mexico Neovolcanic Belt (TMNB) is one of the predominant geographical features of Mexico, and its geological development has been posited as a primary contributor to the biogeographic histories of many upland taxa in central Mexico (e.g. Campbell & Frost, 1993; Darda, 1994; Sullivan et al., 2000; Castoe et al., 2003). However, a recent series of papers have independently demonstrated the considerable influence of this transverse massif on the biogeography and evolution of the lowland fauna on both the Pacific (Mateos, 2005) and the Atlantic coasts of Mexico (Mulcahy & Mendelson, 2000; Hulsey et al., 2004; Zaldı´var-Rivero´n et al., 2004). These papers generally support earlier hypotheses based on observations of consistent disjunctions in the distribution of fishes (Rosen, 1978), and of reptiles and mammals (Pe´rezHigareda & Navarro, 1980), specifically along the Atlantic Versant of Mexico. Therefore, it is suggested that the TMNB may form a common geographical barrier to lowland species in this region. The concept of a massive volcanic chain such as the TMNB acting as a vicariant feature to lowland populations is easily tractable, but the reality is that the TMNB withers to a tiny string of lava-rock strewn hills at its eastern terminus (Fig. 1). It makes final contact with the current coastline, as a series of small fingers of raised lava-rock (the northern-most being just south of the small town of Palma Sola, the southern-most just north of the city Cardel, with a small pocket of suitable habitat occurring in the middle, near the town of El Viejo´n) all in the

state of Veracruz, Mexico (Fig. 2). The imposing backbone of the TMNB began activity during the mid-Miocene to Pliocene (de Cserna, 1989), with its greatest development during the Pliocene in the eastern portion (Ferrusquia-Villafranca, 1993), giving rise to the highest peaks of Mexico. However, the easternmost fingers of the TMNB, where they contact the Gulf of Mexico, are barely noticeable to the casual observer. During much of the Pliocene, sea-level maxima caused by a warmer climate effectively covered the entire coastal plain of eastern Mexico (Bryant et al., 1991). Later during the Pleistocene, sea levels fluctuated with corresponding glacial and inter-glacial cycles (Beard et al., 1982). Glacial-maxima lowered sea levels exposing large sections of the continental shelf, which greatly increased the aerial extent of the coastal plain (Fig. 1), whereas glacial minima raised sea levels, inundating much of the coastal plain by about 300 m (Ewing & Lopez, 1991). The area of the Sierra de Los Tuxtlas (Fig. 1) is of volcanic origin from the late Cenozoic (Ferrusquia-Villafranca, 1993), and may have been isolated from the mainland during the Pliocene and during Pleistocene glacial minima. Farther south, debate continues as to whether a seaway existed at the Isthmus of Tehuantepec, separating northern Mexico from Nuclear Central America (Campbell, 1999). By comparing phylogenetic structure from mtDNA sequences of lowland toads in this region, Mulcahy & Mendelson (2000) demonstrated that the prevailing concept of a single wide-ranging species (Bufo valliceps Wiegmann) should be replaced by recognition of a species pair showing an apparent parapatric distribution in the region: B. valliceps

Figure 1 Map of the study area, showing major geographic features discussed in the text. Grey shading from light to dark indicates elevations from 300–900 m, 900–2100 m and > 2100 (including black) m a.s.l., respectively. Dotted line shows extent of continental shelf, much of which was exposed during lower sea levels of the Pleistocene (from Bryant & Bryant, 1991).

1890

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo)

Figure 2 Map of the eastern terminus of the Trans-Mexico Neovolcanic Belt (TMNB) in central Veracruz. Squares indicate reference towns used in text, circles indicate collecting localities for Bufo nebulifer (Mx8) B. valliceps (Mx9–10) and B. marinus (Mx8– 10) in relation to the eastern-most portion of the TMNB (see Figs 3 & 4 for complete sampling of each taxon). Shaded contours indicate elevation following Fig. 1. Los Tuxtlas is the Pliocene volcanic uplift on the coastal plain south of the eastern terminus of the TMNB (see Discussion).

ranging from central Veracruz, Mexico, southward to Costa Rica; and B. nebulifer Girard, ranging from central Veracruz northward to the southern United States of America. Mulcahy & Mendelson (2000) proposed and tested two historical hypotheses related to the timing of this speciation event: (1) Miocene–Pliocene vicariance associated with the orogeny of the TMNB; and (2) Pleistocene dispersal and vicariance caused by low sea levels initially exposing the continental margin, followed by raised sea levels, which obliterated the coastal plain in this region. Their results supported the Miocene–Pliocene vicariance hypothesis. However, their sampling was insufficient to demonstrate strict parapatry of these two species on the narrow coastal plains on either side of the TMNB. A third species, B. marinus Linnaeus, ranges from the Amazon Basin of South America to the southern tip of Texas, USA, in North America. This species is largely sympatric with B. valliceps, partially sympatric with B. nebulifer, and has an apparently continuous distribution across the eastern terminus of the TMNB. Bufo marinus is of a South American origin (Pauly et al., 2004; Pramuk, 2006), and most likely entered Central America during the Pliocene (Slade & Moritz, 1998). A broad-scale phylogeographic study of B. marinus (Slade & Moritz, 1998) indicated that this complex showed dramatic historical effects of the orogeny of the Andes in South America, and lower levels of genetic divergence between samples from Costa Rica and Mexico. These results suggest that a finer-scale study may show that the biogeographic history of northern B. marinus is more consistent with the Pleistocene dispersal and vicariance hypothesis of Mulcahy & Mendelson (2000).

The three species of toads in our study are ecologically similar in their general reproductive biology and overall natural history. All three species are invasive, ‘weedy’ species that are typically more abundant in secondary, degraded habitats than in undisturbed primary forests (Zug & Zug, 1979; Mendelson, 1994; Lee, 1996; Campbell, 1998; McCranie & Wilson, 2002; Savage, 2002). These attributes would suggest that they are suitably comparable to one another, in order to test hypotheses of historical biogeography and that their invasive, dispersal-prone tendencies would make them a conservative test of the effect of the TMNB on lowland species. In this paper, we use c. 970 base-pairs (bp) of sequence data from two mtDNA genes (cyt b and 16S) from recently collected samples along a geographic transect across the TMNB to address the historical biogeography of the Atlantic Versant lowlands of Mexico. Specifically, we test three main hypotheses regarding the relationships and distributions concerning bufonid toads in this region: (1) the hypothesis that B. valliceps and B. nebulifer are parapatric at the eastern terminus of the TMNB (i.e. as proposed by Mulcahy & Mendelson, 2000); (2) the hypothesis that the sympatric toad B. marinus shows a pattern of more recent dispersal across the TMNB, consistent with the Pleistocene dispersal, followed by the vicariance hypothesis of Mulcahy & Mendelson (2000); and (3) the hypothesis that there is evidence of a phylogeographic signal within B. valliceps and/or B. marinus, which is consistent with a historical seaway across the Isthmus of Tehuantepec. We use both parsimony and likelihood topology tests to compare our phylogenetic results with alternative hypotheses. We also discuss genetic variation, in terms of sequence divergence and the rate of molecular evolution previously calibrated for the B. valliceps/nebulifer group (Mulcahy & Mendelson, 2000), between clades associated with either sides of the TMNB and the Isthmus of Tehuantepec. MATERIALS AND METHODS Taxon sampling We collected sequence data for 19 new B. valliceps and B. nebulifer specimens along the eastern terminus of the TMNB (Fig. 3). For consistency, we continued our numbering system from that of Mulcahy & Mendelson (2000). Ten of the new specimens represent B. nebulifer (locality Mx8, Figs 2 & 3) and were collected just north of the town of Palma Sola, seven represent B. valliceps (locality Mx10, Figs 2 & 3) collected south of Cardel, and two were collected in the middle of the eastern-most reaches of the TMNB, near the town of El Viejo´n, tentatively assigned to B. valliceps (locality Mx9, Figs 2 & 3). We compared our new sequences with those of Mulcahy & Mendelson (2000), and used one representative of each unique haplotype in the phylogenetic analyses. This resulted in a total of 28 specimens of B. nebulifer (11 haplotypes) and 30 specimens of B. valliceps (21 haplotypes; see Results). Multiple specimens from the same locality were individually labelled

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1891

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III

Figure 3 Map of the geographic distribution and sampling for the Bufo valliceps and B. nebulifer used in this study. Squares indicate new sample localities, circles represent those from Mulcahy & Mendelson (2000); grey circles indicate identical haplotypes omitted from phylogenetic analyses. Box indicates area of the TMNB transect highlighted in Fig. 2.

alphabetically (e.g. Mx8a–Mx8j). We also collected sequence data from 23 specimens of B. marinus from Central America and northern Mexico (Fig. 4). Again, for sake of convenience, we kept the same numbering system as the B. valliceps/ B. nebulifer data because many of our sampling localities were the same (e.g. Mx8–Mx10). Outgroup taxa were chosen following phylogenetic hypotheses of Mulcahy & Mendelson (2000; D.G. Mulcahy & J.R. Mendelson, unpubl. data); Pauly et al. (2004) and Pramuk (2006). Outgroup taxa for analyses of B. nebulifer and B. valliceps included B. mazatlanensis, B. campbelli (AY008253–4 of Mulcahy & Mendelson, 2000), and B. coccifer (AY927856 and AY927863 of Mendelson et al., 2005). Outgroup taxa for analyses of B. marinus included B. crucifer, B. schneideri and a specimen of B. marinus from South America (Pramuk, 2006). Slade & Moritz (1998) showed that B. marinus may be paraphyletic with respect to South American samples and B. schneideri, therefore we considered the South American individual of B. marnius as an outgroup. Voucher specimens for this study are deposited at the following institutions: University of Texas at Arlington (UTA); Museo de Zoologia, Facultad de Ciencias, Universidad Autonoma de Mexico (MZFC), and University of Kansas (KU). Voucher numbers, specific locality information and GenBank accession numbers for each specimen are listed in Table 1. 1892

Laboratory protocols We collected sequence data from c. 550 bp region of the ribosomal gene 16S and c. 420 bp region of the protein-coding gene cyt b. The gene-regions have been demonstrated to be accurate markers for resolving recent divergences (Graybeal, 1993, 1997; Mulcahy & Mendelson, 2000). For all specimens, tissue was taken from liver, muscle, skin or toe clips and stored either at )80 C or in 95% ethanol. Extraction of DNA followed standard phenol/chloroform extraction methods (Maniatis et al., 1982). Primers used for PCR amplification and sequence reactions (for cyt b, MVZ43: GAGTCTGCCT[A/ T]AT[T/C]GC[C/T]CA[A/G]AT, 3¢ and MVZ28: CGAGGC[G/C]CC[T/C]GCAAT[A/G]ATAA, 3¢; for 16S, 16Sar: CGCCTGTTTATCAAAAACAT, 3¢ and 16Sbr: CCGGTCTGAACTCAGATCACGT 3¢). Amplifications for PCR followed that of Mulcahy & Mendelson (2000). Sequence reactions were done in both directions using the same primers used for PCR, with BigDyeTM Terminator Cycle Sequencing Kits (Version 2, ABI Part No. 4303152; Foster City, CA, USA) in 10–12 mL reactions following recommended protocols. Sequence-reaction products were cleaned with SephadexTM (Piscataway, NJ, USA) and run on an ABI 377 automated sequencer. Individual sequences were aligned with complimentary strands in SequencherTM 4.1.2 (Gene Codes Corp.,

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo)

Figure 4 Map of the northern distribution and sampling for Bufo marinus (squares) in this study. Sample localities follow the same numbering system from the B. valliceps and B. nebulifer, with the addition of several new localities for B. marinus only. Grey squares indicate identical haplotypes omitted from phylogenetic analyses. Box indicates area of the TMNB transect highlighted in Fig. 2.

Ann Arbor, MI, USA). The protein-coding region of cyt b was translated into amino acid sequence and the G-C content was inspected to verify that authentic mtDNA sequences were obtained, using MacClade 4.0 (Maddison & Maddison, 2000). Phylogenetic analysis To determine if the two gene regions, 16S and cyt b, contained similar phylogenetic signals, we conducted partition-homogeneity tests. We ran 100 replicates with 10 random additions at each replicate between 16S and cyt b. This test determines whether the data sets are giving significantly different signals. Maximum parsimony (MP) and Maximum-likelihood (ML) analyses were conducted in paup* beta version 4.0b10 (Swofford, 1999). All MP analyses were conducted with Tree-bisection reconnection branchswapping, ACCTRAN optimization, with 100 random stepwise additions. Non-parametric bootstrap analyses under the MP criterion were conducted in paup* with 1000 replicates, using random step-wise additions with 25 replicates for each bootstrap replicate. Nodes with indices of 75% or greater were considered well supported. The program Modeltest 3.06 (Posada & Crandall, 1998) was used to select appropriate models of evolution, and was conducted on each gene separately, as well as both genes combined for each data set (B. valliceps/B. nebulifer and B. marinus). The combined

genes model parameters were used in paup* for the ML analyses of both genes combined for 100 step-wise random addition replicates. Three Bayesian Inference (BI) analyses were conducted on combined (gene) data sets for each group using the program Mr. Bayes version 10.3 (Huelsenbeck & Ronquist, 2001). Each gene-region was run under its own model using data partitions. Each analysis was run for 5 · 106 generations, each using four-heated Markov chains (using program defaults) and each was sampled every 100 generations. Posterior-probabilities were estimated by conducting a 50% majority-rules consensus tree after discarding the burn-in trees. Runs were determined to reach stationarity visually by plotting likelihood scores against generations. Nodes with 95% or greater posterior-probabilities were considered significant. We use the rate of molecular evolution of 0.33% per lineage, per million years, (Mulcahy & Mendelson, 2000) to estimate dates of divergence in the B. valliceps/nebulifer complex. However, we could not confidently calibrate a rate in the B. marinus group because the B. valliceps/nebulifer and B. marinus data are not clock-like when compared with other Middle American bufonids (D.G. Mulcahy & J.R. Mendelson, unpubl. data). Therefore, we report average pair-wise sequence divergence between major clades in both groups to estimate the relative age of divergences, but acknowledge those in the B. marinus group may be less accurate.

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1893

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III Table 1 Specimen information for sequence data collected for this study (other specimens used in this study are from Mulcahy & Mendelson, 2000). The first column list the species and locality number from maps (Figs 2–4), and names in phylogenetic trees (Figs 5 & 6). The second column gives general locality followed by museum specimen number and GenBank accession numbers for 16S and cyt b Taxon ID

Locality

Bufo nebulifer Mx8a Mx8b Mx8c Mx8d Mx8e Mx8f Mx8g Mx8h Mx8i Mx8j

Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico:

Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz:

north north north north north north north north north north

B. valliceps Mx9a Mx9b Mx10a Mx10b Mx10c Mx10d Mx10e Mx10f Mx10g

Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico: Mexico:

Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz: Veracruz:

near El Viejo´n near El Viejo´n south of Cardel south of Cardel south of Cardel south of Cardel south of Cardel south of Cardel south of Cardel

of of of of of of of of of of

Palma Palma Palma Palma Palma Palma Palma Palma Palma Palma

Sola Sola Sola Sola Sola Sola Sola Sola Sola Sola

Voucher No.

16S

Cyt b

UTA A-54860 UTA A-54861 UTA A-54862 UTA A-54863 UNAM-JRM 4851 UNAM-JRM 4852 UTA A-54864 UTA A-54865 UNAM-JRM 4855 UTA A-54859

DQ415599 DQ415600 DQ415601 DQ415602 DQ415603 DQ415604 DQ415605 DQ415606 DQ415607 DQ415608

DQ415619 DQ415620 DQ415621 DQ415622 DQ415623 DQ415624 DQ415625 DQ415626 DQ415627 DQ415628

UNAM-JRM 4824 UNAM-JRM 4825 UNAM-JRM 4799 UNAM-JRM 4801 UTA A-54855 UTA A-54856 UTA A-54858 UTA A-54854 UTA A-54857

DQ415609 DQ415610 DQ415611 DQ415612 DQ415613 DQ415614 DQ415615 DQ415616 DQ415617

DQ415629 DQ415630 DQ415631 DQ415632 DQ415633 DQ415634 DQ415635 DQ415636 DQ415637

B. mazatlanensis

Mexico: Sinaloa: near Cosala

UNAM-JRM 4491

DQ415618

DQ415638

B. marinus Mx8a Mx8b Mx8c Mx8d Mx9a Mx9b Mx9c Mx9d Mx9e Mx9f Mx9g Mx10a Mx10b Mx11a Mx11b Mx12 C1 E1 E2 G1 G5 H3 H4

Mexico: Veracruz: north of Palma Sola Mexico: Veracruz: north of Palma Sola Mexico: Veracruz: north of Palma Sola Mexico: Veracruz: north of Palma Sola Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: near El Viejo´n Mexico: Veracruz: south of Cardel Mexico: Veracruz: south of Cardel Mexico: Guerrero: near Atoyac Mexico: Guerrero: near Atoyac Mexico: Sinaloa: near Cosala Costa Rica: Heredia: at Chilamate El Salvador: Ahuachapan: El Imposible El Salvador: Ahuachapan: El Imposible Guatemala: Huehuetenango: near Nenton Guatemala: Izabal: Montanas del Mico Honduras: El Paraiso: Las Manos Honduras: Colon: Quebrada Machin

UTA A-54875 UNAM-JRM 4846 UNAM-JRM 4848 UNAM-JRM 4844 UTA A-54882 UTA A-54873 UNAM-JRM 4834 UNAM-JRM 4835 UTA A-54881 UTA A-54879 UTA A-54877 UTA A-54878 UTA A-54871 UTA A-54869 UTA A-54870 UTA A-54868 toe-clip only KU 289750 KU 289772 UTA A-50876 UTA A-50870 UTA A-50638 USNM 534124

DQ415547 DQ415548 DQ415549 DQ415550 DQ415551 DQ415552 DQ415553 DQ415554 DQ415556 DQ415555 DQ415557 DQ415558 DQ415559 DQ415560 DQ415561 DQ415562 DQ415563 DQ415564 DQ415565 DQ415566 DQ415567 DQ415568 DQ415569

DQ415573 DQ415574 DQ415575 DQ415576 DQ415577 DQ415578 DQ415579 DQ415580 DQ415582 DQ415581 DQ415583 DQ415584 DQ415585 DQ415586 DQ415587 DQ415588 DQ415589 DQ415590 DQ415591 DQ415592 DQ415593 DQ415594 DQ415595

B. crucifer

Brazil: Sao Paulo

USNM 303015

DQ415570

DQ415596

B. marinus

Ecuador: Loja NA Vilcabamba

KU 217482

DQ415571

DQ415597

B. schneideri

Paraguay: San Luis de la Sierra

KU 289057

DQ415572

DQ415598

We tested alternative phylogenetic hypotheses to support or refute particular biogeographic scenarios under both parsimony and likelihood conditions. Constraints were designed to 1894

test putative biogeographic barriers inferred from the unconstrained phylogenetic hypotheses. Alternative topologies were constructed using MacClade and implemented as constraints

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo) in paup* and run for 100 random-addition replicates, in the same manner as the unconstrained analyses (see Appendix 1 for constraint topologies). Under parsimony conditions, alternative hypotheses were tested by finding the most parsimonious trees compatible with the imposed constraints and comparing them with the unconstrained MP trees using the Wilcoxon signed-rank tests (Templeton, 1983; Felsenstein, 1985). One- and two-tailed Wilcoxon signed-rank tests were conducted because one-tailed probabilities are close to the exact but are not always conservative, whereas two-tailed test are always conservative (Felsenstein, 1985). We used the S–H (Shimodaira & Hasegawa, 1999) topology test to compare the constrained and unconstrained ML trees, with one-tailed (RELL) distributions and 1000 bootstrap replicates. This test is preferred over the K–H (Kishino & Hasegawa, 1989) test (Goldman et al., 2000). However, it may be too conservative in some cases (Buckley et al., 2001), and borderline results should be used with caution (Townsend et al., 2004). RESULTS Phylogenetic analyses of B. valliceps and B. nebulifer Of the 19 recently collected individuals included from the B. valliceps/nebulifer data set, 14 represented unique haplotypes when sampled across both genes (i.e. some individuals were identical for the 16S sequence, but had different cyt b sequences, and vice versa). When compared with the data of Mulcahy & Mendelson (2000) this resulted in a total of 32 unique haplotypes for the B. valliceps/nebulifer data set. Two individuals from locality Mx8 were identical (Mx8d and h), two individuals (Mx8i and j) were identical to samples Mx1 to Mx2, and four unique haplotypes were recovered from locality Mx10 (Mx10d, f, g were identical; Mx10e was identical to Mx3c from Los Tuxtlas). All phylogenetic analyses were conducted on a data set consisting of the 32 unique haplotypes and three outgroup taxa (B. coccifer, B. mazatlanensis and B. campbelli). Analysis of the 16S data set (562 characters, 73 variable, 26 parsimony-informative) produced 3767 trees (103 steps each); tree not shown. A strict consensus of these trees contained one clade, consisting of all the samples of B. nebulifer, with the remaining samples in an unresolved polytomy (tree not shown). The cyt b data set (421 characters, 105 variable, 57 parsimony-informative) produced six trees (155 steps); tree not shown. A strict consensus of these trees is generally consistent with the overall combined MP tree (see below). The partition-homogeneity tests revealed the two data sets were significantly different (P ¼ 0.04); however, this test may be too conservative at the a ¼ 0.05 level and combining weakly-incongruent data sets generally provides a more accurate estimation (Cunningham, 1997). Therefore, we combined the two generegions into a single analysis. The MP analysis (100 random additions) of 35 operational taxonomic units (OTUs) resulted in two most-parsimonious trees (262 steps, RI ¼ 0.88, CI ¼ 0.79). The only difference

between the two trees was the placement of a sub-clade containing samples of B. nebulifer from Texas. This clade was placed basal to the remaining samples of B. nebulifer in one tree, and nested among them in another. A strict consensus of the two trees was nearly identical to the ML and BI analyses (see below), and supports two monophyletic clades referable to B. valliceps and B. nebulifer. These two clades were supported by bootstrap values of 97 and 100, respectively (Fig. 5). There was relatively little phylogenetic structure among samples of B. nebulifer; however, two sub-clades within B. valliceps show some suggestion of a phylogenetic break that is geographically consistent with the Isthmus of Tehuantepec. Two of the three samples from Los Tuxtlas, Veracruz (north of the isthmus, samples Mx3a and b), were placed in a clade with all other samples found south of the Isthmus of Tehuantepec, while one (Mx3c) grouped with those north of the isthmus (Fig. 5). Modeltest selected TrN + G (gamma shape parameter) model under the hierarchical likelihood ratio test (hLRT) with base frequencies of A ¼ 0.3016; C ¼ 0.2366; G ¼ 0.1789; T ¼ 0.2828, six substitution types with a rate-matrix of A–G ¼ 7.3187; C–T ¼ 11.4953 (all others equal to 1.0), and a gamma shape parameter: G ¼ 0.0713. The ML analysis run with these parameters produced one tree with a score of )ln L ¼ 2723.60976 (Fig. 5). The results of this analysis were similar to those of the MP analysis, with both B. nebulifer and B. valliceps being monophyletic, and with some support for an isthmian break in B. valliceps (again slightly obscured by samples from Los Tuxtlas, Veracruz). In the Bayesian analyses, the three searches each appeared to reach stationarity near 15,000 generations, therefore the first 2000 trees representing conservatively the first 20,000 generations were discarded as the burn-in process. Average tree scores ()ln L ¼ 2826.1891) were taken from each of the 48,000 trees. The 50% majorityrules consensus trees were largely consistent with the MP, and identical to the ML analyses; results for the first BI run are shown in Fig. 5. The outgroup taxon B. coccifer was removed in Fig. 5 to focus on branch-lengths within the ingroup. Uncorrected sequence divergence between B. valliceps and B. nebulifer ranged from 2.8% to 4.0%; similarly, the corrected differences ranged from 2.8% to 4.1% (Table 2). Corrected sequence divergence (HKY85) for the two clades of B. valliceps on either side of the isthmus ranged from 1.0% to 2.0% (Table 2). It is worth noting that the model of evolution used in the ML phylogenetic analyses of Mulcahy & Mendelson (2000) was the Hasegawa–Kishino–Yano (HKY85; Hasegawa et al., 1985). That model was chosen to account for among-site rate variation in the data (Mulcahy & Mendelson, 2000), prior to the wide use of Modeltest. In the present study, we used the program Modeltest and a different model of evolution was selected for the data (TrN + G; see above). The results from the phylogenetic analyses based on the two different models are very similar; however, the corrected sequence divergences between taxa are very different. For example, the corrected sequence divergence average between B. valliceps and B. nebulifer is 3.3% using the HKY85 model (Table 2), whereas the average difference using the TrN + G

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1895

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III

Figure 5 Bayesian analysis phylogram for the Bufo valliceps/nebulifer data set. Topology was nearly identical to the MP and ML analyses. The outgroup B. coccifer was secondarily removed to focus on branch lengths within the ingroup. Parsimony bootstrap values are shown above and Bayesian posterior probabilities (· 100) are shown below significant clades. Bars on right indicate geographic distribution of haplotypes relative to the TMNB and the Isthmus of Tehuantepec. Table 2 Pair-wise percentage sequence divergence within and between Bufo nebulifer and B. valliceps samples (excluding identical haplotypes); averages followed by ranges in parentheses. Comparisons with HKY85-corrected distances (see text) are shown in italics, others are uncorrected. Comparisons within B. valliceps are shown within and between those found north and south of the Isthmus of Tehuantepec (with Mx3a–b in the south) B. valliceps

B. nebulifer B. valliceps North of Isthmus South of Isthmus

B. nebulifer

B. valliceps

North of Isthmus

South of Isthmus

0.4% (0.1–0.8%) 3.3% (2.8–4.1%) – –

3.2% (2.8–4.0%) 1.0% (0.1–2.1%) – –

– – 0.3% (0.1–0.6%) 1.5% (1.0–2.1%)

– – 1.5% (1.0–2.0%) 0.7% (0.1–1.6%)

model is 5.3% (range ¼ 4.1–7.4%), which changes interpretations using the molecular clock. We used HKY85-corrected distances to compare sequence data between lineages in this study (Table 2) to be consistent with our previous study; however, we present the TrN + G model calculations for a comparison. Using the rate of molecular evolution (0.33% per lineage, per million year) from Mulcahy & Mendelson (2000) 1896

and the HKY85-corrected data, we calculate – from Table 2, 3.3% divided by 2 (per lineage), divided by 0.33 (the rate of evolution) – an average c. 5.0 ma (range ¼ 4.2–6.2 ma using min.–max., Table 2) separation of B. valliceps and B. nebulifer and c. 8.0 ma (range ¼ 6.2–11.2 ma) using the TrN + G corrected data. Using the HKY85 calculations, the divergence at the Isthmus of Tehuantepec in B. valliceps is estimated to be

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo) c. 2.3 ma (range ¼ 1.5–3.2 ma), and c. 2.7 ma (range ¼ 1.7– 4.1 ma) using the TrN + G corrected data (1.8%, range ¼ 1.1–2.7% sequence divergence). Phylogenetic analyses of B. marinus Of the 23 individuals from the B. marinus data set, 19 had unique haplotypes when sampled across both genes. Analysis of the 16S data set (548 characters, 50 variable, 30 parsimonyinformative) produced 4 trees (56 steps); tree not shown. Analysis of the cyt b data set (425 characters, 93 variable, 42 parsimony-informative) produced 70 trees (131 steps), the strict consensus of which was generally consistent with the combined MP tree (see below). The partition-homogeneity tests revealed that the two data sets were not significantly different (P ¼ 0.90) and the two were combined into a single analysis. The combined MP analysis (100 random additions) of 22 OTUs resulted in 50 most-parsimonious trees (189 steps, RI ¼ 0.81, CI ¼ 0.80). Bootstrap values for the major clades are shown in Fig. 6. We found phylogenetic signal associated with the break across the TMNB, however the clades were only weakly supported, and showed very little genetic divergence across this break (see below). One clade, containing the samples north of the TMNB (Mx8) had low bootstrap support (63), and was sister to the clade of samples from the middle

(Mx9; Fig. 6). The remaining samples along the transect (Mx10) were placed in a polytomy with those from Guerrero (Mx11) and Sinaloa (Mx12), Mexico. In the MP analysis, we found support for a break associated with the Isthmus of Tehuantepec (Fig. 6), with the exception of one sample (G1) taken from the head of the Grijalva Valley, in Guatemala. This sample was placed among the westerly samples from Mexico, and the support for this clade was weak (< 50%; Fig. 6); however, the clade found southeast of the Isthmus was wellsupported. For the two genes combined, Modeltest selected TrN + G (gamma shape parameter) as the best-fit model under the hierarchical likelihood ratio test (hLRT) criteria, with base frequencies of A ¼ 0.2855; C ¼ 0.2417; G ¼ 0.1824; T ¼ 0.2904, six substitution types with a rate-matrix of A–G ¼ 7.8102; C–T ¼ 14.1003 (all others equal to 1.0), and a gamma shape parameter: G ¼ 0.1042. The ML analysis run with these parameters produced a tree a score of )ln L ¼ 2355.4517. Our ingroup samples all had relatively short branch lengths, and the results of this analysis were generally similar to that of the MP analysis. We found similar support for a phylogenetic break across the TMNB and the Isthmus of Tehuantepec, with sample G1 placed in the clade of samples from Mexico (west of the isthmus). In this analysis, B. schneideri rendered B. marinus paraphyletic with respect to

Figure 6 Bayesian analysis phylogram for Bufo marinus. Topology was nearly identical to the MP and ML analyses. Parsimony bootstrap values are shown above and Bayesian posterior probabilities (· 100) are shown below for significant clades. Bars on right indicate geographic distribution of haplotypes relative to the TMNB and the Isthmus of Tehuantepec. Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1897

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III our sample from South America and the remaining Mesoamerican samples (not shown). For each gene-region tested separately, Modeltest selected the HKY + G (gamma shape parameter) model under the hierarchical likelihood ratio test (hLRT) criteria. Therefore, the three Bayesian analyses were run, each with a 2-substitutiontype (transition/transversion) plus gamma model, and searches appeared to reach stationarity near 10,000–12,000 generations. The first 2000 trees representing conservatively the first 20,000 generations were discarded as the burn-in process. Average tree scores ()ln L ¼ 2400.1138) were taken from the 48,000 trees of each run. The 50% majority-rules consensus tree and posterior-probabilities for the first run are shown in Fig. 6. The results were largely consistent with those of the MP and ML analyses, with strong support for a monophyletic Mesoamerican clade of B. marinus low support for a break across the TMNB, and strong support for an isthmian break, again with the exception of the position of G1. Sequence divergence for B. marinus across the TMNB ranged from 0.4% to 0.7% (Table 3). Using the difference between north and south of the TMNB, we estimated a date of divergence of c. 0.9 ma (range ¼ 0.6–1.1). Sequence divergence between samples of B. marinus on either side of the Isthmus of Tehuantepec were slightly larger, ranging from 0.8% to 2.2% (Table 3), estimated as c. 2.7 ma (range ¼ 1.2–3.3). Sequence variation was much greater among samples southeast of the Isthmus of Tehuantepec (0.1–2.2%) as compared to variation found north-west of the isthmus (0.1–1.3%), where the greatest variation was found between sample Mx12 and Mx8–10 (0.8– 1.3%); Mx12 differed from Mx11 by 0.8%. The sample from South America was very different from the Mesoamerican samples (6.3–7.4%; Table 3), with an estimated date of divergence ranging from c. 10.9 ma (range ¼ 9.5–11.2). Topology tests for both species groups Although clades associated with two major geographic features (putative barriers) were recovered in both species groups, the overall structure among those clades was strikingly different. The timing of divergence associated with each of the respective geographic barriers was different for each species group. The

North of TMNB Middle South of TMNB

North of Isthmus South of Isthmus

Central America South America

1898

B. valliceps/nebulifer group showed divergence across the TMNB prior to separation across the Isthmus of Tehuantepec, while the B. marinus group showed the opposite (Fig. 7). However, in both groups, clades associated with either side of the isthmus were not reciprocally monophyletic. We tested alternative phylogenies in order to validate the strength of the effect of these barriers on our data. First, in order to test the effect of the TMNB on the B. valliceps/ nebulifer group, we constrained two B. valliceps haplotypes (Mx9a–b) to be in the B. nebulifer clade and conducted the same analyses that were performed on the unconstrained data. This resulted in six equally-parsimonious trees, each of which were significantly different from either of the unconstrained trees and the ML constrained topology was also significantly different (‘TMNB’, Table 4). Therefore we could reject paraphyly of B. valliceps with respect to the samples at locality Mx9. Secondly, to test the affect of an isthmian break in B. valliceps, we constrained haplotypes Mx3a–c to be in the clade north-west of the isthmus, but south of the TMNB (with Mx9–10). This resulted in 12 trees that were each significantly different from the unconstrained trees. However, the ML constrained tree had an )ln L score difference of 21.740, which was not significantly different (‘Isthmus’, Table 4). Here, we reject reciprocally monophyly of B. valliceps haplotypes on either side of the Isthmus of Tehuantepec, but with caution. Thirdly, to test for the alternative scenario recovered in the B. marinus group, we constrained the B. valliceps/ nebulifer group to have the divergence associated with the isthmus to precede that of the TMNB; i.e., enforcing the topology of Fig. 7b on the B. valliceps/nebulifer group, consequently rendering B. valliceps paraphyletic: (B. valliceps S. of Isthmus (B. valliceps N of Isthmus + B. nebulifer)). This analysis resulted in four MP trees that were significantly different from either of the unconstrained trees based on the parsimony analysis (‘MNB vs. Isthmus’, Table 4). However, the constrained ML topology had a likelihood score that differed by 27.37, which was not significantly different based on the S–H topology test (Table 4). Based on these results, we tentatively reject the hypothesis that the B. valliceps/nebulifer group shows structure similar to that of B. marinus associated with the TMNB and the Isthmus of Tehuantepec, because the

North of TMNB

Middle

South of TMNB

0.1% (0.1–0.2%) 0.4% (0.3–0.5%) 0.6% (0.4–0.7%)

0.4% (0.3–0.5%) 0.2% (0.1–0.3%) 0.7% (0.5–0.8%)

0.6% (0.4–0.7%) 0.7% (0.5–0.8%) 0.62%

North of Isthmus

South of Isthmus

0.6% (0.1–1.3%) 1.8% (0.8–2.2%)

1.7% (0.8–2.2%) 1.1% (0.1–2.2%)

Central America

South America

1.2% (0.1–2.2%) 7.2% (6.7–7.9%)

6.8% (6.3–7.4%) –

Table 3 Pair-wise percentage sequence divergence within and between Bufo marinus samples (excluding identical haplotypes); averages followed by ranges in parentheses. Comparisons with HKY85-corrected distances are shown in italics, others are uncorrected. Samples from north and south of the TMNB represent localities Mx8 and Mx10, respectively, whereas those in the middle are from locality Mx9. There is only one withincomparison of those south of the TMNB. Comparisons made with those north and south of the Isthmus of Tehuantepec include Mx11–12 as north of the Isthmus

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo)

Figure 7 Area cladograms recovered for the two species groups in this study. (a) Area-cladogram for the Bufo valliceps/nebulifer species group, showing that the TMNB divergence occurred prior to the separation of B. valliceps haplotypes at the Isthmus of Tehuantepec. (b) Area-cladogram for B. marinus, showing the opposite pattern of the B. valliceps/nebulifer species group.

ML test was marginally non-significant (Buckley et al., 2001; Townsend et al., 2004). Likewise, we conducted similar topology tests for the B. marinus group. Here, we constrained samples north of the TMNB (Mx8–9) to be monophyletic, and those south of the TMNB and north of the isthmus (Mx10–12, and G1) to be

monophyletic and obtained 25 trees, none of which were significantly different from the 50 unconstrained MP trees. The ML constrained topology was also not significantly different (‘TMNB’, Table 4). Therefore we could not reject a separation of haplotypes across the TMNB; however, we note that the level of divergence across the TMNB in B. marinus haplotypes is much less than that observed in the B. valliceps/nebulifer group (0.6% vs. 3.3%, respectively; Tables 2 & 3). Secondly, we constrained clades on either side of the Isthmus of Tehuantepec to be monophyletic by forcing G1 with the samples originally recovered in the clade south of the isthmus (G3, E1– 2, H3–4, C1). This resulted in 35 MP trees and an ML tree with a likelihood score that had a difference of 7.08; none of these trees was significantly different from the unconstrained trees (‘Isthmus’, Table 4). Therefore, we could not reject an effect of the Isthmus of Tehuantepec on B. marinus haplotypes. Lastly, we constrained the B. marinus group to show the same overall topological structure as the B. valliceps/nebulifer group (i.e. enforcing the topology of Fig. 7a on B. marinus: (Mx8 (Mx9– 12, G1) + (G3, E1–2, H3–4, C1)). The results of both the MP and ML tests were significantly different from the unconstrained trees (‘TMNB vs. Isthmus’, Table 4). Therefore we could reject the hypothesis that B. marinus shows an overall pattern similar to that of the B. valliceps/nebulifer group.

Table 4 Results of the topology tests. Tree scores and number of equally parsimonious trees recovered in the unconstrained phylogenies are shown for both groups in the left column. Columns to the right show results for the constrained phylogenies, with the range of summary statistics (see text for explanation of constraint topologies and Appendix 1 for actual constraints enforced). Unconstrained phylogenies

Constrained phylogenies

Bufo valliceps/nebulifer

TMNB

Isthmus

TMNB vs. Isthmus

(6 trees, 274 steps)

(12 trees, 271 steps)

(4 trees, 280 steps)

Parsimony

n 

Zà

P§

n

Z

P

n

Z

P

(2 trees, 262 steps)

11 12

3.21 3.46

0.0007** 0.0003**

11

2.18

0.0147**

17 19

4.02 3.38

0.0001** 0.0007**

Likelihood

Constrain

Diff.

P

Constrain

Diff.

P

Constrain

Diff.

P

)ln L ¼ 2723.6098

2744.3887

20.78

0.035*

2745.3506

21.740

0.085

2750.9784

27.37

0.057

Bufo marinus

(25 trees, 190 steps)

(35 trees, 192 steps)

(20 trees, 194 steps)

Parsimony

n

Z

P

n

Z

P

n

Z

P

(50 trees, 189 steps)

1 3 5

1.00 0.58 0.45

0.1587 0.2819 0.3274

5 7

1.34 1.13

0.0899 0.1284

5 7

2.24 1.89

0.0127** 0.0294*

Likelihood

Constrain

Diff.

P

Constrain

Diff.

P

Constrain

Diff.

P

)ln L ¼ 2355.4517

2359.5382

4.08

0.303

2362.5321

7.08

0.143

2374.5949

19.14

0.045*

 Number of characters differing in minimum number of changes on paired topologies. àNormal approximation for Wilcoxon signed-ranks tests. §Asterisks indicate a significant difference between the overall shortest tree and the constrained topologies. One asterisk denotes significance using the one-tailed probability and two asterisks denote significance using the two-tailed probability. One-tailed probabilities are shown, two-tailed are twice the values. Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1899

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III DISCUSSION Our data support a role for two major geographic features in shaping the historical biogeography in two independent lineages of lowland toads. However, the sequence of events appears different for each group (Fig. 7). The TMNB has been shown to represent a barrier to gene flow in other lowland taxa (Hulsey et al., 2004; Zaldı´var-Rivero´n et al., 2004), while the Isthmus of Tehuantepec has been more typically invoked as a barrier for upland taxa (Marshall & Liebherr, 2000; Sullivan et al., 2000). In our study, we confirmed that the B. valliceps/ nebulifer species boundary lies along the eastern terminus of the TMNB, consistent with the Pliocene-vicariance hypothesis of Mulcahy & Mendelson (2000) and identified subsequent phylogeographic structure in B. valliceps associated with the Isthmus of Tehuantepec. In contrast, B. marinus showed phylogeographic structure along either side of the Isthmus of Tehuantepec, prior to that of the TMNB. At the TMNB, the signal in the B. marinus data is more consistent with the ‘Pleistocene-dispersal, followed by vicariance hypothesis’ of Mulcahy & Mendelson (2000). We evaluated the veracity of these contrasting patterns between the B. valliceps/nebulifer and B. marinus groups by using alternative topology tests. We imposed the B. marinus phylogenetic structure (Fig. 7b) on the B. valliceps/nebulifer group, and vice versa, and the resulting topologies were significantly different, based on the MP analyses, from the shortest trees recovered in each group (Table 4). The ML topology test recovered marginally nonsignificant results for the B. valliceps/nebulifer group with respect to the isthmus (Table 4). However, this result is in conflict with a well-supported node based on the MP and Bayesian analyses (97 and 96 respectively; Fig. 5), which suggests that the ML-based (S–H) topology test may be more conservative than the MP topology test (Wilcoxon signedrank). Therefore, we reject the hypothesis that these sympatric species groups, with similar natural histories, show the same phylogeographic structure. Previous molecular studies invoking vicariance for related taxa distributed on either side of the TMNB consisted of data sampling far to the north and south, and inferred clade boundaries to be at the eastern terminus of the TMNB (Mulcahy & Mendelson, 2000; Hulsey et al., 2004; Zaldı´varRivero´n et al., 2004). Our additional samples of B. valliceps and B. nebulifer across the TMNB (Figs 2 & 3) demonstrate the strict parapatry on either side of this now seemingly minor geographic barrier. Our sample size in the middle of the TMNB is relatively small (n ¼ 2), such that future work may reveal these two species to be sympatric in the small area of suitable habitat near the town of El Viejo´n (Mx9, Fig. 2). Based on the fact that species of Bufo are frequently interfertile (Blair, 1963; Masta et al., 2002), it would not be surprising to us if these species were found to hybridize in this area. If this were the case, we may not be able to detect the presence of both species with mtDNA as our only genetic marker. Nonetheless, our data strongly indicate that these northern and southern clades represent distinct evolutionary species, showing a 1900

precisely parapatric distribution across this ancient geographic barrier. Our increased sampling south of the TMNB greatly improved the phylogenetic support for the monophyly of B. valliceps, which was previously lacking (Mulcahy & Mendelson, 2000). We attribute this difference to the increased sampling between the TMNB and the Isthmus of Tehuantepec, and the phylogenetic signal associated with the isthmus. The sampling of Mulcahy & Mendelson (2000) contained only three samples of B. valliceps north of the isthmus (Mx3a–c), all from Los Tuxtlas. These three haplotypes were scattered across the distribution of other haplotypes of B. valliceps, with Mx3c recovered in the most basal position (Mulcahy & Mendelson, 2000). With increased sampling (this study), Mx3c was placed in a clade of haplotypes found north of the Isthmus of Tehuantepec (Fig. 5), a clade estimated to be of late Pliocene–Pleistocene origin. However, samples Mx3a–b remain nested among haplotypes in the clade south of the isthmus. We were able to reject the alternative hypothesis that samples north of the isthmus formed a monophyletic clade (Table 4). At this time, we can not discriminate between incomplete lineage sorting vs. a recent dispersal for the occurrence of haplotypes Mxa–b at Los Tuxtlas. However, because the Los Tuxtlas landscape is of late Pliocene volcanic origin (Ferrusquia-Villafranca, 1993) – the approximate divergence time for the B. valliceps isthmian clades – we favour a more recent dispersal explanation for the occurrence of haplotypes from both clades at Los Tuxtlas. Our analyses of the B. marinus data do not support a hypothesis that the Pliocene uprising of the TMNB caused a major vicariant event among populations of this species, as it did in the B. valliceps/nebulifer group. The data from B. marinus do contain phylogeographic structure associated with the TMNB; however, the branch lengths are very shallow (Fig. 6) and the sequence divergences are low (Table 3). Based on these results, we suggest that B. marinus dispersed across this barrier during lower sea levels associated with glacial maxima during the Pleistocene. Our rate of evolution estimated from the B. valliceps/nebulifer group corroborates this hypothesis. If we calibrate a rate for B. marinus based on a Pliocene (5.4 ma) divergence across the TMNB, consistent with that of the B. valliceps/nebulifer group (Mulcahy & Mendelson, 2000; this study), we obtain a rate of c. 1.6%, per lineage, per million years, which is twice the expected rate for ectotherms (Tan & Wake, 1995), and far beyond that estimated for other bufonid mtDNA (Graybeal, 1993; Macey et al., 1998). Bufo marinus is known to be of South American origin (i.e. all of its closest relatives are included in clades restricted to South America; Slade & Moritz, 1998; Pauly et al., 2004; Pramuk, 2006). Our results are consistent with a model in which B. marinus dispersed into Central America during the early Pliocene, and later into southern Mexico, subsequent to the uplift of the TMNB. A fossil record referred to B. marinus from the Miocene of Kansas, USA, (Wilson, 1968) refutes this scenario. However, recent concern has questioned the reliability of precise taxonomic assignment based on fragmentary fossil data in anurans (Bever, 2005; see also Pauly et al., 2004), and the

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo) specific designation of this record remains questionable (Holman, 2003). Data from both B. valliceps and B. marinus show phylogenetic structure associated with the Isthmus of Tehuantepec. While our sampling is not fully adequate to test for an effect of a Pliocene seaway across the Isthmus of Tehuantepec, our data do provide two independent lines of evidence to suggest that it may have been present. Our estimated dates of divergence from clades on either side of the isthmus were approximately 2–3 ma, in both species groups. This suggests a common barrier of some type existed for lowland taxa at the Isthmus of Tehuantepec during the mid-Pliocene. Beard et al. (1982) suggest sea-level fluctuations associated with continental glaciation began approximately 3 ma, indicating high sea levels prior to the first glacial cycle beginning c. 2.8 ma. This corresponds to our timing of divergence in both B. valliceps and B. marinus at the isthmus. Sullivan et al. (2000) also found sequence divergence, in harvest mice (Reithrodontomys), consistent with a late Pliocene trans-isthmus marine barrier. We were able to reject monophyletic lineages on either side of the isthmus for B. valliceps; however, we were not able to reject this alternative in B. marinus. Bufo marinus is abundant in the lowlands of the isthmus and indeed in all adjacent lowland areas (Duellman, 1960; J.R. Mendelson, pers. obs.). Bufo valliceps however, is common along the Atlantic lowlands northeast of the isthmus, absent along the Pacific Coast of Oaxaca and Guerrero to the southwest of the isthmus, and is inexplicably uncommon along the Pacific Coastal Plain of Chiapas and Guatemala to the southeast of the isthmus (Mendelson, 2001). This biogeographical contrast between the two species suggests that B. marinus may disperse more readily into new regions than does B. valliceps. The exception to the isthmian clades in our B. marinus data (sample G1) may be an artifact of poor sampling in this area, or may represent a recent dispersal from the northwest into the Grijalva Valley, similar to that observed in B. valliceps (cf. Fig. 3). Further sampling of multiple taxa, including B. valliceps and B. marinus should be conducted across this region, in order specifically to test hypotheses associated with a Pliocene seaway across the Isthmus of Tehuantepec, and subsequent dispersal routes across this barrier. The three species of toads in our study are similar in having widespread, lowland distributions, and by having great tendencies to disperse into a variety of habitats. Despite these evident similarities in their natural histories, our data indicate that they have experienced quite different histories with respect to the formation of the TMNB and the possible inundation of the Isthmus of Tehuantepec. These differences are most likely attributable to matters of timing – the ancestor of B. nebulifer + B. valliceps predates the major orogeny of the TMNB in this area, while B. marinus is likely to have arrived long after it represented a significant barrier to dispersal. Most taxa examined to date show species-level divergence across the TMNB (Hulsey et al., 2004; Zaldı´var-Rivero´n et al., 2004). Here, we have shown low levels of sequence divergence across this barrier within a single species–B. marinus. Pe´rez-Higareda

& Navarro (1980) point out some 16 species of reptiles that straddle this barrier and many others that show subspecific delineations on either side, yet these remain to be tested with molecular data. Both B. valliceps and B. marinus show low levels of intraspecific variation, possibly associated with a midPliocene barrier across the Isthmus of Tehuantepec. The question remains as to why neither B. nebulifer nor B. valliceps have secondarily dispersed over the TMNB, concurrent with B. marinus. Perhaps our limited sampling of the area in the middle (or our choice of molecular marker) was not sufficient to detect sympatry in these species. Nevertheless, we were able to pin-point a putative area of contact between these two species for future investigation. In a logical sequence then, our results demonstrate an ancient vicariant event, a more recent dispersal event, and the area of possible ongoing biological influences maintaining species boundaries (viz., between B. nebulifer and B. valliceps). ACKNOWLEDGEMENTS We would like to thank the following people, institutions and curators for providing tissue samples and locality data for this study: Jonathan A. Campbell, Kevin de Queiroz, Oscar FloresVillela, Carl J. Franklin, Eli B. Greenbaum, John H. Malone, Jenny B. Pramuk, Mahmood Sasa M., John E. Simmons, Eric N. Smith, Larry D. Wilson, UTACV, KUNHM, USNM and UNAM. We thank Paul Wolf, Mike Pfrender, Eric O’Neill, and Chris Feldman for assistance in the lab, and Chris Feldman and Bob Macey for assistance in data analysis, Karen R. Lips and Ron M. Bonett for field assistance, and Ted Papenfuss for use of his Mac (G5) for phylogenetic analyses. USU Department of Biology and an Undergraduate Research and Creative Opportunities Grant (URCO), awarded to B. H. Morrill, provided funding the data collection. Fieldwork along the transect of the TMNB was funded by a Theodore Roosevelt Memorial Grant from the American Museum of Natural History awarded to D. G. Mulcahy and a National Geographic grant to J. R. Mendelson and K. R. Lips. REFERENCES Beard, J.H., Sangree, J.B. & Smith, L.A. (1982) Quaternary chronology, paleoclimate, depositional sequences, and eustatic cycles. American Association of Petroleum Geologists Bulletin, 66, 158–169. Bever, G.S. (2005) Variation in the ilium of North American Bufo (Lissamphibia; Anura) and its implications for specieslevel identification of fragmentary anuran fossils. Journal of Vertebrate Paleontology, 25, 548–560. Blair, W.F. (1963) Evolutionary relationships of North American toads of the genus Bufo: a progress report. Evolution, 17, 1–16. Bryant, W.R. & Bryant, J.R. (1991) Bathymetric Chart: Gulf of Mexico Region. The Gulf of Mexico Basin: the geology of North America (ed. by A. Salvador), Vol. J, Plate 1. Geological Society of America, Boulder, CO.

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1901

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III Bryant, W.R., Lugo, J., Cordova, C. & Salvador, A. (1991) Physiography and bathymetry. The Gulf of Mexico Basin: the geology of North America (ed. by A. Salvador), Vol. J, pp. 13– 30. Geological Society of America, Boulder, CO. Buckley, T.R., Simon, C., Shimodaira, H. & Chambers, G.K. (2001) Evaluating hypotheses on the origin and evolution of the New Zealand alpine cicadas (Maoricicada) using multiple-comparison tests of tree topology. Molecular Biology and Evolution, 18, 223–234. Campbell, J.A. (1998) Amphibians and reptiles of northern Guatemala, the Yucatan and Belize. University of Oklahoma Press, Norman. Campbell, J.A. (1999) Distribution patterns of amphibians in Middle America. Distribution patterns of amphibians: a global perspective (ed. by W.W. Duellman), pp. 111–209. John Hopkins University Press, Baltimore. Campbell, J.A. & Frost, D.R. (1993) Anguid lizards of the genus Abronia: revisionary notes on the species of Nuclear Central America and adjacent Mexico, descriptions of four additional species, with a phylogenetic hypothesis for the genus and an identification key. Bulletin of the American Museum of Natural History, 216, 1–121. Castoe, T.A., Chippindale, P., Campbell, J.A., Ammerman, L.K. & Parkinson, C.L. (2003) Evolutionary relationships of the Middle American Jumping Pitvipers (Genus Atropoides) and phylogeography of the Atropoides nummifer complex. Herpetologica, 49, 421–432. de Cserna, Z. (1989) An outline of the geology of Mexico. The geology of North America (ed. by A. Salvador), Vol. A, pp. 223–264. Geological Society of America, Boulder, CO. Cunningham, C.W. (1997) Can three incongruence tests predict when data should be combined? Molecular Biology and Evolution, 14, 733–740. Darda, D.M. (1994) Allozyme variation and morphological evolution among Mexican salamanders of the genus Chiropterotriton (Caudata: Plethodontidae). Herpetologica, 50, 164–187. Duellman, W.E. (1960) A distributional study of amphibians of the Isthmus of Tehuantepec, Mexico. University of Kansas Publication, Museum of Natural History, 13, 19–72. Ewing, T.E. & Lopez, R.F. (1991) Principle structural features. The Gulf of Mexico Basin: the geology of North America (ed. by A. Salvador), Vol. J, Plate 2. Geological Society of America, Boulder, CO. Felsenstein, J. (1985) Confidence limits on phylogenies with a molecular clock. Systematic Zoology, 34, 152–161. Ferrusquia-Villafranca, I. (1993) Geology of Mexico: a synopsis. Biological diversity of Mexico: origins and distribution (ed. by T.P. Ramamoorthy, R. Bye and A. Lot), pp. 3–107. Oxford Press, New York. Goldman, N., Anderson, J.P. & Rodrigo, A.G. (2000) Likelihood-based tests of topologies in phylogenetics. Systematic Biology, 49, 652–670. Graybeal, A. (1993) The phylogenetic utility of cytochrome b: lessons from bufonid frogs. Molecular Phylogenetics and Evolution, 2, 256–269. 1902

Graybeal, A. (1997) Phylogenetic relationships of bufonid frogs and tests of alternate macroevolutionary hypotheses characterizing their radiation. Zoological Journal of the Linnaean Society, 119, 297–338. Hasegawa, M,, Kishino, H. & Yano, T. (1985) Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution, 21, 160–174. Holman, J.A. (2003) Fossil frogs and toads of North America. Indiana University Press, Bloomington, IN. Huelsenbeck, J.P. & Ronquist, F. (2001) MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics, 17, 754–755. Hulsey, C.D., Garcia de Leon, F.J., Johnson, Y.S., Hendrickson, D.A. & Near, T.J. (2004) Temporal diversification of Mesoamerican cichlid fishes across a major biogeographic boundary. Molecular Phylogenetics and Evolution, 31, 754–764. Kishino, H. & Hasegawa, M. (1989) Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea. Journal of Molecular Evolution, 29, 170–179. Lee, J.C. (1996) The amphibians and reptiles of the Yucata´n Peninsula. Comstock Publishing Co., Ithaca, New York. Macey, R.J., Schulte, J.A., II, Larson, A., Fang, Z., Wang, Y., Tuniyev, B.S. & Papenfuss, T.J. (1998) Phylogenetic relationships of toads in the Bufo bufo species group from the eastern escarpment of the Tibetan Plateau: A case of vicariance and dispersal. Molecular Phylogenetics and Evolution, 9, 80–87. Maddison, W.P. & Maddison, D.R. (2000) MacClade, analysis of phylogeny and character evolution, Version 4.0. Sinauer, Sunderland, MA. Maniatis, T., Fristch, E.F. & Sambrook, J. (1982) Molecular cloning: a laboratory manual. Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY. Marshall, C.J. & Liebherr, J.K. (2000) Cladistic biogeography of the Mexican transition zone. Journal of Biogeography, 27, 203–216. Masta, S.E., Sullivan, B.K., Lamb, T. & Routman, E.J. (2002) Molecular systematics, hybridization, and phylogeography of the Bufo americanus complex in Eastern North America. Molecular Phylogenetics and Evolution, 24, 302–314. Mateos, M. (2005) Comparative phylogeography of livebearing fishes in the genera Poeciliopsis and Poecilia (Poeciliidae: Cyprinodontiformes) in central Mexico. Journal of Biogeography, 32, 775–780. McCranie, J.R. & Wilson, L.D. (2002) The amphibians of Honduras. Society for the Study of Amphibians and Reptiles, Ithaca, New York. Mendelson, J.R., III (1994) A new species of toad (Anura: Bufonidae) from the lowlands of eastern Guatemala. Occasional Papers of the Museum of Natural History, University of Kansas, 166, 1–21. Mendelson, J.R., III (2001) A review of the Guatemalan toad Bufo ibarrai (Anura: Bufonidae), with distributional and taxonomic comments of Bufo valliceps and Bufo coccifer. Mesoamerican herpetology: systematics, natural history, and conservation (ed. by J.D. Johnson, R.G. Webb and O. Flores-

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

Biogeography of lowland toads (Bufo) Villela), pp. 10–19. The University of Texas at El Paso, El Paso, TX. Mendelson, J.R., III, Williams, B.L., Sheil, C.A. & Mulcahy, D.G. (2005) Systematics of the Bufo coccifer complex (Anura: Bufonidae) of Mesoamerica. Scientific Papers, University of Kansas Natural History Museum, 38, 1–27. Mulcahy, D.G. & Mendelson, J.R., III (2000) Phylogeography and speciation of the morphologically variable, widespread species Bufo valliceps, based on molecular evidence from mtDNA. Molecular Phylogenetics and Evolution, 17, 173–189. Pauly, G.B., Hillis, D.M. & Cannatella, D.C. (2004) The history of a nearctic colonization: molecular phylogenetics and biogeography of the nearctic toads (Bufo). Evolution, 58, 2517–2535. Pe´rez-Higareda, G. & Navarro, L.D. (1980) The faunistic districts of the low plains of Veracruz, Mexico, based on reptilian and mammalian data. Bulletin of Maryland Herpetological Society, 16, 54–69. Posada, D. & Crandall, K.A. (1998) Modeltest: testing the model of DNA substitution. Bioinformatics, 14, 817–818. Pramuk, J.B. (2006) Phylogeny of South American Bufo (Anura: Bufonidae) inferred from combined evidence. Zoological Journal of the Linnaean Society, 146, 407–452. Rosen, D.E. (1978) Vicariant patterns and historical explanation in biogeography. Systematic Zoology, 27, 159–188. Savage, J.M. (2002) The amphibians and reptiles of Costa Rica. The University Chicago Press, Chicago, IL. Shimodaira, H. & Hasegawa, M. (1999) Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Molecular Biology and Evolution, 16, 1114–1116. Slade, R.W. & 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. Sullivan, J., Arellano, E. & Rogers, D.S. (2000) Comparative phylogeography of Mesoamerican Highland Rodents: concerted versus independent response to past climatic fluctuations. The American Naturalist, 155, 755–768. Swofford, D.L. (1999) PAUP*: Phylogenetic Analysis Using Parsimony (*and other methods), Beta Version 4.0b10. Sinauer, Sunderland, MA. Tan, A.-M. & Wake, D.B. (1995) MtDNA phylogeography of the California newt Taricha torosa (Caudata, Salamandridae). Molecular Phylogenetics and Evolution, 4, 383– 394. Templeton, A.R. (1983) Phylogenetic inference from restriction endonuclease cleavage site maps with particular reference to the evolution of humans and the apes. Evolution, 37, 221–244. Townsend, T.M., Larson, A., Louis, E. & Macey, R. (2004) Molecular phylogenetics of Squamata: the position of snakes, amphisbaenians, and dibamids, and the root of the squamate tree. Systematic Biology, 53, 735–757. Wilson, R.L. (1968) Systematic and faunal analysis of a Lower Pliocene vertebrate assemblage from Trego County, Kansas.

Contributions from the Museum of Paleontology, University of Michigan, 2, 75–126. Zaldı´var-Rivero´n, A., Leon-Regagnon, V. & Nieto-Montes de Oca, A. (2004) Phylogeny of the Mexican coastal leopard frogs of the Rana berlandieri group based on mtDNA sequences. Molecular Phylogenetics and Evolution, 30, 38–49. Zug, G.R. & Zug, P.B. (1979) The marine toad, Bufo marinus: a natural history resume of native populations. Smithsonian Contributions to Zoology, 284, 1–58.

BIOSKETCHES Daniel G. Mulcahy is a PhD candidate in the Department of Biology, at Utah State University (USU). His graduate research involves biogeography and systematics of Neotropical and western North American amphibians and reptiles. Joseph R. Mendelson III received his PhD in 1997 from The University of Kansas, and has been conducting research and conservation studies on Neotropical amphibians and reptiles for 16 years. He is now Curator of Herpetology at Zoo Atlanta, and Adjunct Associate Professor of Biology at Utah State University. Benson H. Morrill is a PhD student in L. Rickords’ molecular biology laboratory at Utah State University. He was an undergraduate at USU and completed his honours thesis as part of this study. His academic interests are herpetology and genetics.

Editor: Brett Riddle

APPENDIX 1 TOPOLOGIES USED IN CONSTRAINT TEST. Constrained topologies are labelled following Table 4 for the three tests for each group: TMNB, Isthmus, and TMNB vs. Isthmus. For the Bufo valliceps/nebulifer data, the following translation is used: 1 B. coccifer, 2 B. mazatlanensis, 3 B. campbelli, 4 Bneb Lou1, 5 Bneb Tx2, 6 Bneb Tx4, 7 Bneb Mx1, 8 Bneb Mx8b, 9 Bneb Mx8c, 10 Bneb Mx8d, 11 Bneb Mx8e, 12 Bneb Mx8a, 13 Bneb Mx8f, 14 Bneb Mx8g, 15 Bval Mx3a, 16 Bval Mx3b, 17 Bval Mx3c, 18 Bval Mx5, 19 Bval Mx6, 20 Bval Mx7, 21 Bval B1, 22 Bval B2, 23 Bval H1, 24 Bval H2, 25 Bval G1, 26 Bval G2, 27 Bval G3, 28 Bval G4a, 29 Bval G4b, 30 Bval Mx10a, 31 Bval Mx10b, 32 Bval Mx10c, 33 Bval Mx10d, 34 Bval Mx9a, 35 Bval Mx9b, with the following topologies were used for the constrained phylogenetic analyses: TMNB ¼ (1, ((2, ((4–14, 34–35), (15– 33))), 3)); Isthmus ¼ (1, ((2, ((4– 14), ((18–29,), (15–17, 30–35)))), 3)); TMNB vs. Isthmus ¼ (1, ((2, ((18–29), ((15–35), (4–14)))), 3)); For the B. marinus data, the following translation is used: 1 B. crucifer, 2 B. schneideri, 3 B. marinus S. A., 4 Bmar Mx11b, 5 Bmar

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd

1903

D. G. Mulcahy, B. H. Morrill and J. R. Mendelson III Mx9f, 6 Bmar Mx9a, 7 Bmar Mx9d, 8 Bmar G3, 9 Bmar Mx9b, 10 Bmar Mx11a, 11 Bmar C1, 12 Bmar Mx12, 13 Bmar Mx8a, 14 Bmar H3, 15 Bmar E1, 16 Bmar Mx9c, 17 Bmar H4, 18 Bmar G1, 19 Bmar Mx8b, 20 Bmar Mx8c, 21 Bmar Mx10b, 22 Bmar Mx10a, with the following topologies were used for the

1904

constrained phylogenetic analyses: TMNB ¼ (1, (2, (3, (((5–7, 9, 13, 16, 19, 20), (4, 10, 12, 18, 21, 22)), (8, 11, 14, 15, 17))))); Isthmus ¼ (1, (2, (3, ((4–7, 9, 10, 13, 16, 19, 20, 22, 12, 21), (8, 11, 14, 15, 17, 18))))); TMNB vs. Isthmus ¼ (1, (2, (3, ((13, 19, 20), ((4–7, 9, 10, 12, 16, 22, 18, 21), (8, 11, 14, 15, 17)))))).

Journal of Biogeography 33, 1889–1904 ª 2006 The Authors. Journal compilation ª 2006 Blackwell Publishing Ltd