The role of mitochondrial introgression in illuminating the evolutionary history of Nearctic treefrogs

bs_bs_banner

Zoological Journal of the Linnean Society, 2014, 172, 103–116. With 4 figures

The role of mitochondrial introgression in illuminating the evolutionary history of Nearctic treefrogs ROBERT W. BRYSON JR1,2*, BRIAN TILSTON SMITH3, ADRIAN NIETO-MONTES DE OCA4, URI OMAR GARCÍA-VÁZQUEZ4 and BRETT R. RIDDLE1 1

School of Life Sciences, University of Nevada Las Vegas, Las Vegas, NV 89154-4004, USA Department of Biology and Burke Museum of Natural History and Culture, University of Washington, Seattle, WA 98195-1800, USA 3 Museum of Natural Science, Louisiana State University, Baton Rouge, LA 70803, USA 4 Museo de Zoología Alfonso L. Herrera, Facultad de Ciencias, Universidad Nacional Autónoma de México, Circuito exterior s/n, Cd. Universitaria, México 04510, Distrito Federal, Mexico 2

Received 17 January 2014; revised 30 April 2014; accepted for publication 1 May 2014

Inferring the evolutionary history of lineages often becomes difficult when gene histories are in conflict with each other. Introgression, for example, can cause DNA sequences from one species to be more similar to sequences of a different species and lead to incongruence amongst gene trees. However, incorporating congruent and incongruent locus-specific phylogenetic estimates with the geographical distribution of lineages may provide valuable insight into evolutionary processes important to speciation. In this study, we investigated mitochondrial introgression within the Hyla eximia group to better understand its role in illuminating the evolutionary history and phylogeography of these treefrogs. We reconstructed and compared the matrilineal history of the Hyla eximia group with estimates of evolutionary history inferred from nuclear genes. We tested for introgression within the mitochondrial and nuclear genes using a posterior predictive checking approach. Reconstructions of the species tree based on the mitochondrial DNA (mtDNA) and nuclear DNA data were strongly discordant. Introgression between lineages was widespread in the mtDNA data set (145 occurrences amongst 11 of the 16 lineages), but uncommon in the nuclear genes (12 occurrences amongst four of the 16 lineages). Nonetheless, the geographical structuring of mtDNA within species provides valuable information on biogeographical areas, ancient areas of hybridization, and unique histories of lineages within the H. eximia group. These results suggest that the combination of nuclear, mitochondrial, and spatial information can provide a more complete picture of ‘how evolutionary history played out’, particularly in cases where mitochondrial introgression is known to occur. © 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116. doi: 10.1111/zoj.12169

ADDITIONAL KEYWORDS: Hylidae – Mexico – mitochondrial capture – phylogeography – speciation.

INTRODUCTION Coupling processes of diversification with spatial patterns of genetic structure is a central goal of phylogeographical research. However, achieving this goal can be complicated when gene histories are in conflict with each other. Two common reasons for

*Corresponding author. E-mail: [email protected]

genealogical discordance amongst loci are introgression and incomplete lineage sorting (Moore, 1995; Funk & Omland, 2003). Recent methodological advances using coalescent-based models can incorporate stochastic lineage sorting into phylogenetic estimation (Carstens & Knowles, 2007; Edwards, Liu & Pearl, 2007; Heled & Drummond, 2010). However, introgression, particularly deep in time, remains a confounding factor in phylogeographical research (Joly, McLenachan & Lockhart, 2009).

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

103

104

R. W. BRYSON JR ET AL.

Introgression involves the integration of DNA from one species into another species through the backcrossing of hybrids with pure parental types of one or both species. Introgression can cause DNA sequences from one species to be more similar to sequences of a different species and lead to incongruence amongst gene trees. Introgression can be relatively easy to distinguish after recent hybridization, which results in nonmonophyletic species in gene trees and shared sequences in areas of sympatry (Joly et al., 2009). Ancient introgression, however, may be more difficult to detect. As time since hybridization increases, introgressed sequences are more likely to sort and the geographical signals of introgression are more likely to erode (Funk & Omland, 2003; Linnen & Farrell, 2007). Identifying ancient introgression events becomes challenging because genetic lineages will have become reciprocally monophyletic in formerly hybridizing species and direct evidence for donor taxa is erased (Linnen & Farrell, 2007). The maternally inherited mitochondrial genome often introgresses more rapidly than biparentally or paternally inherited components of the nuclear genome (Ballard & Whitlock, 2004; Chan & Levin, 2005; but see Haldane’s rule, e.g. Carling & Brumfield, 2008). The smaller effective population size of mitochondrial DNA (mtDNA) compared with nuclear loci may facilitate this process, such that even low levels of introgression may be sufficient to establish an introgressed mitochondrial haplotype in a foreign population (Takahata & Slatkin, 1984; Ballard & Whitlock, 2004; Chan & Levin, 2005). This process can lead to the complete replacement of the mtDNA in the recipient species by mitochondrial genome capture (e.g. Rabosky et al., 2009; Reid, Demboski & Sullivan, 2012; Tang et al., 2012; Willis, Farias & Ortí, 2014). Rapid genetic drift in small populations and adaptive introgression can both accelerate the rate of mitochondrial capture (Ballard & Whitlock, 2004; Toews & Brelsford, 2012). Ancient hybridizations and mitochondrial gene capture followed by reproductive isolation may result in the establishment of a new evolutionary trajectory in the introgressed maternal genome over time and the development of new lineage-specific mutations (Bryson et al., 2010). As a result, the geographical distribution of lineages with a captured mtDNA genome may also be informative. For example, the horned lizards Phrynosoma platyrhinos and Phrynosoma goodei appear to have the captured mitochondrial genome of Phrynosoma mcalli (Leaché & McGuire, 2006), yet several studies inferred phylogeographical patterns in these species using only mtDNA (Jones, 1995; Mulcahy et al., 2006; Luxbacher & Knouft, 2009; Jezkova, 2010). If hybridization and mtDNA capture was followed by reproductive isolation, then mutations in the mtDNA of P. platyrhinos and P. goodei should develop inde-

pendent of P. mcalli. Through time this would result in genetic structuring in the mitochondrial genomes of P. platyrhinos and P. goodei that reflect their evolutionary history postdating introgression. Indeed, marked genetic structuring in the mtDNA of P. platyrhinos and P. goodei corresponds well to historical processes across the distribution of these lizards (Mulcahy et al., 2006). Phylogeographical inference from introgressed mtDNA has provided similar insight into the evolutionary histories of a number of other taxa, including fish (Willis et al., 2007; Gross, 2012), mammals (Shields et al., 2000; Good et al., 2008), and invertebrates (Zakharov et al., 2009; Marková et al., 2013). Mitochondrial introgression and gene capture have been documented within Nearctic hylid treefrogs in the Hyla eximia group (Bryson et al., 2010; Klymus et al., 2010; Klymus & Gerhardt, 2012). The H. eximia group is broadly distributed across the moderate to high elevations of south-western North America (Fig. 1), ranging from southern Utah and Colorado south across most of mainland Mexico and into Guatemala (Duellman, 2001). The taxonomic composition of this group has changed several times over the past 50 years (Duellman, 1970, 2001; Hua et al., 2009), but the core group is comprised of Hyla arenicolor, an undescribed species from the Balsas Basin of Mexico (Bryson et al., 2010; Klymus & Gerhardt, 2012), and the morphologically and ecologically similar species Hyla arboricola, Hyla euphorbiacea, Hyla eximia, Hyla plicata, Hyla walkeri, and Hyla wrightorum (Eliosa-León, 2002a, b; Smith et al., 2007; Hua et al., 2009; Bryson et al., 2010; Klymus & Gerhardt, 2012). Although H. arenicolor and the Balsas Basin species occupy rocky canyons and streams, the remaining species inhabit similar highland habitats and are generally associated with mixed pine−oak forests. In this study, we expanded the investigation of mtDNA introgression within the H. eximia group across a larger geography, more populations, and more species to better understand its role in revealing the evolutionary history and phylogeography of these treefrogs. We expected that the introgression of a mitochondrial genome into another species should not only result in mitochondrial/nuclear DNA discordance but also restart the progression of mtDNA evolutionary divergence from that moment onward. As a result, we predict that incorporating congruent and incongruent locus-specific phylogenetic estimates with the geographical distribution of lineages will provide valuable insight into evolutionary history. We tested this expectation by reconstructing and comparing the matrilineal phylogeny of the H. eximia group with the phylogenetic history inferred from nuclear genes. We confirmed introgression between lineages using a posterior predictive checking model and data from both the mtDNA and nuclear DNA (nuDNA).

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

INTROGRESSION IN NEARCTIC TREEFROGS 110° W

105° W

100° W

95° W

110° W

40° N

40° N

35° N

35° N

H. arenicolor

105° W

100° W

105

95° W

H. wrightorum

30° N

30° N

25° N

25° N

H. plicata 20° N

H. eximia

20° N

N

N

Balsas Basin

15° N

H. arboricola H. euphorbiacea

15° N 0 125 250

500

750

0 125 250

1000 km

500

750

1000 km

H. walkeri

Figure 1. Distribution of Hyla eximia group treefrogs (thick black lines). Black dots indicate localities sampled in this study (listed in Table S1). The Balsas Basin lineage represents an undescribed species within the group (Bryson et al., 2010; Klymus & Gerhardt, 2012).

MATERIAL AND METHODS TAXON

SAMPLING AND LABORATORY METHODS

We collected tissue samples from 128 frogs of the currently recognized species in the H. eximia group from throughout their distributions (Fig. 1; Table S1). We also included samples from H. arboricola, recently removed from the synonymy of H. eximia (Eliosa-León, 2002a, b). We did not include samples of Hyla bocourti, a possible member of the H. eximia group (Duellman, 2001) that may be extinct (Acevedo & Young, 2004). We used Hyla cinerea and Hyla versicolor as outgroups (Pyron & Wiens, 2011). We augmented sampling with tissues obtained from museum collections, indicated in Table S1. Following methods specified in Bryson et al. (2010), we extracted genomic DNA from toe or liver tissue and sequenced three mitochondrial DNA gene regions, including part of the NADH dehydrogenase subunit 4 and its flanking tRNAs (ND4, 759 bp) and the complete ATPase subunits 8 and 6 (ATPase8, 165 bp; ATPase6, 684 bp; with an overlap of 10 bp between the two genes). We also sequenced three nuclear genes, including 520 bp of proopiomelanocortin A (POMC), 169 bp of β-crystallin (cryB), and 517 bp of tyrosinase (Tyr1). Primer sequences for Tyr1 were from Bossuyt & Milinkovitch (2000). We edited and manually aligned forward

and reverse sequences for each individual using SEQUENCHER v. 4.2 (Gene Codes Corporation, Ann Arbor, MI). For cryB data, which contained numerous indels, we performed an additional sequence alignment with MAFFT v. 6 (Katoh et al., 2002; Katoh & Toh, 2008) using default settings and the G-INS-i algorithm. We computationally determined the gametic phase of heterozygous variants using PHASE v. 2.1.1 (Stephens & Donnelly, 2003). For each nuDNA data set, separate runs of 400 iterations each were carried out, accepting results with a probability threshold of 0.7 or higher. All polymorphic sites with a probability less than 0.70 were coded in both alleles with the appropriate IUPAC ambiguity code. To test for recombination within nuclear loci, we performed three independent recombination tests (RDP, GENECOV, and MaxChi) in the program RDP v. 3.44 (Martin et al., 2010) using the default settings for all three methods.

PHYLOGEOGRAPHICAL

ESTIMATION

We generated a mtDNA phylogeny using MrBayes v. 3.1 (Ronquist & Huelsenbeck, 2003) to examine geographical structure within the H. eximia group. Geographical lineages within H. arenicolor inferred in a previous study from either mtDNA or amplified fragment length polymorphism data were nearly

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

106

R. W. BRYSON JR ET AL.

identical (Klymus & Gerhardt, 2012), which suggests that mtDNA is a useful marker for broadly delimiting geographical lineages within the H. eximia group. We used MrModeltest v. 2.1 (Nylander, 2004) to select the best-fit model of evolution for each gene region, based on Akaike information criteria. Analyses in MrBayes consisted of four runs (nruns = 4) conducted with three heated and one cold Markov chain sampling every 100 generations for 5 000 000 generations. We used a mean branch length exponential prior of 100 and heating temperature of 0.02, which improved mixing and convergence. We verified that the average standard deviation of split frequencies was below 0.01 at the end of the analysis and visualized output parameters using TRACER v. 1.4 (Rambaut & Drummond, 2007) to further ascertain convergence and stationarity. We discarded all samples obtained during the first 1 000 000 (25%) generations as burn-in. We used a general mixed Yule coalescent model (GMYC) to delineate geographically distinct clusters of samples (hereafter referred to as ‘lineages’). The GMYC model identifies independent evolutionary clusters by detecting a threshold value at the transition from interspecific to intraspecific branching patterns. Although this model has been recently criticized for potentially over-splitting species with pronounced genetic structure (Satler, Carstens & Hedin, 2013), it remains an effective way to estimate genetically distinct geographical clusters at the phylogeographical (intraspecific) scale. We used SPLITS (Pons et al., 2006) in the R command-line package (R Development Core Team, 2011) to run a likelihood implementation of the GMYC model. We performed both single and multiple threshold versions of the model. The multiple threshold variant tests whether allowing the species-coalescent transition to vary across the tree significantly improves the fit to the model (Monaghan et al., 2009). We generated an ultrametric tree for use in SPLITS using BEAST v. 1.7.2 (Drummond & Rambaut, 2007). We used separate models of evolution for each of the four gene regions and we ran the analyses for 40 000 000 generations, with samples retained every 1000 generations, using a Yule tree prior. Results were displayed in TRACER to confirm acceptable mixing and likelihood stationarity, appropriate burn-in, and adequate effective sample sizes above 200 for all estimated parameters. After discarding the first 4 000 000 generations (10%) as burn-in, the parameter values of the samples from the posterior distribution were summarized on the maximum clade credibility tree using TreeAnnotator v. 1.7.2 (Drummond & Rambaut, 2007).

SPECIES

TREES AND DIVERGENCE TIMES

We used *BEAST (Heled & Drummond, 2010), a part of the BEAST v. 1.6.2 package (Drummond & Rambaut,

2007), to estimate time-calibrated species trees. Because the mtDNA data were suspected to support a different topology as a result of mitochondrial introgression and gene capture (Bryson et al., 2010; Klymus et al., 2010; Klymus & Gerhardt, 2012), we analysed the mtDNA and nuDNA data sets separately. In our mtDNA species-tree reconstructions, we included all samples (N = 123) and the two outgroups with the exception of samples of H. arenicolor from the Grand Canyon of Arizona, which contain introgressed mtDNA indistinguishable from H. wrightorum (Bryson et al., 2010; Klymus et al., 2010; Klymus & Gerhardt, 2012). In our nuDNA species-tree analyses, we included two to four individuals (N = 35) from each phylogeographical lineage inferred from our mtDNA data and the two outgroups. We estimated best-fit models of evolution for each gene using MrModeltest. In each species-tree analysis we used relaxed uncorrelated lognormal clocks for each gene and a Yule speciation prior. For mtDNA analyses, we unlinked substitution and clock models for each gene region. We ran analyses for 100 000 000 generations, with samples retained every 1000 generations. Results were displayed in TRACER to confirm acceptable mixing and likelihood stationarity, appropriate burnin, and adequate effective sample sizes. After discarding the first 10% of generations as burn-in, the parameter values of the samples from the posterior distribution were summarized on the maximum clade credibility tree using TreeAnnotator. We repeated this burn-in and visualization procedure for each of the three nuclear gene trees coestimated by *BEAST as part of the nuDNA species-tree reconstruction. To calibrate our mtDNA species tree, we used the mean rate estimated from the fossil-calibrated chronogram of Hyla in Bryson et al. (2010). This mean rate, estimated for the same mtDNA genes used here, was 9.63 × 10−3 substitutions/site/million years (1.926% change between lineages per million years). We gave the mtDNA rate a lognormal distribution and specified a logarithmic standard deviation of 0.2 around this mean to produce a wide upper and lower 95% prior credibility interval of 6.79 × 10−3 and 1.312 × 10−2, similar to the 95% posterior credibility interval in Bryson et al. (2010). To calibrate our nuDNA species tree, we estimated clock rates relative to a mtDNA rate. Using MEGA5 (Tamura et al., 2011), we calculated mean betweengroup uncorrected p-distances from mtDNA data and each of the three nuclear genes for three paired groupings in our study: H. arenicolor from the Sonoran Desert/ H. arenicolor from the Chihuahuan Desert + Central Mexican Plateau, H. plicata/H. euphorbiacea, and the Balsas Basin species/H. eximia from the southern Sierra Madre Occidental + Central Mexican Plateau. For each nuclear gene, genetic distances between each paired grouping were divided by the corresponding mtDNA

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

INTROGRESSION IN NEARCTIC TREEFROGS distances, and the average of these three numbers calculated to obtain a comparative mutation rate relative to mtDNA. These relative rates were multiplied by the mtDNA mutation rate from above (9.63 × 10−3). Using this method, we obtained the following rates for substitutions/site/million years: cryB = 4.43 × 10−3, POMC = 2.99 × 10−3, and Tyr1 = 1.35 × 10−3. The calculated rate for the Tyr1 gene was similar to the Tyr1 rate estimated for the Hyloidae superfamily of frogs (which includes Hyla) using fossil data (1.95 × 10−3, Fouquet et al., 2012). This similarity is important because it suggests that our divergence dates may still be informative despite a history of introgression in the H. eximia group (see below). In fact, reconstructing the nuclear species tree with the Tyr1 substitution rate from Fouquet et al. (2012), with rates for cryB and POMC estimated relative to Tyr1 by *BEAST, yields generally similar divergence dates across the tree (e.g. dates are within 2 000 000 years of each other at the base of the tree, within 1 000 000 years at the base of the H. eximia group, and identical at the base of the H. arenicolor clade). As with our mtDNA species-tree analyses, we used a lognormal distribution and a relatively wide logarithmic standard deviation of 0.2 for the prior distribution of each nuclear gene’s mutation rate.

lated data sets to determine how well the model fits the data. If the observed pairwise distance is smaller than 95% of the simulated value, then a model involving no hybridization can be rejected and it can be concluded that the inaccuracy of the model is because of the presence of hybrid sequences (Joly et al., 2009). We conducted two separate analyses. First, we inferred a species tree in *BEAST using both mtDNA and nuDNA and the same priors as specified above. We used the posterior distribution of species trees, population sizes, and branch lengths to simulate gene trees and sequence data sets in JML under the coalescent with no migration. We then assessed whether introgression could be detected for any pair of individuals in any locus at a probability of < 0.05. However, this species-tree estimate may be skewed as a result of conflict between the mtDNA and nuclear gene trees owing to mtDNA introgression (e.g. Eckert & Carstens, 2008). Accordingly, we repeated the analysis using a species tree reconstructed in *BEAST using only the nuDNA, and followed the assumption that this second analysis should result in a more accurate estimate of introgression within the nuclear loci.

RESULTS GENETIC

INTROGRESSION Although previous studies (Bryson et al., 2010; Klymus et al., 2010; Klymus & Gerhardt, 2012) have suggested mtDNA introgression between H. arenicolor and other species of the H. eximia group, this hypothesis has not been formally tested using a statistical framework. Further, limited evidence suggests that nuclear introgression between H. arenicolor and the H. eximia group might also be possible (Bryson et al., 2010). As our species-tree analyses assume that the only process creating discordance amongst loci is genetic drift within species (Liu & Pearl, 2007), incorporating gene trees that have topologies resulting from introgression could alter the parameter estimates for the model and skew the species-tree topology (Willis et al., 2014). We used a posterior predictive checking approach (Joly et al., 2009) implemented in JML v. 1.01 (Joly, 2012) to test for introgression between lineages for all loci. This program tests whether the minimum genetic distance between the sequences of two individuals is smaller than those simulated under a model involving no hybridization such as incomplete lineage sorting (Joly et al., 2009). JML uses as input the posterior of species trees with branch lengths and population sizes from *BEAST to generate a gene tree from which DNA sequences are then simulated. A test quantity, the minimum pairwise distance between sequences of two species, is estimated from the observed and simu-

107

DATA

The full mtDNA data set (1598 bp) contained 482 parsimony-informative sites. The nuclear loci contained less variation than the mtDNA (parsimonyinformative sites, POMC: 76/520 bp; cryB: 45/169 bp; Tyr1: 63/517 bp). We were unable to obtain ATPase sequence data for two individuals (MX435 H. walkeri and MX436 H. eximia). Models of sequence evolution selected for gene partitions were GTR + I + G (ND4, ATPase6) and HKY + G (ATPase8, tRNAs) in the mtDNA gene tree data set, and GTR + I + G (ND4, ATPase6, POMC), HKY + G (ATPase8, tRNAs), GTR + I (cryB), and HKY + I (Tyr1) for the species-tree data sets. None of the three methods used in RDP detected recombination in the nuclear genes. All aligned sequences were deposited in the Dryad repository: http://dx.doi.org/10.5061/dryad.3rc02.

PHYLOGEOGRAPHICAL

ESTIMATION

We identified two major mtDNA clades within the H. eximia group (Fig. 2). One clade contained most of the H. arenicolor samples. The second clade contained the Balsas Basin species, H. arenicolor from the south-western region of the Central Mexican Plateau, H. arenicolor from the Grand Canyon of Arizona, and the six morphologically and ecologically similar species H. eximia, H. wrightorum, H. plicata, H. euphorbiacea, H. arboricola, and H. walkeri.

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

108

R. W. BRYSON JR ET AL. MX135 Hyla cinerea MX136 Hyla versicolor

MX294 MX177 MX371 MX260 MX149 MX293 MX368 MX347 MX414 MX415 MX292 MX159 MX370 MX150 MX223 MX157 MX248 MX178 MX329 MX330 MX384 MX158 MX376 MX378 MX361 MX222 MX247 MX261 MX354 MX155 MX249 MX257 MX362 MX156 MX251 MX256 MX151 MX250 MX258 MX259 MX152 MX252

arenicolor SD arenicolor CP

arenicolor CD

.73

arenicolor CMP

MX154 MX254 MX359 MX348 MX358 MX255 MX349 MX138 MX243 MX436 MX122 MX147 MX31 eximia MX126 MX340 MX39 MX121 MX124 MX146 MX387 arenicolor MX388 MX389 Grand Canyon MX390 MX32 MX123 MX125 MX341 MX360 MX153 MX176 MX253 MX357 MX246 MX355 MX356 MX234 MX174 MX139 MX33 MX140 MX132 eximia MX141 MX173 eximia MX134 MX142 MX353 MX143 euphorbiacea MX343 euphorbiacea MX344 MX235 MX236 MX237 MX70 MX71 MX144 MX145 MX326 MX435 MX232 MX233 MX29 MX148 MX130 MX175 MX241 MX38 MX405 MX40 MX129 MX37 MX34 MX128 MX65 MX170 MX242 MX403 MX131 MX172 MX245 MX345 MX35 MX240 MX68 MX239 MX127 MX238

Balsas eximia DGO

wrightorum

0.1 sub./ site

.66

arenicolor CMPsw

plicata

.72

.91

.66

euphorbiacea W euphorbiacea E walkeri arboricola eximia OrS

eximia OcS

eximia CMP

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

INTROGRESSION IN NEARCTIC TREEFROGS

109

Figure 2. Phylogeographical structure within Hyla eximia group treefrogs inferred from Bayesian analyses of 1598 bp of mitochondrial DNA. Lineages were inferred using a single-threshold general mixed Yule coalescent model. All major nodes that received ≥ 0.95 Bayesian posterior probability are depicted with black dots. Bayesian posterior probability support values at all other nodes are provided. Abbreviations: CD, Chihuahuan Desert; CMP, Central Mexican Plateau; CMPsw, south-western Central Mexican Plateau; CP, Colorado Plateau; DGO, Durango; E, eastern; OcS, southern Sierra Madre Occidental; OrS, southern Sierra Madre Oriental; SD, Sonoran Desert; W, western. ◀

110° W

105° W

100° W

95° W

40° N

110° W

105° W

100° W

95° W

40° N

Balsas Basin

ex

ex

arenicolor CMPsw arenicolor SD

35° N

35° N

arenicolor CP arenicolor CD

eu

arenicolor CMP 30° N

30° N

25° N

25° N

ex

wrightorum arboricola plicata euphorbiacea E euphorbiacea W

20° N

walkeri

20° N

eximia OrS eximia DGO

N

N

eximia OcS eximia CMP

15° N

15° N 0 125 250

500

750

1000 km

0 125 250

500

750

1000 km

Figure 3. Geographical distribution of inferred mitochondrial lineages within Hyla eximia group treefrogs. Lineages are colour-coded to correspond to haplotypes plotted on the map. Circled dots indicate localities of samples used in multilocus species-tree reconstructions. Five mismatches between species designation (H. eximia, abbreviated ex; Hyla euphorbiacea, abbreviated eu) and mitochondrial lineage are noted. Abbreviations: CD, Chihuahuan Desert; CMP, Central Mexican Plateau; CMPsw, south-western Central Mexican Plateau; CP, Colorado Plateau; DGO, Durango; E, eastern; OcS, southern Sierra Madre Occidental; OrS, southern Sierra Madre Oriental; SD, Sonoran Desert; W, western.

The single and multiple threshold GMYC models were not significantly different from each other (X2 = 6.315, df = 3, P = 0.097) so results from the single threshold model were used to infer phylogeographical lineages. The 16 independent evolutionary clusters estimated by the GMYC model were generally concordant with taxonomy and geography (Figs 2, 3). Hyla wrightorum, H. plicata, H. walkeri, and H. arboricola each formed single lineages. One haplotype of H. eximia (MX31) from Durango was nested within the H. wrightorum lineage (Fig. 2). Two haplotypes of H. euphorbiacea and two haplotypes of H. eximia grouped with H. plicata haplotypes within the H. plicata lineage (Fig. 2). All haplotypes of H. arenicolor from the Grand Canyon of Arizona grouped with H. wrightorum, consistent with previous studies (Klymus et al., 2010; Klymus & Gerhardt, 2012). Pronounced genetic structure was detected in H. arenicolor, H. eximia, and H. euphorbiacea,

and each species contained several distinct geographically delimited lineages (Figs 2, 3). Five lineages of H. arenicolor were identified from well-defined biogeographical regions, including the Sonoran Desert and associated highlands (arenicolor SD), Chihuahuan Desert and associated highlands (arenicolor CD), Colorado Plateau (arenicolor CP), Central Mexican Plateau (arenicolor CMP), and south-western region of the Central Mexican Plateau (arenicolor CMPsw). Haplotypes of H. arenicolor from the Grand Canyon of Arizona were indistinguishable from H. wrightorum. Hyla eximia formed four lineages distributed in southwestern Durango (eximia DGO), the southern Sierra Madre Oriental (eximia OrS), across the Central Mexican Plateau (eximia CMP), and in the southernmost Sierra Madre Occidental and two sites to the south (Tapalpa, Jalisco) and south-east (Nuevo San Joaquin, Querétaro; eximia OcS). This Nuevo San

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

110

R. W. BRYSON JR ET AL.

Joaquin, Querétaro locality was well within the distribution of the Central Mexican Plateau lineage (Fig. 3); otherwise, the distributions of the four lineages of H. eximia were non-overlapping. The two lineages of H. euphorbiacea corresponded to regions of the Sierra Madre del Sur in eastern (euphorbiacea E) and western (euphorbiacea W) Oaxaca.

SPECIES

TREES AND DIVERGENCE TIMES

Reconstructions of the species tree based on the mtDNA and nuDNA data were strongly discordant (Fig. 4). Phylogenetic relationships amongst the lineages varied between the two species trees to such an extent that relationships between only four of the 16 lineages were similar in both trees (the sister relationship between H. arenicolor from the Chihuahuan Desert and Central Mexican Plateau, and the sister relationship between H. eximia from the southern Sierra Madre Occidental and the Central Mexican Plateau). Although nodal support across both species trees was generally weak, even the well-supported nodes in each tree supported different relationships. Estimated divergence dates in the mtDNA species tree were noticeably older and 95% posterior credibility intervals considerably wider compared with the nuDNA species tree. For example, the basal divergence within the H. eximia group was estimated to have occurred 15.8 Mya (24.3–8.0 Mya) in the mtDNA species tree. However, the estimated date of this divergence event in the nuDNA species tree was many millions of years later, around 6.4 Mya (8.8– 4.2 Mya). Based on mean estimated dates of divergence, 47% of the divergences within the H. eximia group in the mtDNA species tree occurred during the Pleistocene as compared to 67% in the nuDNA species tree. Each of the three nuclear gene trees contained limited phylogenetic resolution and relatively few strongly supported nodes, especially towards the base of the trees (Fig. S1). However, samples of several lineages consistently grouped together, such as H. walkeri, H. wrightorum, and Balsas Basin. Lineages of H. arenicolor formed a monophyletic group in the Tyr1 and cryB gene trees, and samples of H. arboricola formed a monophyletic group in the POMC and cryB gene trees.

INTROGRESSION Results from JML suggested that introgression between lineages was widespread in the mtDNA data set. The analysis detected 145 instances of potential introgression involving samples from all lineages except the Balsas Basin and H. arenicolor from the Sonoran and Chihuahuan Deserts, Colorado Plateau, and Central Mexican Plateau (Table S2). The nuDNA data set, however, contained relatively few occurrences of potential introgression. Based on analyses using the species

tree reconstructed from both mtDNA and nuDNA, there was no strong evidence for introgression in the POMC or Tyr1 data sets. Introgression in the cryB data set was detected between H. arboricola and H. euphorbiacea sequences. JML results inferred from the species tree estimated using only the nuDNA suggested no significant introgression in the Tyr1 and cryB and only four occurrences in the POMC data set (involving one sequence of H. eximia from the southern Sierra Madre Oriental and one sequence of H. plicata, and one sequence of H. eximia from the southern Sierra Madre Oriental and three sequences of H. arboricola).

DISCUSSION H. EXIMIA

MTDNA INTROGRESSION IN THE

GROUP

Although previous studies have provided evidence for mitochondrial introgression and gene capture between H. arenicolor and H. eximia, and H. arenicolor and H. wrightorum (Bryson et al., 2010; Klymus et al., 2010; Klymus & Gerhardt, 2012), results from our study suggest that these phenomena are widespread amongst nearly all species in the H. eximia group. Within the six morphologically and ecologically similar species H. eximia, H. wrightorum, H. plicata, H. euphorbiacea, H. arboricola, and H. walkeri, two lineages in particular appear to have captured the mtDNA genome of a different lineage. Three H. eximia from three different localities in south-western Durango form a genetically distinct monophyletic lineage that is strongly supported (1.0 posterior probability, PP) as sister to H. wrightorum based on analyses of mtDNA (Figs 2, 4). However, analyses of the nuDNA alternatively suggest that all lineages of H. eximia form a monophyletic group (0.91 PP), consistent with call and morphological data (Duellman, 2001). Two samples of H. arboricola from the Sierra Madre del Sur in Guerrero are placed as sister to H. eximia from the southern Sierra Madre Oriental in the mtDNA gene tree (1.0 PP, Fig. 2) and mtDNA species tree (0.72 PP, Fig. 4). However, the nuDNA species tree strongly supports (1.0 PP) the placement of H. arboricola in a clade with two geographically proximate treefrogs, H. plicata from the Trans-Mexican Volcanic Belt and H. euphorbiacea from the Sierra Madre del Sur (Figs 3, 4). The gene histories of the mtDNA and nuDNA in H. eximia group treefrogs clearly appear to be in conflict.

EVOLUTIONARY

HISTORY OF THE

H. EXIMIA

GROUP

Our time-calibrated species tree reconstructed from nuDNA (Fig. 4) provides a reasonable estimate of the evolutionary history of the H. eximia group as introgression within the nuclear genes appears to be low. Unfortunately, many of the nodes across this tree are not well supported, perhaps because of the reduced

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

INTROGRESSION IN NEARCTIC TREEFROGS

111

cinerea versicolor arenicolor SD arenicolor CP .72

arenicolor CD arenicolor CMP

.90 .89

Balsas Basin wrightorum eximia DGO arenicolor CMPsw

.58

plicata euphorbiacea E

.73

euphorbiacea W

.55

walkeri

.77

.69

eximia OcS eximia CMP

.89

mtDNA

arboricola

.62

eximia OrS

.72

25

20

15

10

5

0

cinerea versicolor Balsas Basin wrightorum

Hyla arenicolor

.86

arenicolor CMPsw .57

Balsas Basin

arenicolor SD arenicolor CP

<.50 <.50 .78

arenicolor CD arenicolor CMP arboricola

Hyla wrightorum

plicata .91

euphorbiacea W

.67

Hyla eximia

euphorbiacea E walkeri eximia OrS

.64

eximia DGO

.91

nuDNA

.84 .52

25 OLIGO

20

15 MIOCENE

10

eximia OcS eximia CMP

5 0 PLIO PLEIS

Time (millions of years ago) Figure 4. Species-tree reconstructions for Hyla eximia group treefrogs estimated from mitochondrial DNA (mtDNA) or nuclear DNA (nuDNA) using *BEAST. Bars indicate 95% posterior credibility intervals of divergence dates. All major nodes that received ≥ 0.95 Bayesian posterior probability are depicted with black dots. Bayesian posterior probability support values at all other nodes are provided. Inset shows oscillograms of the advertisement calls (90 s) of four H. eximia group species. Note call similarities between Hyla wrightorum and H. eximia and call differences between H. wrightorum and Hyla arenicolor (data from Klymus, Humfeld & Gerhardt, 2012). Abbreviations: CD, Chihuahuan Desert; CMP, Central Mexican Plateau; CMPsw, south-western Central Mexican Plateau; CP, Colorado Plateau; DGO, Durango; E, eastern; OcS, southern Sierra Madre Occidental; Oligo, Oligocene; OrS, southern Sierra Madre Oriental; Pleis, Pleistocene; Plio, Pliocene; SD, Sonoran Desert; W, western.

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

112

R. W. BRYSON JR ET AL.

phylogenetic signal in the three nuclear genes (Fig. S1). As a result it is difficult to confidently discern the relationships of many lineages. Nonetheless, several patterns within the species tree are evident. First, three clades appear to be reasonably supported (> 0.90 PP): H. eximia; H. arenicolor; and H. arboricola + H. euphorbiacea + H. plicata. The compositions of the former two clades match current taxonomic designations based on mating calls and morphology (Duellman, 2001; Klymus et al., 2010; Klymus & Gerhardt, 2012). The latter clade is a composite of three species distributed in southern Mexico (Fig. 3) across three mountain ranges with shared biotic affinities (Paniagua & Morrone, 2009). Second, speciation within each of these three clades occurred nearly entirely within the Pleistocene. Estimated divergence dates between lineages of H. arenicolor are dramatically younger than those based on an mtDNA gene-tree estimate (Bryson et al., 2010). As a result, diversification in these and other H. eximia group treefrogs was probably more heavily impacted by Pleistocene glacial–interglacial cycles than older Neogene events. Lastly, the inclusion of H. wrightorum with H. arenicolor (0.86 PP) and the basal placement of the Balsas Basin lineage within the H. eximia group (0.78 PP) are both contrary to relationships based on call and morphology data (Klymus & Gerhardt, 2012, and references therein). These relationships suggest that either morphology and mating calls are plastic or that a more complex history of hybridization within the H. eximia group exists that was undetected by our sampling and JML analyses. Nearly identical mating calls and similar morphologies between H. wrightorum and H. eximia (Klymus & Gerhardt, 2012; Fig. 4), both inferred as nonsister species in our species tree, suggest that the latter explanation may be true. Further, the lack of significant evidence for hybridization inferred from our JML analyses simply suggests that the program could not distinguish whether any incongruence was the result of introgression or incomplete lineage sorting, which is not evidence that hybridization was completely absent (Joly et al., 2009). Future studies with more nuclear loci may reveal an additional layer of complexity within the evolutionary history of the H. eximia group.

IS

THE MTDNA GENOME IN THE

H. EXIMIA

GROUP STILL INFORMATIVE?

Although it seems clear that mtDNA introgression and gene capture severely limit the utility of this marker for reconstructing phylogenetic relationships within the H. eximia group, can mtDNA still provide insight into the evolutionary history of these widespread treefrogs? At a minimum, the general concordance between the geographical distribution of mitochondrial lineages and biogeographical regions suggests that mtDNA

still retains important information about the phylogeographical structure of these treefrogs across south-western North America. Previous studies have demonstrated that at least four of the five maternal lineages of H. arenicolor and the Balsas Basin species inhabit distinct biogeographical areas (Barber, 1999; Bryson et al., 2010; Klymus et al., 2010; Klymus & Gerhardt, 2012), including the Sonoran and Chihuahuan Deserts, Colorado Plateau, Central Mexican Plateau, and Balsas Basin. One lineage of H. arenicolor is distributed across the south-western corner of the Central Mexican Plateau, a region also inhabited by a lineage of H. eximia (Fig. 3) and a lineage of Mexican gophersnake (Bryson, García-Vázquez & Riddle, 2011a). Most of the remaining nine lineages within the H. eximia group are also allopatrically distributed within distinct biogeographical regions (Fig. 3). Hyla wrightorum is distributed broadly across the Sierra Madre Occidental and adjacent northern highlands, H. plicata is distributed almost exclusively within the higher elevations of the Trans-Mexican Volcanic Belt, and H. walkeri is found along the Chiapan Highlands. Hyla arboricola occurs in the western Sierra Madre del Sur of Guerrero, a region characterized by high biotic endemism (Fa & Morales, 1993). Of the three remaining lineages of H. eximia, one inhabits the Central Mexican Plateau and one is distributed along the southern Sierra Madre Oriental, a highland region that contains other distinct phylogeographical lineages (Durish et al., 2004; Bryson et al., 2011b, c). The third lineage of H. eximia occurs in the humid pine−oak forests of south-western Durango. Although this area lies within the Sierra Madre Occidental, the wetter climate on the Pacific versant supports a humid pine−oak forest that differs from the drier, cooler pine−oak forests along the interior (Webb, 1984) and contains several narrow endemics (Tanner, Dixon & Harris, 1972; Bryson et al., 2007). Lastly, two lineages of H. euphorbiacea are distributed in regions of the Sierra Madre del Sur in eastern and western Oaxaca. Similar genetic lineages have been found in a number of highland taxa (Parra-Olea et al., 2005; Puebla-Olivares et al., 2008; Ruiz-Sanchez & Specht, 2013). The mitochondrial genome of H. eximia group treefrogs may also point to ancient areas of hybridization and episodes of gene capture from extinct ancestors. For example, results from JML suggest that H. arboricola from the Sierra Madre del Sur in Guerrero and H. eximia from the southern Sierra Madre Oriental have a history of hybridization. However, these two areas are currently isolated from each other by the arid Balsas Basin drainage (Fig. 3). Hybridization provides evidence for a recent connection of these regions by ephemeral corridors of pine−oak woodlands during Pleistocene glacial periods (Bryson et al., 2011b). Similarly, H. arenicolor from the Grand Canyon of Arizona and

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

INTROGRESSION IN NEARCTIC TREEFROGS H. wrightorum from the pine forests of northern Arizona are geographically isolated yet share a similar mitochondrial genome, suggesting hybridization occurred in this region during the Pleistocene (Bryson et al., 2010). The deep mtDNA divergence of the Sonoran Desert lineage of H. arenicolor from other H. arenicolor (Figs 2, 4), a pattern not recovered with nuDNA (Fig. 4), provides evidence for a unique history of this lineage. One possible explanation is that the Sonoran Desert lineage possesses the introgressed mitochondrial genome of a divergent extinct lineage of H. arenicolor (Klymus & Gerhardt, 2012) and thus represents an ‘mtDNA fossil’ (Bossu & Near, 2009). This may also be the case for the Balsas Basin lineage. Such a scenario has been posited previously for hylid frogs (Holloway et al., 2006; Klymus & Gerhardt, 2012) and fish (Bossu & Near, 2009; Nevado et al., 2011), and highlights the usefulness of introgressed mitochondrial genomes when inferring evolutionary history. Although important information on the evolutionary history of lineages can be retained in mitochondrial genomes impacted by introgression and gene capture, it is worth noting that this utility may be influenced by the causes of introgression. If introgression is strongly driven by selection (e.g. through sex-biased hybridization or adaptive selection on mtDNA), then inferences may be different than if it is purely neutral (see Toews & Brelsford, 2012). Although the causes of mitochondrial introgression in the H. eximia group are unknown, they may be a result of sex-biased asymmetry (Klymus et al., 2010; Klymus & Gerhardt, 2012) or adaptive introgression (Bryson et al., 2010). As argued by Toews & Brelsford (2012), future research is needed to identify the drivers of mitochondrial introgression to further address questions about how demographics and natural selection might influence the evolutionary history of species groups prone to introgression.

CONCLUSIONS This study adds to a growing list of species groups that demonstrate phylogenetic discord resulting from mtDNA introgression and capture (e.g. Bossu & Near, 2009; Leaché, 2010; Reid et al., 2012; Ruane et al., 2014). However, by incorporating the congruent and incongruent locus-specific phylogenetic estimates with the geographical distribution of mtDNA lineages, we gain valuable insight into phylogeography and evolutionary processes important to speciation. It is because of the very fact that introgression is usually limited in time and space to discrete geographical regions that the phylogeography is amenable to reconstruction through careful detective work on introgressed mtDNA genomes. We therefore suggest that including mtDNA in studies of species groups prone to introgression can sometimes provide a more complete picture of ‘how evo-

113

lutionary history played out’ – complete with divergence and bouts of gene flow – rather than dismissing mtDNA altogether simply because it can return ‘incorrect’ phylogenies. These ‘incorrect’ phylogenies turn out to be more reflective of incorrect assumptions on the part of researchers (see Rubinoff & Holland, 2005; Edwards & Bensch, 2009) than they are of fatal problems with the underlying data.

ACKNOWLEDGEMENTS We dedicate this study to the late Fernando MendozaQuijano. Without his enthusiasm and generous support through the years, this study would not have been possible. For providing tissues, we thank the following: Ambrose Monell Cryo Collection and D. R. Frost (AMNH), D. Dittmann and R. Brumfield (LSUMNS), C. Cicero (MVZ), D. Lazcano (Universidad Autónoma de Nuevo León), J. Alvarado-Díaz and A. Quijada-Mascareñas (Universidad Michoacana de San Nicolás de Hidalgo, Michoacán), J. A. Campbell, C. Franklin, and E. N. Smith (University of Texas at Arlington), P. Barber, H. EliosaLeón, C. Grunwald, R. Howe, J. Jones, J. Jungles, K. Klymus, F. R. Mendoza-Paz, the late F. MendozaQuijano, E. Mociño-Deloya, E. O’Neill, J. ReyesVelasco, and J. J. Wiens. For assistance in the field, we are grateful to E. Enderson, M. R. Graham, C. Grunwald, R. W. Hansen, C. Harrison, M. I. Ingrasci, J. Jones, G. Quijano-Manila, R. Mendoza-Paz, the late F. MendozaQuijano, S. Neiswenter, S. Ruane, M. Price, B. Thomason, M. Torocco, and G. Weatherman. Permission provided by D. M. Dennis to use frog drawings and K. Klymus to use oscillograms in Figure 4. This project was funded in part through grants from the AMNH (Theodore Roosevelt Memorial Fund), SWAN (Howard McCarley Student Research Award), UNLV (Graduate and Professional Student Association, and Major Research Instrumentation grant DBI-0421519), CONACyT grant no. 154093 to A. Nieto-Montes de Oca, and NSF grant DEB-0102383 to J. A. Campbell and O. Flores Villela. Permits to collect in Mexico were issued by SEMARNAT to RWB, UOGV, the late Fernando Mendoza-Quijano, and UTA and MZFC. All handling of animals followed animal use protocols approved by the University of Nevada at Las Vegas Animal Care Committee (R7011105-203). For additional support and advice, we thank J. Chaves, D. Lazcano, J. Klicka, K. Klymus, J. E. McCormack, and J. J. Wiens. Finally, we thank the Associate Editor and two anonymous reviewers for insightful comments on previous versions of this manuscript.

REFERENCES Acevedo M, Young B. 2004. Hyla bocourti. The IUCN Red List of Threatened Species, Version 2014.1. Available at: http:// www.iucnredlist.org.

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

114

R. W. BRYSON JR ET AL.

Ballard JWO, Whitlock MC. 2004. The incomplete naturalhistory of mitochondria. Molecular Ecology 13: 729– 744. Barber PH. 1999. Phylogeography of Hyla arenicolor (Cope) based on mitochondrial sequence data. Molecular Ecology 8: 547–562. Bossu CM, Near TJ. 2009. Gene trees reveal repeated instances of mitochondrial DNA introgression in orangethroat darters (Percidae: Etheostoma). Systematic Biology 58: 114– 129. Bossuyt F, Milinkovitch MC. 2000. Convergent adaptive radiations in Madagascan and Asian ranid frogs reveal covariation between larval and adult traits. Proceedings of the National Academy of Sciences, USA 97: 6585– 6590. Bryson RW, García-Vázquez UO, Riddle BR. 2011a. Phylogeography of Middle American gophersnakes: mixed responses to biogeographical barriers across the Mexican Transition Zone. Journal of Biogeography 38: 1570–1584. Bryson RW, Murphy RW, Graham MR, Lathrop A, Lazcano-Villareal D. 2011b. Ephemeral Pleistocene woodlands connect the dots for highland rattlesnakes of the Crotalus intermedius group. Journal of Biogeography 38: 2299– 2310. Bryson RW, Murphy RW, Lathrop A, Lazcano-Villareal D. 2011c. Evolutionary drivers of phylogeographical diversity in the highlands of Mexico: a case study of the Crotalus triseriatus species group of montane rattlesnakes. Journal of Biogeography 38: 697–710. Bryson RW, Nieto-Montes de Oca A, Jaeger JR, Riddle BR. 2010. Elucidation of cryptic diversity in a widespread Nearctic treefrog reveals episodes of mitochondrial gene capture as frogs diversified across a dynamic landscape. Evolution 64: 2315–2330. Bryson RW, Pastorini J, Burbrink FT, Forstner MRJ. 2007. A phylogeny of the Lampropeltis mexicana complex (Serpentes: Colubridae) based on mitochondrial DNA sequences suggests evidence for species-level polyphyly within Lampropeltis. Molecular Phylogenetics and Evolution 43: 674–684. Bryson RW, Smith BT, Nieto-Montes de Oca A, GarcíaVázquez UO, Riddle BR. 2014. Data from: The role of mitochondrial introgression in illuminating the evolutionary history of Nearctic treefrogs. Dryad Digital Repository. doi: 10.5061/dryad.3rc02. Carling MD, Brumfield RT. 2008. Haldane’s rule in an avian system: using cline theory and divergence population genetics to test for differential introgression of mitochondrial, autosomal, and sex-linked loci across the Passerina bunting hybrid zone. Evolution 62: 2600–2615. Carstens BC, Knowles LL. 2007. Estimating phylogeny from gene tree probabilities in Melanoplus grasshoppers despite incomplete lineage sorting. Systematic Biology 56: 400–411. Chan KMA, Levin SA. 2005. Leaky prezygotic isolation and porous genomes: rapid introgression of maternally inherited DNA. Evolution 59: 720–729. Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7: 214.

Duellman WE. 1970. The hylid frogs of Middle America. Monograph of the Museum of Natural History, The University of Kansas 1: 1–753. Duellman WE. 2001. The hylid frogs of middle America, 2nd edn. Lawrence, KS: Society for the Study of Amphibians and Reptiles. Durish ND, Halcomb KE, Kilpatrick CW, Bradley RD. 2004. Molecular systematics of the Peromyscus truei species group. Journal of Mammalogy 85: 1160–1169. Eckert AJ, Carstens BC. 2008. Does gene flow destroy phylogenetic signal? The performance of three methods for estimating species phylogenies in the presence of gene flow. Molecular Phylogenetics and Evolution 49: 832–842. Edwards SV, Bensch S. 2009. Looking forwards or looking backwards in avian phylogeography? A comment on Zink and Barrowclough 2008. Molecular Ecology 18: 2930–2933. Edwards SV, Liu L, Pearl DK. 2007. High resolution species trees without concatenation. Proceedings of the National Academy of Sciences, USA 104: 5936–5941. Eliosa-León HR. 2002a. Variación geográfica en Hyla eximia (Anura: Hylidae). MSc Thesis, Universidad Nacional Autónoma de México. Eliosa-León HR. 2002b. Variación geográfica en Hyla eximia (Anura: Hylidae). Boletin de la Sociedad Herpetologica Mexicana 10: 59–60. Fa JE, Morales LM. 1993. Patterns of mammalian diversity in Mexico. In: Ramamoorthy TP, Bye R, Lot A, Fa J, eds. Biological diversity in Mexico: origins and distribution. New York: Oxford University Press, 319–361. Fouquet A, Noonan BP, Rodrigues MT, Pech N, Gilles A, Gemmell NJ. 2012. Multiple Quaternary refugia in the eastern Guiana Shield revealed by comparative phylogeography of 12 frog species. Systematic Biology 61: 461– 489. Funk DJ, Omland KE. 2003. Species-level paraphyly: frequency, causes, and consequences, with insights from animal mitochondrial DNA. Annual Review of Ecology and Systematics 34: 397–423. Good JM, Hird S, Reid N, Demboski JR, Steppan SJ, Martin-Nims TR, Sullivan J. 2008. Ancient hybridization and mitochondrial capture between two species of chipmunks. Molecular Ecology 17: 1313–1327. Gross JB. 2012. The complex origin of Astyanax cavefish. BMC Evolutionary Biology 12: 105. Heled J, Drummond A. 2010. Bayesian inference of species trees from multilocus data. Molecular Biology and Evolution 27: 570–580. Holloway AK, Cannatella DC, Gerhardt HC, Hillis DM. 2006. Polyploids with different origins and ancestors form a single sexual polyploidy species. American Naturalist 167: E88–E101. Hua X, Fu C, Li J, Nieto-Montes de Oca A, Wiens JJ. 2009. A revised phylogeny of Holarctic treefrogs (genus Hyla) based on nuclear and mitochondrial DNA sequences. Herpetologica 65: 246–259. Jezkova T. 2010. Reconstructing species responses to past climatic changes using niche modeling and genetic data. PhD Thesis, University of Nevada Las Vegas.

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

INTROGRESSION IN NEARCTIC TREEFROGS Joly S. 2012. JML: testing hybridization from species trees. Molecular Ecology Resources 12: 179–184. Joly S, McLenachan PA, Lockhart PJ. 2009. A statistical approach for distinguishing hybridization and incomplete lineage sorting. American Naturalist 174: e54–e70. Jones KB. 1995. Phylogeography of the desert horned lizard (Phrynosoma platyrhinos) and the short-horned lizard (Phrynosoma douglassi): patterns of divergence and diversity. PhD Thesis, University of Nevada Las Vegas. Katoh K, Misawa K, Kuma K, Miyata T. 2002. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Research 30: 3059– 3066. Katoh K, Toh H. 2008. Improved accuracy of multiple ncRNA alignment by incorporating structural information into a MAFFT-based framework. BMC Bioinformatics 9: 212. Klymus KE, Gerhardt HC. 2012. AFLPs resolve intraspecific relationship and infer genetic structure among lineages of the canyon treefrog, Hyla arenicolor. Molecular Phylogenetics and Evolution 65: 654–667. Klymus KE, Humfeld SC, Gerhardt HC. 2012. Geographic variation in male advertisement calls and female preference of the wide-ranging canyon treefrog, Hyla arenicolor. Biological Journal of the Linnean Society 107: 219–232. Klymus KE, Humfeld SC, Marshall VT, Cannatella D, Gerhardt HC. 2010. Molecular patterns of differentiation in canyon treefrogs (Hyla arenicolor): evidence for introgressive hybridization with the Arizona treefrog (H. wrightorum) and correlations with advertisement call differences. Journal of Evolutionary Biology 23: 1425–1435. Leaché AD. 2010. Species trees for spiny lizards (genus Sceloporus): identifying points of concordance and conflict between nuclear and mitochondrial data. Molecular Phylogenetics and Evolution 54: 162–171. Leaché AD, McGuire JA. 2006. Phylogenetic relationships of horned lizards (Phrynosoma) based on nuclear and mitochondrial data: evidence for a misleading mitochondrial gene tree. Molecular Phylogenetics and Evolution 39: 628– 644. Linnen CR, Farrell BD. 2007. Mitonuclear discordance is caused by rampant mitochondrial introgression in Neodiprion (Hymenoptera: Diprionidae) sawflies. Evolution 61: 1417– 1438. Liu L, Pearl DK. 2007. Species trees from gene trees: reconstructing Bayesian posterior distributions of a species phylogeny using estimated gene tree distributions. Systematic Biology 56: 504–514. Luxbacher AM, Knouft JH. 2009. Assessing concurrent patterns of environmental niche and morphological evolution among species of horned lizards (Phrynosoma). Journal of Evolutionary Biology 22: 1669–1678. Marková S, Dufresne F, Manca M, Kotlík P. 2013. Mitochondrial capture misleads about ecological speciation in the Daphnia pulex complex. PLoS ONE 8: e69497. Martin DP, Lemey P, Lott M, Moulton V, Posada D, Lefeuvre P. 2010. RDP3: a flexible and fast computer program for analyzing recombination. Bioinformatics 26: 2462– 2463.

115

Monaghan MT, Wild R, Elliot M, Fujisawa T, Balke M, Inward DJG, Lees DC, Ranaivosolo R, Eggleton P, Barraclough TG, Vogler AP. 2009. Accelerated species inventory on Madagascar using coalescent-based models of species delineation. Systematic Biology 58: 298–311. Moore WS. 1995. Inferring phylogenies from mtDNA variation: mitochondrial gene trees versus nuclear gene trees. Evolution 49: 718–726. Mulcahy DG, Spaulding A, Mendelson JR, Brodie ED Jr. 2006. Phylogeography of the Flat-tailed Horned lizard (Phrynosoma mcallii) and systematics of the P. mcallii– platyrhinos mtDNA complex. Molecular Ecology 15: 1807– 1826. Nevado B, Fazalova V, Backeljau T, Hanssens M, Verheyen E. 2011. Repeated unidirectional introgression of nuclear and mitochondrial DNA between four congeneric Tanganyikan cichlids. Molecular Biology and Evolution 28: 2253–2267. Nylander JAA. 2004. MRMODELTEST v2. Program distributed by the author. Uppsala: Evolutionary Biology Centre, Uppsala University. Paniagua LL, Morrone JJ. 2009. Do the Oaxacan Highlands represent a natural biotic unit? A cladistic biogeographical test based on vertebrate taxa. Journal of Biogeography 36: 1939–1944. Parra-Olea G, García-París M, Papenfuss TJ, Wake DB. 2005. Systematics of the Pseudoeurycea bellii (Caudata: Plethodontidae) species complex. Herpetologica 61: 145–158. Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, Kamoun S, Sumlin WD, Vogler AP. 2006. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Systematic Biology 55: 595– 609. Puebla-Olivares F, Bonaccorso E, Espinosa de los Monteros A, Omland KE, Llorente-Bousquets JE, Peterson AT, Navarro-Sigüenza AG. 2008. Speciation in the emerald toucanet (Aulacorhynchus prasinus) complex. The Auk 125: 39–50. Pyron RA, Wiens JJ. 2011. A large-scale phylogeny of Amphibia including over 2,800 species, and a revised classification of extant frogs, salamanders, and caecilians. Molecular Phylogenetics and Evolution 61: 543–583. R Development Core Team. 2011. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing. Available at: http://www.R-project.org Rabosky DL, Talaba AL, Donnellan SC, Lovette IJ. 2009. Molecular evidence for hybridization between two Australian desert skinks, Ctenotus leonhardii and Ctenotus quattuordecimlineatus (Scincidae : Squamata). Molecular Phylogenetics and Evolution 53: 368–377. Rambaut A, Drummond AJ. 2007. Tracer v1.5. Available at: http://beast.bio.ed.ac.uk/Tracer Reid NM, Demboski JR, Sullivan J. 2012. Phylogeny estimation of the radiation of western North American chipmunks (Tamias) in the face of introgression using reproductive protein genes. Systematic Biology 61: 44–62. Ronquist F, Huelsenbeck JP. 2003. MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574.

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017

116

R. W. BRYSON JR ET AL.

Ruane S, Bryson RW, Pyron RA, Burbrink FT. 2014. Coalescent species delimitation in milksnakes (genus Lampropeltis) and impacts on phylogenetic comparative analyses. Systematic Biology 63: 231–250. Rubinoff D, Holland BS. 2005. Between two extremes: mitochondrial DNA is neither the panacea nor the nemesis of phylogenetic and taxonomic inference. Systematic Biology 54: 952–961. Ruiz-Sanchez E, Specht CD. 2013. Influence of the geological history of the Trans-Mexican Volcanic Belt on the diversification of Nolina parviflora (Asparagaceae: Nolinoideae). Journal of Biogeography 40: 1336–1347. Satler JD, Carstens BC, Hedin M. 2013. Multilocus species delimitation in a complex of morphologically conserved trapdoor spiders (Mygalomorphae, Antrodiaetidae, Aliatypus). Systematic Biology 62: 805–823. doi:10.1093/sysbio/syt041. Shields GF, Adams D, Garner G, Labelle M, Pietsch J, Ramsay M, Schwartz C, Titus K, Williamson S. 2000. Phylogeography of mitochondrial DNA variation in brown bears and polar bears. Molecular Phylogenetics and Evolution 15: 319–326. Smith SA, Nieto Montes de Oca A, Reeder TW, Wiens JJ. 2007. A phylogenetic perspective on elevational species richness patterns in Middle American treefrogs: why so few species in lowland tropical rainforests? Evolution 61: 1188–1207. Stephens M, Donnelly P. 2003. A comparison of Bayesian methods for haplotype reconstruction from population genotype data. American Journal of Human Genetics 73: 1162– 1169. Takahata N, Slatkin M. 1984. Mitochondrial DNA gene flow. Proceedings of the National Academy of Sciences, USA 81: 1764–1767. Tamura K, Peterson D, Peterson N, Stecher G, Nei M,

Kumar S. 2011. MEGA5: Molecular Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolution 28: 2731–2739. Tang Q-Y, Liu S-Q, Yu D, Liu H-Z, Danley PD. 2012. Mitochondrial capture and incomplete lineage sorting in the diversification of balitorine loaches (Cypriniformes, Balitoridae) revealed by mitochondrial and nuclear genes. Zoologica Scripta 41: 233–247. Tanner WW, Dixon JR, Harris HS. 1972. A new subspecies of Crotalus lepidus from western Mexico. Great Basin Naturalist 32: 16–24. Toews DPL, Brelsford A. 2012. The biogeography of mitochondrial and nuclear discordance in animals. Molecular Ecology 21: 3907–3930. Webb RG. 1984. Herpetogeography in the Mazatlán-Durango region of the Sierra Madre Occidental, Mexico. In: Seigel RA, Hunt LE, Knight JL, Malaret L, Zuschlag NL, eds. Vertebrate ecology and systematics – a tribute to Henry S. Fitch. Lawrence: Museum of Natural History, University of Kansas, 217–241. Willis SC, Farias IP, Ortí G. 2014. Testing mitochondrial capture and deep coalescence in Amazonian cichlid fishes (Cichlidae: Cichla). Evolution 68: 256–268. doi:10.1111/ evo.12230. Willis SC, Nunes MS, Montana CG, Farias IP, Lovejoy NR. 2007. Systematics, biogeography, and evolution of the Neotropical peacock basses Cichla (Perciformes: Cichlidae). Molecular Phylogenetics and Evolution 44: 291–307. Zakharov EV, Lobo NF, Nowak C, Hellmann JJ. 2009. Introgression as a likely cause of mtDNA paraphyly in two allopatric skippers (Lepidoptera: Hesperiidae). Heredity 102: 590–599.

SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s web-site: Figure S1. Ultrametric trees for the proopiomelanocortin A (POMC), tyrosinase (Tyr1), and β-crystallin (cryB) genes coestimated with the nuclear DNA species tree in *BEAST. Nodes that received ≥ 0.95 Bayesian posterior probability support are depicted with black dots. Table S1. Voucher data for treefrogs used in this study. Table S2. Occurrences of potential introgression estimated by JML from the mitochondrial DNA species-tree output with P-values < 0.05.

SHARED DATA Data deposited in the Dryad Digital Repository (Bryson et al., 2014).

© 2014 The Linnean Society of London, Zoological Journal of the Linnean Society, 2014, 172, 103–116

Downloaded from https://academic.oup.com/zoolinnean/article-abstract/172/1/103/3797036 by guest on 30 November 2017