## Abstract

Mutation frequencies can be modeled as a Poisson random field (PRF) to estimate speciation times and the degree of selection on newly arisen mutations. This approach provides a quantitative theory for comparing intraspecific polymorphism with interspecific divergence in the presence of selection and can be used to estimate population genetic parameters. Although the original PRF model has been extended to more general biological settings to make statistical inference about selection and divergence among model organisms, it has not been incorporated into phylogeographic studies that focus on estimating population genetic parameters for nonmodel organisms. Here, we modified a recently developed time-dependent PRF model to independently estimate genetic parameters from a nuclear and mitochondrial DNA data set of 22 sister pairs of birds that have diverged across a biogeographic barrier. We found that species that inhabit humid habitats had more recent divergence times and larger effective population sizes than those that inhabit drier habitats, and divergence time estimated from the PRF model were similar to estimates from a coalescent species-tree approach. Selection coefficients were higher in sister pairs that inhabited drier habitats than in those in humid habitats, but overall the mitochondrial DNA was under weak selection. Our study indicates that PRF models are useful for estimating various population genetic parameters and serve as a framework for incorporating estimates of selection into comparative phylogeographic studies.

A range of analytical approaches have been developed to estimate divergence times and effective population sizes from DNA data (Takahata *et al.* 1995; Yang 1997, 2002; Hickerson *et al.* 2003; Rannala and Yang 2003; Hey and Nielsen 2004; Heled and Drummond 2010). These types of models have been applied to an array of empirical systems and used to reconstruct population history (*e.g.*, Jennings and Edwards 2005; Won and Hey 2005; Lee and Edwards 2008; Hurt *et al.* 2009; Hailer *et al.* 2012), but the extension of these analyses to incorporate selection has been limited. A potential means of including selection estimates into studies examining population or phylogeographic history is the implementation of approaches that compare polymorphic and fixed nucleotide patterns within and between species (Hudson *et al.* 1987). Using this method, McDonald and Kreitman (1991) proposed a statistical test for homogeneity of synonymous and nonsynonymous DNA mutations among closely related species and demonstrated that an excess of nonsynonymous differences among fruit fly species was the result of the fixation of beneficial mutations.

The McDonald–Kreitman test was extended and incorporated into a quantitative model used to compare intraspecific polymorphism with interspecific divergence by modeling mutation frequencies as a Poisson random field (PRF) (Sawyer and Hartl 1992). The PRF models the number of synonymous and nonsynonymous fixed differences and polymorphic sites as independent Poisson random variables. The expected values of the Poisson counts are derived by approximating discrete time, discrete state Markov chains, using diffusion processes. Maximum likelihood and Bayesian inference have been used to estimate synonymous and nonsynonymous mutation rates, the selection coefficient of nonsynonymous mutations, and interspecific divergence times by comparing theoretical expected values with observed counts in PRF models (Sawyer and Hartl 1992; Bustamante *et al.* 2002).

The first developed PRF model assumed equal fitness for amino acid changes on the same locus, constant population sizes, random mating, no migration between species, genic selection, independence among nucleotide sites, and species at mutation–selection–drift equilibrium after divergence (Sawyer and Hartl 1992). These assumptions are likely violated in natural populations and several studies have attempted to make the PRF model more biologically realistic by relaxing model assumptions. Extensions of the original model include incorporating migration in an island model of subdivision (Wakeley 2008), a random effects model (Sawyer *et al.* 2003, 2007), relaxation of the genic selection assumption (Williamson *et al.* 2004), and a time-inhomogeneous model (Williamson *et al.* 2005; Boyko *et al.* 2008). The assumption of mutation–selection–drift equilibrium has been integrated into a time-dependent PRF model (Amei and Sawyer 2010) and shown to yield robust divergence time and selection estimates (Amei and Sawyer 2012). The greater biological realism of PRF models may make them attractive for analyses of comparative phylogeographic studies that focus on estimating population genetic parameters in nonmodel organisms.

Phylogeographic studies often use multilocus data sets composed of mitochondrial DNA (mtDNA) and nuclear DNA to estimate population parameters (Brito and Edwards 2009). The fast mutation rate and maternal inheritance of mtDNA make it a presumably more informative marker than nuclear DNA for identifying population structure in species (Zink and Barrowclough 2008). However, mtDNA and nuclear DNA can exhibit discordant phylogeographic patterns due to sex-biased dispersal (Peters *et al.* 2012), introgression (Bryson *et al.* 2010; Reid *et al.* 2012), and/or incomplete lineage sorting (Ballard and Whitlock 2004). Further, the utility of mtDNA as a phylogenetic marker has been questioned because of evidence of selective sweeps causing mtDNA genetic diversity to be incongruent with species abundance estimates (Bazin *et al.* 2006; Nabholz *et al.* 2008, 2009). Since phylogeographic studies often do not measure selection, the usage of PRF models to estimate divergence times, effective population sizes, and selection coefficients provides a means of incorporating the informativeness of mtDNA and at the same time assessing the strength of selection on mtDNA markers.

Comparative phylogeography allows for the testing of hypotheses that have shaped the diversification history of multiple codistributed taxa. There are two general analytical approaches for estimating and inferring comparative phylogeographic patterns (Hickerson *et al.* 2010). The first one is a model-based approach that compares the likelihoods of alternative diversification scenarios and attempts to find the model that best explains a phylogeographic pattern. The other approach relies on estimating phylogeographic parameters from phylogenetic trees and coalescent modeling and then assesses the level of congruence in parameters among codistributed species. Both of these methodological frameworks have been used to evaluate the historical assembly of bird faunas from a range of systems, including the Neotropics (Burney and Brumfield 2009; Smith *et al.* 2013), the Afrotropics (Voelker *et al.* 2013), Australia (Huang *et al.* 2011; Dolman and Joseph 2012), montane forests in Mexico (Barber and Klicka 2010), the Malay Archipelago (Lim *et al.* 2011), North American deserts (Zink *et al.* 2001), and the Caribbean (Sly *et al.* 2011). An overall goal of both analytical approaches is to evaluate whether patterns inferred from neutral genetic variation are consistent with the predictions of geological changes or climatic fluctuations. While comparative phylogeographic studies serve as a powerful approach for detecting shared evolutionary histories among species, the overall role of selection in the historical assembly of these bird faunas remains unclear.

In this study, we used a recently developed PRF inference method (Amei and Sawyer 2012) and applied it to a published comparative phylogeographic data set. We modified this PRF model by applying a two-step approach that uses both maximum likelihood and a Bayesian framework to independently estimate genetic parameters in 22 sister pairs of birds that have diverged across the Isthmus of Panama, a biogeographic barrier dominated by lowland rainforest that connects North and South America. Based on the mammalian fossil record, it has been proposed that during periods of climate change this rainforest corridor contracted, allowing drier habitat taxa to more readily disperse between continents (Webb and Rancy 1996). We tested this hypothesis by using 22 sister pairs of birds that have been delimited into climatic assemblages of either humid rainforest habitat or drier and more open habitat (Smith *et al.* 2012). If episodes of contracting and expanding rainforest were responsible for speciation, then humid and dry habitat taxa should exhibit temporal differences in divergence times. Similar to the pattern in the mammalian fossil record, we would expect more recent divergence times for humid habitat taxa. Previous work using a multispecies coalescent approach inferred that these humid and dry habitat assemblages of taxa may have diverged across the Isthmus of Panama during alternating episodes of contracting and expanding rainforest (Smith *et al.* 2012). However, taxa likely experience differing selective pressures and demographic histories in dry and humid habitats, and the magnitude of these differences is unclear. Here, we use a PRF model to estimate divergence times, effective population sizes, and selection coefficients to test whether humid and dry habitat assemblages have differing evolutionary histories. Understanding the timing of speciation between assemblages of taxa that inhabit humid and dry habitats may help explain the processes that generated high diversity in tropical regions.

## Methods

### Data

We used a published multilocus data set of sister pairs of birds that exhibit genetic breaks across the Isthmus of Panama (Smith *et al.* 2012). The original data set was composed of 29 sister pairs, but some of the pairs exhibited considerable genetic structure that was generated after divergence across the Isthmus of Panama. This genetic structure violates the random mating assumption of the PRF model, so we used only 22 of the 29 pairs that do not have deep genetic structure in populations on each side of the Isthmus of Panama. Based on the results of a climatic assemblage assignment test, 13 pairs were assigned to the humid habitat group and 9 were assigned to a drier habitat group, but there was greater climatic variance among and within pairs assigned to the dry habitat group (Smith *et al.* 2012). Taxon names and corresponding climatic assemblages are listed in Supporting Information, Table S1, and further details on the taxa and data are available in Smith *et al.* (2012). For each species pair, we used sequence alignments of five non-protein-coding markers from the nuclear genome, aconitase 1(ACO1), β-fibrinogen intron 5 (βFib5), myoglobin intron 2, eukaryotic translation elongation factor 2 (EEF2), and ornithine decarboxylase (ODC), and one protein-coding gene from mitochondrial genome, NADH dehydrogenase subunit 2 (ND2), to generate McDonald–Kreitman tables for all loci in the program DnaSP (Librado and Rozas 2009). The number of alleles across the six loci ranged from 1 to 18.

### Model

We followed a two-step procedure to estimate species divergence time and infer selection on nonsynonymous mutations. First, assuming that selective effects of segregating mutations of noncoding genes are neutral, we used Markov chain Monte Carlo simulations to estimate divergence time for each of the 22 species pairs by applying a time-dependent PRF model to DNA alignments of the five noncoding genes. Second, given an estimated divergence time of a species pair, we inferred the strength and direction of selection by comparing theoretical values with observed counts of the corresponding McDonald–Kreitman (M-K) table generated from the single coding gene.

In the time-dependent PRF model, the traditional 2 × 2 M-K table was generalized to a 2 × 3 contingency table that included two different types of polymorphisms (Amei and Sawyer 2010, 2012). Adopting the conventional notations, we used to represent the number of synonymous sites that are fixed between a pair of species in a sample of aligned DNA sequences from the two species, to represent the number of synonymous sites that are polymorphic in only one species sample, and to be the number of synonymous sites that are polymorphic in both samples. For each species pair, we assumed that the two species have the same haploid effective population size they diverged generations ago, and the synonymous mutation rate of sites in a noncoding gene in the two lineages is per gene per generations. Now, suppose that we have *m* aligned DNA sequences from one species and *n* from the other species. Under the assumption that nucleotide sites evolve independently, the counts and are distributed as independent Poisson random variables with means expressed by functions of the genetic parameters and (explicit expressions are in File S1). We applied a Bayesian approach by treating the two model parameters and as random variables with prior distributions uniform and gamma, respectively. Given the observed counts and at each locus *i*, the full-likelihood function based on the sampling formulas (Equations S1–S3 in File S1) and the prior distributions is given by (1)where *M* is the number of loci ( in our sample), and are gamma and uniform densities, andNumerical solutions for integrations are obtained by using the Crank–Nicolson method for integrals involving the transition density (defined in File S1) and Gauss–Legendre quadrature for all other types (Press *et al.* 2007). All sites in the five noncoding genes are treated as silent sites to calculate observed values of and for Setting the hyperparameters and to a small value () and *T* to a fixed large value, posterior distributions of the model parameters and are obtained by MCMC simulations with convergence ensured by trace plots and Gelman–Rubin statistics (Gelman 1996). For each species pair, Markov chain Monte Carlo simulations are implemented 120,000 times, the first 20,000 burn-in iterations are disregarded, and model parameters are estimated from 10,000 samples taken every 10 iterating steps. Based on the derivation of the PRF model, is an aggregated synonymous mutation rate at locus *i*, per generations (Amei and Sawyer 2010). We divided this quantity by the total number of synonymous sites at locus *i* to get the rate of substitution per synonymous site per generations and denote it by Based on our estimates, the synonymous mutation rates of the five loci are similar; thus we used the mean of the five loci, as a final estimate of the synonymous mutation rate per synonymous site per generations. A substitution rate of per site per million years (MY), based on a published intron rate for birds (Axelsson *et al.* 2004; Ellegren 2007), was used to express the generations in a unit of MY and the haploid effective population size ( for diploid population size) was obtained by assuming a generation time of 1 year. Consequently, the species divergence time in units of MY is given by multiplying the model-estimated by

Next, we compared the ND2 M-K table with theoretical expected values to infer selection. Let and be the amino acid nonsynonymous counterparts of the previously defined and Under the same model assumptions, the variables are independent Poisson random variables and their expected values are given explicitly in Amei and Sawyer (2012). The mathematical expressions of these Poisson means are similar to those given in Equations S1–S3 in File S1 but contain the selection coefficient γ. The magnitude and sign of γ describe the strength and direction of the selection intensity of amino acid nonsynonymous mutations. We assumed that polymorphic or fixed mutations have the same selection coefficients within a gene and estimated γ by setting the ratio of the expected number of nonsynonymous polymorphic sites to the expected number of nonsynonymous fixed differences to equal the ratio of the numbers of observed nonsynonymous polymorphic sites and fixed differences. That is, (2)Since the mutation parameter (nonsynonymous mutation rate) cancels in the ratio of the expected counts, Equation 2 becomes an equation of γ and Using the estimated from step 1, the numerical solution of the selection coefficient γ is obtained by the Newton–Raphson method and it can be shown that the obtained value is the same as its maximum-likelihood estimate (Sawyer and Hartl 1992). In addition to the point estimate of γ, we applied a 95% credible interval (CI) of estimated from step 1, to Equation 2 to obtain an interval estimate of γ. Specifically, for each of the 22 species pairs, the lower and upper bounds of the interval estimate of γ were calculated by substituting the upper and lower bounds of the 95% CI of the corresponding

## Results

### Effective population size

Haploid effective population size ( for diploid population size) for the 22 bird species pairs ranged from 440,000 to 3,640,000 (Table S1). Median (middle line of the box) values are 900,000 for dry habitat pairs and 1,550,000 for humid habitat pairs (Figure 1). Overall, estimates are smaller in dry habitat pairs () than in humid habitat pairs (). One exception was *Troglodytes aedon*, which has a much larger than the other dry habitat pairs.

### Divergence times

We report medians and 95% credible intervals for divergence times estimated from the time-dependent PRF model, using the five noncoding nuclear loci. Median estimates of the posterior distributions ranged from 0.79 to 4.83 MY ago for dry habitat pairs and 0.36 to 2.81 MY ago for humid habitat pairs (Table S1). The side-by-side boxplot in Figure 1 shows that the interquartile ranges (height of the box) of the are 3.12 MY ago for the 9 dry habitat bird species and 1.14 MY for the 13 humid habitat species. Although the interquartile ranges of the estimates differ by almost 2 MY between the two habitat groups, the median values of the two groups are quite close. Median values were 1.88 MY ago for dry habitat pairs and 1.13 MY ago for humid habitat pairs. The median value for all pairs from the PRF model was 1.23 MY ago (Table S1) in comparison to 1.92 MY ago from the *BEAST divergence times generated in Smith *et al.* (2012). For comparison, the *BEAST estimates in Smith *et al.* (2012) were generated using the following settings: relaxed uncorrelated molecular clocks based on published substitution rates (Axelsson *et al.* 2004), a Yule process on the species-tree prior, and best-fit models of sequence evolution. The PRF model estimates were typically younger than the *BEAST time estimates, but were not significantly different (two tailed Mann–Whitney test, *u* = 0.939, *P* = 0.358). The 95% credible intervals from the PRF model and *BEAST overlapped for all species except for *Xenops minutus* (Table S1; Figure 2). Median divergence times estimated from the PRF model between the dry habitat group and the humid habitat group were not statistically significant (two sided Mann–Whitney test, *u* = 82, *P* = 0.125).

### Selection coefficients

For selection coefficients we report values from the protein-coding mitochondrial gene ND2. The selection parameter γ of nonsynonymous mutation was scaled in terms of the haploid effective population size that is, where *s* is the same selection coefficient per generation. For each of the 22 species pairs, we obtained a point and an interval estimate of the selection coefficient *s* by applying the ND2 M-K table as well as the estimated values of and from step 1. Selection coefficients varied across the 22 sister pairs (Table 1) and were different between the two climatic assemblages. Among the point estimates of the *s* of the 22 species pairs, 19 have their changing from to and the other 3 have negative selection coefficients of and respectively. For the majority of the species pairs, selection coefficients were positive for nonsynonymous changes, but the strength of selection was weak. There are 11 species pairs whose selection coefficients were < 14 pairs’ selection coefficients were < and 17 had Overall, selection coefficients were relatively larger both in magnitude and in variation for dry habitat birds as opposed to small and consistent selection for humid habitat birds. A comparison of selection coefficients between habitat groups (Figure 1) indicated that the estimated *s* values in dry habitat bird pairs ranged from to with an interquartile range of In contrast, the *s* values for the 13 humid habitat pairs ranged from to with a much smaller interquartile range of Two species pairs had larger *s* values (*s* = *Campylorhynchus chiapensis* and *C. griseus*, a dry habitat pair; and *s* = *Cyanocompsa cyanoides*, a humid habitat pair) than the other species pairs. To analyze the effect of the biased divergence time estimate on the selection coefficient, we generated an interval estimate for *s* (Table 1; Figure 3). Except for the abovementioned two species pairs, mean lengths of the interval estimates of the selection coefficient are small, ranging from 1.59 for the dry habitat group to 0.93 for the humid habitat group. There are 6 species pairs whose interval estimates contain zero and the rest are higher than the neutral line (Figure 3).

## Discussion

We developed a two-step approach to evaluate divergence times and selection patterns in humid and dry assemblages of birds that are distributed across the Isthmus of Panama. Our study detected temporal differences in divergence times, effective population sizes, and selection coefficients between the taxa that inhabit humid and drier habitats. There was greater variance in all three parameter estimates among the dry habitat pairs that had older divergence times, smaller effective populations, and larger selection coefficients than the humid habitat pairs. The PRF model as a comparative phylogeographic tool provides an analytical framework similar to tree-based methods that yield independent estimates of population genetic parameters that can be used to assess congruence among codistributed species. Overall, we found that the PRF model is useful for comparative phylogeographic studies interested in both the timing of divergence and the role of selection in shaping the evolutionary history of communities.

### Estimating divergence time with a PRF model

Our estimates of species divergence times indicated that dry habitat birds began diverging as early as 4.8 MY ago and continued to 0.8 MY ago. In contrast, speciation of humid habitat birds started much later and went through less variation with an estimated of 2.81 MY ago to a latest time of 0.36 MY ago. The older divergence times in dry habitat birds correspond to the pattern observed in the mammalian fossil record that drier habitat taxa were the first to disperse across the Isthmus of Panama (Webb and Rancy 1996), but the median divergence times among groups were not significantly different. This finding is consistent with the conclusions of several comparative phylogeographic studies that recovered multiple episodes of divergence across barriers within a single assemblage of birds (Barber and Klicka 2010; Lim *et al.* 2011; Sly *et al.* 2011). Also, the assignment of sister pairs into discrete climatic groups may have limited our ability to detect significant differences among groups. Climatic preferences across the 22 sister pairs represent a gradient from humid- to dry-adapted species, and binning the pairs into discrete groups may mask more nuanced diversification patterns.

The divergence times estimated using the PRF model detected qualitatively similar speciation patterns to the results based on a species-tree approach (Smith *et al.* 2012), although the PRF model-based estimates were more recent (Figure 2). The observed differences between the models are likely attributable to *BEAST having more parameters because it incorporates more complex models of sequence evolution and a relaxed molecular clock (Drummond and Rambaut 2007; Heled and Drummond 2010). The incorporation of variable mutation rates (Yang 1997, 2002; Burgess and Yang 2008) or a relaxed molecular clock (Drummond *et al.* 2006) in our model would more accurately capture the mutational process, but the similar results between *BEAST and the PRF model along with simulation studies (Amei and Sawyer 2012) suggest that our model provides accurate parameter estimates. Despite the current limitations of the PRF model, our approach yields reconstructions of population history comparable to patterns estimated from other models and recovered similar complex patterns to those that are often observed in empirical systems (*e.g.*, Lim *et al.* 2011).

### Comparative population sizes

The estimation of effective population size in divergence time models accounts for the impact of ancestral population size on divergence time estimation (Edwards and Beerli 2000; Takahata *et al.* 1995) and provides the size of extant populations. We found that humid habitat pairs had larger median than dry habitat pairs (Figure 1). This finding is consistent with the differences in the distributions of humid and dry tropical biomes in Central and South America. Dry biomes have smaller geographic distributions and are highly fragmented in comparison to humid biomes that are more continuously distributed from Mexico to Bolivia (Prado and Gibbs 1993). The larger geographic area that humid biomes encompass presumably allows birds occurring in these habitats to have larger population sizes. Our inferred estimates of are consistent with the observation that genetic diversity is correlated with range size (Frankham 1996).

### Selection on mtDNA

We inferred that the ND2 gene showed evidence of positive selection in most of the sister pairs (Figure 3). The magnitude of the estimated selection coefficients was < for most pairs and only three pairs had negative selection coefficients with absolute values < There were two pairs in our data set, each from one climatic group, estimated to have much larger positive selection coefficients than the other pairs. Based on our results it is not directly apparent why these sister pairs would be under greater selection than the other species. One possible explanation is that nonsynonymous mutations are being fixed at an accelerated rate in these taxa because they have much smaller (Smith and Klicka 2013). Additional analyses would be necessary to investigate the role of because the PRF model assumes equal for both lineages. However, in the case of *C. chiapensis* (from the *C. chiapensis*/*C. griseus* pair), the species likely has a small because it has a small geographic range along the coast of Mexico. The observed differences in selection coefficients between dry and humid habitat birds require further inquiry. For example, species that inhabit drier habitats may be under greater selection and larger multilocus data sets may detect more apparent evidence of selection observed in mtDNA, such as adaptive selective sweeps (Ballard and Whitlock 2004; Bazin *et al.* 2006).

Selection coefficients were sensitive to the uncertainty surrounding divergence time. To evaluate the influence divergence time had on the selection coefficient, we estimated *s* using different point estimates of the lower and upper bounds of the 95% CI of When the time estimate was larger than the median estimate, the selection coefficient decreased, and, vice versa, the selection coefficient increased when divergence time was smaller than the median time estimate. There were six sister pairs whose selection coefficient interval estimates contained zero and the rest were above the neutral line. Thus, for the majority of sister pairs the direction of selection remained the same when we used the interval estimate (Table 1; Figure 3). But, for the pairs that contained zero in their interval estimate we could not distinguish between selection and neutral evolution (Table 1; Figure 3). The effective population size and the scaled selection coefficient γ are estimated independently and do not bias the γ-value. However, an over- or underestimated would decrease or increase the nonscaled selection coefficient *s*, due to the relation Estimated selection coefficients were influenced by the uncertainty surrounding other parameters, but the general patterns remained robust.

The strength of selection on mtDNA has important implications for using mtDNA in phylogeographic studies (Edwards and Bensch 2009; Galtier *et al.* 2009). Previous work has found that selection on mtDNA makes it a poor indicator of population size (Bazin *et al.* 2006; Nabholz *et al.* 2009) or species history (Bensch *et al.* 2006) or is uncorrelated with life history traits (Nabholz *et al.* 2008), but other work suggests that selection is not strong enough to bias phylogenetic inference (Zink 2005). Our results that mtDNA was under weak positive selection differ from previous work showing purifying selection reduced mtDNA genetic variability (Bazin *et al.* 2006; Stewart *et al.* 2008), but evidence of positive selection on mtDNA is not a unique finding. A comparison of mainland–island birds found that the mtDNA from large mainland populations showed evidence of positive selection (Wright *et al.* 2009). The variances in the magnitude and type of selection between our study and previous work is unclear, but the discrepancies may be attributed to the use of different mtDNA genes (ND2 *vs.* cytochrome *b*), locus length, or calculating selection at differing phylogenetic scales. The overall weak selection we recovered suggests that mtDNA might reflect population history more accurately than previous studies have suggested, but our analyses cannot rule out that even low levels of selection can reduce mtDNA diversity levels within populations, thus limiting the utility of the marker.

### Conclusions

Phylogeographic studies, particularly those on birds, have focused on estimating parameters from neutrally evolving loci to infer patterns of diversification, while neglecting the role of selection (Hickerson *et al.* 2010). We demonstrated that the PRF model presented here provides an analytical framework for the incorporation of estimates of selection along with standard population genetic parameters into a comparative phylogeographic study. Our approach provided similar divergence time patterns to those estimated using a multispecies coalescent framework, despite the PRF model assuming a strict molecular clock, a single mutation rate across loci and sites, and the same effective population size for ancestral and daughter populations. These limiting factors can be incorporated into more complex PRF models and are poised to become increasingly useful as genetic data sets become larger. The development of next-generation sequencing technologies now makes it feasible to generate thousands of loci for phylogeographic studies (*e.g.*, Smith *et al.* 2013) and PRF models can efficiently estimate parameters from large multilocus data sets (Amei and Sawyer 2012). The incorporation of selection into phylogeographic studies will yield a more comprehensive view on genetic patterns across species ranges and the processes responsible for shaping their evolution.

## Acknowledgments

We thank J. Klicka for his help generating the multilocus data set. We also thank D. Fletcher, R. Brumfield three anonymous reviewers, and N. Rosenberg for providing comments and feedback that greatly improved the manuscript.

## Footnotes

*Communicating editor: N. A. Rosenberg*

- Received September 20, 2013.
- Accepted October 11, 2013.

- Copyright © 2014 by the Genetics Society of America