Heredity (2005) 94, 217–228

& 2005 Nature Publishing Group All rights reserved 0018-067X/05 $30.00
www.nature.com/hdy

Local adaptation in the rock pocket mouse

(Chaetodipus intermedius): natural selection and
phylogenetic history of populations
HE Hoekstra, JG Krenz1 and MW Nachman
Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ 85721, USA

Elucidating the causes of population divergence is a central Mantel tests, we show that there is no correlation between
goal of evolutionary biology. Rock pocket mice, Chaeotdipus color variation and mtDNA phylogeny, suggesting that
intermedius, are an ideal system in which to study pelage coloration has evolved rapidly. At a finer geographical
intraspecific phenotypic divergence because of the extensive scale, high levels of gene flow between neighboring melanic
color variation observed within this species. Here, we and light populations suggest the selection acting on color
investigate whether phenotypic variation in color is correlated must be quite strong to maintain habitat-specific phenotypic
with local environmental conditions or with phylogenetic distributions. Finally, we raise the possibility that, in some
history. First, we quantified variation in pelage color (n ¼ 107 cases, migration between populations of pocket mice
mice) and habitat color (n ¼ 51 rocks) using a spectro- inhabiting different lava flows may be responsible for similar
photometer, and showed that there was a correlation melanic phenotypes in different populations. Together, the
between pelage color and habitat color across 14 sampled results suggest that color variation can evolve very rapidly
populations (R2 ¼ 0.43). Analyses of mtDNA sequences from over small geographic scales and that gene flow can both
these same individuals revealed strong population structure hinder and promote local adaptation.
in this species across its range, where most variation (63%) Heredity (2005) 94, 217–228. doi:10.1038/sj.hdy.6800600
was partitioned between five geographic regions. Using Published online 3 November 2004

Keywords: adaptation; Chaetodipus; color; gene flow; phenotypic variation; phylogeography

Introduction local conditions. In other cases, pleiotropic effects of

otherwise beneficial mutations may limit their spread.
A central goal of developmental and evolutionary Here, we focus on the extent to which phenotypic
biology is to explain the morphological diversity variation is correlated with local environmental condi-
observed across species. As differences among species tions versus phylogenetic history.
must initially occur as intraspecific polymorphism, The rock pocket mouse, Chaetodipus intermedius,
understanding the causes of intraspecific variation can provides an excellent system to study geographic
provide information about the origin of species-level variation in phenotype within a single species, and
differences. One ongoing debate concerns the relative allows us to explore this variation in light of the
roles of deterministic evolutionary processes and histo- underlying genetic structure of this species. In particular,
rical contingency in shaping the outcome of evolution we have been interested in the evolution of color
(Travisano et al, 1995). It has long been appreciated that differences in response to local environmental condi-
strong selection can lead to local adaptation, provided tions. C. intermedius lives exclusively in rocky habitat
that it is not swamped by gene flow (Haldane, 1948; across the southwestern deserts, and thus C. intermedius
Slatkin, 1985; Lenormand, 2002). On the other hand, habitat is largely discontinuous through most of its
populations may sometimes be constrained by their range. Historically, C. intermedius comprises 10 subspe-
evolutionary history. For example, young populations in cies (Benson, 1933; Dice and Blossom, 1937; Weckerly,
new environments may not have had time to adapt to 1983). Several subspecies have been described based on
dramatic color differences on small isolated lava flows
Correspondence: HE Hoekstra, Current address: Ecology, Behavior and
(Benson, 1933); these mice have extremely dark coats and
Evolution Section, Division of Biological Sciences, University of uniformly melanic hairs. Non-lava-dwelling populations
California at San Diego, 9500 Gilman Drive, MC0116, Muir Biology also show variation in coat color, which often closely
Bldg, Rm 4256, La Jolla, CA 92093-0116, USA. resembles the substrate color on which the mice live
E-mail: hoekstra@ucsd.edu (Benson, 1933; Dice and Blossom, 1937).
1
Current address: JG Krenz, United States Department of Agriculture, Previous work identified the genetic basis of melanism
Agricultural Research Service-Aquaculture Genetics, Oregon State
University-Hatfield Marine Science Center, 2030 SE Marine Science Dr,
in a single population of mice inhabiting the Pinacate
Newport, OR 97365, USA. lava flow in southern Arizona (Nachman et al, 2003).
Received 14 April 2004; accepted 3 September 2004; published Four amino-acid changes in the coding region of the
online 3 November 2004 melanocortin-1 receptor (Mc1r) are perfectly associated
 Local adaptation in pocket mice
HE Hoekstra et al

218
with a coat color polymorphism in this population. A Methods
focused study on the Pinacate region suggested that
strong selection is maintaining Mc1r allele and coat color Phenotypic variation
frequencies across short geographic distances in the face Sampling: A total of 107 rock pocket mice were collected
of high countervailing gene flow (Hoekstra et al, 2004). using Sherman live traps from 14 localities across the
Interestingly, mutations in the Mc1r gene are not species range in Arizona, New Mexico and northern
involved in additional melanic populations in New Mexico (Figure 1, Table 1). To quantify variation in local
Mexico (Hoekstra and Nachman, 2003). While habitat, rocks were collected in areas neighboring the
genes unlinked to those that underlie the phenotype trap-lines at each site. Localities were chosen to maximize
of interest cannot provide a direct test of whether phenotypic and environmental variation, rather than to
melanism has arisen independently in New Mexico, cover the species range. In addition, several localities
they can provide information about population structure, represent paired sampling sites, where substrate color
levels of gene flow and timing of colonization, which differed dramatically over short geographical distances.
may have implications for the evolution of melanic Liver, kidney and spleen samples were taken from each
populations. individual and frozen at 801C. Voucher specimens were
More generally, the extensive phenotypic variation prepared and deposited in the Zoological Collections of
observed in C. intermedius raises several interesting the Department of Ecology and Evolutionary Biology at
questions. First, are morphological differences reflected the University of Arizona.
in the genetic structure of the species? How much
gene flow occurs between populations which differ
in phenotype? Finally, what is the time scale over Phenotypic and environmental variation: The reflectance
which this phenotypic variation has evolved? To of both mouse coat color and corresponding rocks was
address these questions, we have quantified measured as a percentage of a white standard using a
phenotypic variation in coat color and local substrate USB2000 spectrophotometer (Ocean Optics) with a dual
color across 14 populations of C. intermedius. We used a deuterium/halogen light source. A standard reflection
spectrophotometer to quantify color variation among probe with a 200m receptor fiber was held at an uniform
populations and sequenced two mtDNA genes to distance from the surface at 901 to capture both diffuse
characterize the population structure and extent of gene and spectral reflectance. Measurements were taken at
flow between populations relative to phenotypic differ- 1518 points from 300 to 800 nm, and thus included the
entiation. UV spectrum.

Figure 1 Phenotypic variation across the range of C. intermedius in the southwestern US. Photographs represent the typical dorsal coloration
of individuals from each of 14 collecting locales indicated by circles. Filled circles represent lava flows and open circles are nonvolcanic rocky
regions. The border between Arizona and New Mexico roughly represents the interface of the Sonoran and Chihuahuan deserts, respectively.

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

219
Table 1 Collecting localities
Locale County N Latitude/longitude Subspecies Specimen numbers

Black Tank Lava BLK Coconino Co., AZ 3 35135.40 /111138.10 C. i. crinitus HEH 615-7
White Hills WHT Yavapai Co., AZ 10 34139.30 /111154.00 C. i. umbrosus HEH 618-29
Tinajas Altas TIN Yuma Co., AZ 9 32118.80 /114102.90 C. i. phasma HEH 635-9, 644-5, 648-9
Tule Mts. TUL Yuma Co., AZ 10 32110.50 /113140.20 C. i. mearnsi HEH 546-8, 550-3, 589, 592, 603
Pinacate Lava PIN Yuma Co., AZ 10 32105.90 /113128.30 C. i. pinacate HEH 534-42, MWN 1434
Pinacate Mexico MEX Sonora, MX 3 32106.50 /113122.50 C. i. pinacate HEH 512-4
Avra Valley AVR Pinal Co., AZ 10 32123.80 /111108.50 C. i. intermedius MWN 1358-67
Tumamoc Lava TUM Pima Co., AZ 3 32113.00 /111100.40 C. i. nigrimontis HEH 509-511
Portal POR Cochise Co., AZ 5 31154.90 /109103.70 C. i. intermedius MWN 1394-5, 1399-1401
Afton Hill AFT Dona Ana Co., NM 6 32109.60 /107110.40 C. i. intermedius HEH 515-20
Kenzin Lava KNZ Dona Ana Co., NM 10 32103.40 /106157.90 C. i. rupestris HEH 504-6, 522-4, 526, 557, 559, 687
Fra Cristobol Mts. FRA Sierra Co., NM 10 33119.90 /107107.10 C. i. intermedius HEH 571, 574-582
Armendaris Lava ARM Socorro Co., NM 9 33131.00 /106160.00 (un-named) HEH 521, 561-8
Carrizozo Lava CAR Lincoln Co., NM 9 33141.50 /105150.90 C. i. ater MWN 1337-45

In all, 107 nonmolting adults and 51 rocks were maximum-likelihood (ML) algorithms. NJ trees were
analyzed from all 14 collecting sites. For each animal, generated using the transition/transversion (ti/tv) rate
10 measurements were made from the dorsal surface of and gamma distribution shape parameter (g) estimated
the mouse and averaged to produce a general descrip- from Modeltest v3.06 (Posada and Crandall, 1998).
tion of the dorsal coat color. Similarly, 10 measurements Parsimony genealogies were generated with transitions
from the exposed rock surface were taken and averaged. and transversions weighted equally, and also with
Analysis followed that described in Hoekstra and Nach- transversions given a weight two, five and 10 times
man (2003). Here, we report total reflectance (relative to a more than that of transitions; the observed number of
pure white standard) to characterize overall dorsal transitions (n ¼ 250) in our data was 4.6 times greater
pelage and rock surface coloration. than the observed number of transversions (n ¼ 55).
Heuristic searches were performed with stepwise
Genetic variation addition for initial trees and the tree-bisection-
Genomic DNA was prepared using a tissue extraction kit reconnection method of branch swapping. The
(Qiagen DNeasy Tissue kit). MtDNA genes, COIII consensus parsimony tree was used as the starting tree
(783 bp) and ND3 (345 bp), were amplified from each in the estimation of ML trees. Hierarchical likelihood-
individual. In some populations, the intergenic tRNA- ratio tests implemented in Modeltest were used to
Gly was included in the analysis, so the total genetic determine the best-fit model of nucleotide substitution
region analyzed was 1176 bp in length. Primer sequences to estimate phylogenetic trees and genetic distances.
and PCR conditions are reported in Hoekstra et al (2004). Topologies were explored using a branch-and-bound
PCR products were cleaned using spin columns (Qiagen) method. Confidence in the branching structure was
and were sequenced on an ABI 3700. Sequences were assessed by performing 1000 bootstrap replicates.
assembled and aligned in Sequencher (GeneCodes) and
checked by eye. Outgroup sequences were obtained from Population structure and gene flow: To test for
the sister species, C. penicillatus (GenBank AY259036) and significant population structure among populations and
C. baileyi (AY259035). Sequences have been deposited in biogeographical regions, analyses of molecular variance
Genbank (accession numbers AY694010-AY694094, see (AMOVA, Excoffier et al, 1992) were performed in
Appendix A1). ARLEQUIN (Schneider et al, 2000). Pairwise FST
estimates were permuted 1000 times. A one-factor
Summary statistics: The sequence alignment was AMOVA was employed to assess the degree of
imported into DnaSP v3.0 (Rozas and Rozas, 1999) to population structure over all populations. Broad-scale
calculate intra- and interspecific genetic variability. In patterns of regional diversity were examined by pooling
each population, the number of segregating sites and the samples within geographical regions identified as having
number of unique haplotypes were counted. The average a recent common history based on phylogenetic analyses.
number of pairwise differences, p (Nei and Li, 1979), These pooled samples also correlated to known
and diversity based on the number of segregating sites, biogeographical regions. We clustered populations at
y (Waterson, 1975), were calculated. To check for three levels: (1) Sonoran versus Chihuahuan deserts
deviations from neutral expectations in the frequency (Sonoran ¼ BLK þ WHT þ TIN þ TUL þ PIN þ MEX þ
spectrum of polymorphisms, significance values were AVR þ TUM and Chichuahuan ¼ POR þ AFT þ KNZ þ
calculated for Tajima’s D statistic (Tajima, 1989). FRA þ ARM þ CAR), (2) five clusters based on fine-
scale biogeographic regions (northern Arizona ¼ BLK þ
Phylogenetic analysis: To estimate the phylogenetic WHT, southern Arizona ¼ TIN þ TUL þ PIN þ MEX,
relationships among haplotypes and among popu- central Arizona ¼ TUM þ AVR, southern New Mexico ¼
lations, gene genealogies were constructed. Gene trees POR þ AFT þ KNZ, central New Mexico ¼ FRA þ ARM þ
were generated in PAUP* v4.02 (Swofford, 1999) using CAR) and (3) northern Arizona versus all other
neighbor-joining (NJ), maximum parsimony (MP) and populations (northern Arizona ¼ BLK þ WHT).

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

220
We tested for isolation by distance between popula- gram MDIV (Nielsen and Wakeley, 2001). This Markov
tions of C. intermedius by estimating pairwise effective chain Monte Carlo (MCMC) method allows us to
migration rates Nm ¼ [(1/FST)1]/2 between all possible determine whether shared polymorphisms are the result
pairs of populations, where N is the female effective of recurrent gene flow, recent common ancestry or both.
population size, m is the female migration rate and FST is We used a finite-sites mutation model (HKY) following
a measure of population structure (Slatkin, 1993). We also Palsboll et al (2004).
calculated FST following the method suggested by
Rousset (1997), which is a modification of Slatkin’s FST
using an isolation by distance approach, and results were
consistent using both methods. When populations are in Results
migration–drift equilibrium and isolated by distance, the
effective migration rate should be negatively correlated Phenotypic variation
with interpopulation geographic distance. A total of 107 mice and 51 rocks from 14 populations
across the species range were included in the spectro-
photometric measurements (Figure 2). The percent
Correlation between phenotype, genotype and geography reflectance from the dorsal mouse coat ranged from
To test for a significant correlation between phenotypic, 4.02 to 11.96%. The darkest six populations represented
genetic and geographical distance among populations of mice inhabiting lava flows, indicated by an asterisk in
C. intermedius, we performed Mantel tests using ARLE- Figure 2.
QUIN. Mantel procedures can test for an association Reflectance from the substrate was significantly higher
between two matrices using randomization (Manly, than mouse coat reflectance, ranging from 9.18 to 41.8%.
1986). Specifically, the parameter is compared to a There is a clear difference between nonvolcanic rocks (eg
distribution obtained when the matrix is repeatedly granite, gneiss and limestone) and rocks of volcanic
randomized, and the null hypothesis of no association is origins (t ¼ 7.43; Po0.0001). The mean percent reflec-
rejected when the parameter exceeds a given significance tance for nonvolcanic rocks was 7.93% (SD ¼ 0.31) and
level. In this case, we can test whether patterns of color for volcanic rocks was 4.49% (SD ¼ 0.35).
variation are significantly associated with geography There was a significant positive correlation between
(implying local adaptation) and/or genetic distances (the mean rock color and mean mouse color among the 14
influence of shared evolutionary history). We also populations sampled (R2 ¼ 0.429, Po0.05; Figure 2).
conducted partial Mantel tests, which hold one matrix However, there are two notable exceptions. First, one
constant and test for an association between the remain- population from northern Arizona, Black Tank Lava
ing two matrices. As the appropriateness of partial (BLK), had relatively light-colored mice but the rock was
Mantel tests have recently been called into question, we very dark and volcanic in origin. Second, the rock
use the results from the partial Mantel tests only to substrate from the White Hills (WHT) in northern
support results from the Mantel tests (Raufaste and Arizona is extremely light-colored (the lightest reported
Rousset, 2001; Castellano and Balletto, 2002; Rousset, here); however, reflectance from the corresponding mice
2002). was unremarkable.

Population divergence and local adaptation

Genetic variation
We employed a molecular clock to generate estimates of The mtDNA loci surveyed in this study exhibit high
divergence times between clades. This clock is based on
levels of polymorphism in terms of number of haplo-
silent site divergence of mtDNA genes in rodents of types, segregating sites and nucleotide diversity as
approximately 2% per million years (Wilson et al, 1985).
measured by p and y (Table 2). Across the species, there
We compared average pairwise genetic differences were 67 haplotypes and 136 polymorphic sites in the
between: (1) individuals from the northern Arizona
combined 1111 bp of the COIII and ND3 genes (Figure 3).
populations and all other individuals and (2) individuals The mean value of nucleotide diversity was also high,
from the central Arizona populations and individuals
p ¼ 0.014 and y ¼ 0.024. These estimates are similar to
from New Mexico. The latter comparison provides an nucleotide levels observed in other small mammals
age estimate for the colonization of New Mexico by
(Nachman et al, 1994).
pocket mice, and therefore a maximum age for the
melanic populations on the three New Mexico lava
flows. Intrapopulation genetic variation: Between three and 10
To test the role of local adaptation in phenotypic mice were surveyed from each of 14 populations
evolution, we examined the extent of gene flow (Figure 1; Table 2). Highest levels of intrapopulation
occurring between each lava population and the nearest diversity were in the southern Arizona populations,
population on light-colored rocks. We made the follow- including Tinajas Altas Mts (TIN; p ¼ 0.0071), Cabeza
ing population pairwise comparisons: ARM and FRA, Prieta Mts at Tule Well (TUL; p ¼ 0.0098) and the Pinacate
KNZ and AFT, PIN and TUL, and TUM and AVR. For lava flow (PIN and MEX; p ¼ 0.0100 and 0.0120,
each of the four comparisons, we calculated FST and Nm respectively). The nucleotide diversity in populations
using DnaSP. Similarly, we measured gene flow (FST and from southern Arizona is two to 20 times higher than
Nm) between the three lava-dwelling populations in seen in any other population surveyed here. Nucleotide
New Mexico to explore the role of migration in diversity ranged from p ¼ 0.0006 to 0.0049 in other
generating melanic phenotypes on each of the lava populations. Tajima’s D was slightly negative for all
flows. For all seven pairwise comparisons, we also used populations but Carrizozo (CAR), but no values were
a coalescent-based approach implemented in the pro- significantly different from zero.

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

221

Figure 2 Spectrophotometry measurements for C. intermedius and their corresponding rock habitat for 14 populations. Color is measured as a
percent reflectance of a white standard. Localities are aligned from darkest to lightest mice. Asterisks indicate the mice and rocks that
correspond to volcanic habitats. Bars represent one standard deviation. Sample sizes are given below.

Table 2 Summary statistics of mtDNA loci by population

Population N bp Sa Nhapb p y Tajima’s D

All populations 107 1111 136 67 0.0143 0.0237 1.303

Northern Arizona
BLK* 3 1111 1 2 0.0006 0.0006
WHT 10 1111 4 5 0.0012 0.0013 0.340

Southern Arizona
TIN 9 1176 23 7 0.0071 0.0072 0.043
TUL 10 1111 38 9 0.0098 0.0121 0.902
PIN* 10 1111 39 10 0.0100 0.0124 0.949
MEX* 3 1111 20 3 0.0120 0.0120

Central Arizona
AVR 10 1176 15 5 0.0049 0.0048 0.123
TUM* 3 1111 6 2 0.0036 0.0036

Southern New Mexico

POR 5 1111 6 4 0.0023 0.0026 0.668
AFT 6 1176 3 3 0.0010 0.0011 0.447
KNZ* 10 1176 11 6 0.0021 0.0033 1.586

Central New Mexico

FRA 10 1111 9 5 0.0019 0.0029 1.443
ARM* 9 1111 12 6 0.0036 0.0040 0.446
CAR* 9 1111 7 4 0.0028 0.0023 0.831
Number of segregating nucleotide sites.
a

Number of unique haplotypes.

b

Asterisk (*) indicates habitats of volcanic origin.

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

222

Figure 3 Aligned haplotypes and polymorphic sites. The positions of the sites are indicated on the top of the table. For each site, the
consensus nucleotide is given; dots indicate identity to the consensus. The 71 haplotypes are arranged by frequency within each population,
which is given on the right. Haplotypes are unique to collecting locales with four exceptions: shared alleles in PIN-TUL, KNZ-AFT and ARM-
FRA. Site 816 is the only polymorphic site in the intergenic region and was not sequenced in all individuals.

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

223
Modeltest. The topologies were rooted with the sister
species C. penicillatus and C. bayleyi. Here, we present the
NJ topology (Figure 4). Parsimony and ML trees showed
similar topologies and shared the following features: the
two northern Arizona populations, BLK and WHT, were
reciprocally monophyletic and together formed the sister
group to all other populations, suggesting that C.
intermedius may have originated in northern Arizona.
On the other hand, all mice from New Mexico group
together as a recently derived clade. Within New Mexico,
however, none of the individual populations formed
monophyletic groups. Owing to the derived position of
these populations, this pattern may be due to incomplete
lineage sorting. Alternatively, continuous gene flow
among New Mexico populations could also produce
this pattern. Finally, individuals inhabiting lava flows
were scattered across the topology.
The multidimensional scaling (MDS) plot (Figure 5),
based on pairwise FST measures between populations,
displayed distinct patterns of population clustering and
was very similar to the clustering observed with the
phylogenetic tree (Figure 4). Each of the five geographic
regions was highly clustered, and there was a clear
demarcation between Sonoran and Chihuahuan desert
populations. Both the phylogenetic tree and the MDS
plot separated the northern Arizona populations (BLK þ
WHT) from all other populations. In addition, popula-
tions from all five geographic regions clustered indepen-
dently. Particularly strong clustering was observed
among the three populations within southern New
Mexico and among the four populations within southern
Arizona. The correlation between the genetic distance (as
measured by FST) and the two-dimensional plot was high
(r ¼ 0.94).

Population structure: Haplotype diversity was signi-

ficantly partitioned among populations and geographic
regions in C. intermedius. In the one-factor AMOVA of all
populations, most of the variation was explained by
differences among populations (FST ¼ 0.69). Additional
AMOVA analyses revealed that partitioned variation
between the two biogeographic regions, the Sonoran and
Chihuahuan deserts, was significant (FCT ¼ 0.21,
Po0.004), but that most of the variation was
partitioned among populations (FST ¼ 0.65, Po105).
Variation partitioned among the five geographic
regions was high (FCT ¼ 0.64, Po105), as was the
variation partitioned between the two northern Arizona
populations (WHT þ BLK) compared to all other popula-
tions (FCT ¼ 0.34, Po105).

Correlation between phenotype, genotype and geography

Using Mantel tests, we partitioned the data among
genetic, geographic and phenotypic components
(Table 3). There is a significant correlation between
genetic (mtDNA) variation and geographic distance,
where geographic distance explains 40% of genetic
variation (P ¼ 0.002). In addition, geographic distance
was correlated with phenotype (r ¼ 0.24), although the
Figure 3 Continued correlation between geography and genetic variation was
stronger. However, there is no correlation between
Phylogenetic relationships: Phylogenies were construc- genetic variation and phenotype (r ¼ 0.08, P ¼ 0.29), and
ted using the HKY85 þ G model of nucleotide evolution when we control for the effects of geography, the
and the estimates of ti/tv ¼ 6.05 and g ¼ 0.21 from correlation is further weakened (r ¼ 0.00, P ¼ 0.48).

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

224

Figure 4 NJ phylogeny of C. intermedius populations. Topology is rooted with C. penicillatus and C. baileyi. Asterisks indicate individuals
inhabiting lava. Geographic regions are indicated on the right. Bootstrap support is indicated under the internal branches.

Population divergence and local adaptation

To test the role of local adaptation in phenotypic
evolution, we examined the extent of gene flow
occurring between lava populations and the nearest
population inhabiting light-colored rocks. In this study,
there were four pairs of populations for which light and
melanic mice were found in close geographic proximity
(Table 4). We found substantial gene flow in several of
these comparisons. The population migration rate was
highest between PIN and TUL (Nm 460) and also high
in KNZ and AFT (Nm ¼ 8). Both of these estimates of
gene flow are higher than those that typically lead to
population differentiation. This pattern can also be
observed in the phylogenetic tree (Figure 4), where
individuals from these populations are intermingled
along the tips. However, surprisingly there was some
population structure between ARM and FRA as well as
Figure 5 MDS plot of the 14 populations based on a matrix of
TUM and AVR (Nm ¼ 5 and 2, respectively), suggesting
pairwise FST measures. Abbreviations for populations follow those that these lava populations are at least partially isolated
given in Table 1. Diagonal line separates Arizona Sonoran desert from nearby neighboring light-colored populations. In
populations and New Mexico Chihuahuan desert populations. addition, we used a MCMC method to jointly estimate

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

225
gene flow and time of divergence in order to determine if migrated from one lava flow to a neighboring (perhaps
shared polymorphisms are due to recurrent gene flow or younger) lava flow (Table 4). We found the level of gene
simply recent ancestry of populations (Table 5). In all flow was similar between CAR and KNZ compared to
four pairwise comparisons, we were able to reject the ARM and KNZ (NmB0.5). However, migration rates
hypothesis of no recurrent gene flow (M ¼ 0). In all cases, were higher between the geographically closer popula-
estimates of gene flow were lower when recent ancestry tions of CAR and ARM (Nm ¼ 1.2). Using the MCMC
was considered (two- to 10-fold), but the rank order method, we found similar patterns to those revealed by
based on Nm remained the same. FST-based statistics (Table 5). Migration rates were lower
In addition, we estimated levels of gene flow connect- between KNZ and both ARM (Nm ¼ 0.14) and CAR
ing lava-dwelling populations to examine the hypothesis (Nm ¼ 0.02). Migration rates were higher between CAR
that melanic mice evolved once in New Mexico and and ARM (Nm ¼ 0.40), and we were able to reject a
model of recent divergence and no gene flow (M ¼ 0)
between these two populations. All MCMC estimates of
gene flow were lower than those based on FST. However,
Table 3 Mantel and partial Mantel test for significant correlation we noted that 95% credibility estimates are very large
between genetic, geographic and phenotypic variation in 14 and encompass FST estimates of migration rate.
populations
Matrix comparison r P-value
Discussion
MtDNA–Geography 0.40 0.002
MtDNA–Phenotype 0.08 0.285 We documented substantial phenotypic variation across
Geography–Phenotype 0.24 0.030 the range of C. intermedius. This color variation was
MtDNA–Geography_(Phenotype) 0.40 0.001 not, however, correlated with phylogeny, suggesting
MtDNA–Phenotype_(Geography) 0.00 0.476 that history is not responsible for the present distribu-
Genetic distances between populations were measured using tion of phenotypic variation. Although there is sub-
pairwise FST. Geographic distances were based on the central stantial genetic structure between geographic regions,
collecting site for each population. Phenotypic distance (Manhattan high levels of gene flow (or recent ancestry) connect
distance) was measured using the program NSYS. Significance populations within a region. Despite these high levels of
values were determined by comparing the observed z-statistic and
the expected z-statistic, generated by a randomized distribution
local gene flow, color variation is strongly correlated with
from 1000 permutations. Parentheses indicate which factor is habitat color in most populations, suggesting that natural
removed from the correlation. selection for substrate matching is strong in this species.

Table 4 Pairwise estimates of gene flow between the four paired dark and light rock populations and separately the three New Mexico lava
populations
Population comparison FST Nm Distancea (km)

Lava population and neighboring light rock population

PIN–TUL 0.008 60.43 16
TUM–AVR 0.172 2.41 28
KNZ–AFT 0.057 8.35 7
ARM–FRA 0.098 4.61 7

New Mexico lava populations harboring melanic mice

CAR–ARM 0.296 1.19 118
KNZ–CAR 0.510 0.49 222
KNZ–ARM 0.517 0.47 158
Approximate geographic distances, based on the center of each collecting area, are given for each pairwise comparison.
a

Table 5 Pairwise estimates of female effective population size (y ¼ 2Nfu), female migration rate (M ¼ Nfm) and divergence time (T ¼ t/Nf),
based on analysis of two mtDNA genes (COIII and ND3) using the program MDIV (Nielsen and Wakeley, 2001)
Population comparison y M T

Lava population and neighboring light rock population

PIN–TUL 13.49 (9.38–23.09) 8.12 (3.24–19.08) 0.04
TUM–AVR 2.36 (1.53–6.25) 1.38 (0.58–9.14) 0.32
KNZ–AFT 2.65 (0.96–5.78) 1.50 (0.68–9.42) 0.26
ARM–FRA 3.04 (1.85–6.28) 1.46 (0.08–6.58) 0.98

New Mexico lava populations harboring melanic mice

CAR–ARM 3.65 (2.12–7.26) 0.40 (0.04–6.06) 0.60
KNZ–CAR 3.16 (1.99–6.85) 0.02 (0.00–1.16) 1.30
KNZ–ARM 4.61 (2.09–9.48) 0.14 (0.04–1.84) 1.10
A uniform prior (0,10) was assumed for y, M and T. Markov chain was simulated for 5  106 generations, where 5  105 generations were used
for a burn-in time. For y and M, values with the highest likelihood scores are given, followed by 95% credibility intervals within parentheses.

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

226
Adaptive phenotypic variation the northern populations and the remaining populations.
Spectrophotometric measurements reveal a strong corre- These results suggest that C. intermedius may have
lation between substrate color and the dorsal pelage of originated in northern Arizona and expanded to
rock pocket mice (R2 ¼ 0.43). Rock substrate from lava southern Arizona and eventually east into southern
flows show significantly lower reflectance than rocks and central New Mexico. The average pairwise
from other areas. The dorsal pelage of mice inhabiting divergence between these two northern Arizona
these lava flows has correspondingly lower reflectance. populations and the remaining populations was 0.035.
One distinct exception is the population of Black Tank Assuming a 2% per million year clock, C. intermedius
(BLK) in northern Arizona in which the rock has similar may have expanded from northern Arizona southward
reflectance to other lava flows but the mice are similar in about 1.5 million years ago. Interestingly, this
reflectance to other nonlava mice. This anomaly likely corresponds to the age of the Pinacate lava flow,
reflects the fact that the BLK lava flow is an extension of suggesting that pocket mice were likely present in the
the relatively young Sunset Crater (o800 years old); southern Arizona region during the formation of the
consequently, there may have been insufficient time for Pinacate lava flow.
melanic mice to evolve in this population. In addition,
there are no populations of melanic mice nearby from Correlation between phenotype, genotype and geography
which melanic migrants may invade. Are morphological patterns reflected in the genetic
In other cases, migration from neighboring popula- structure of C. intermedius? In this case, removing the
tions may be responsible for the presence of melanic (or geographic component of the phenotypic variation
dark) phenotypes in some areas. For example, the among populations resulted in no correlation between
Tumamoc Hill population, which occurs on olivine phylogeny and color variation (r ¼ 0.00). The phenotypic
basalt, is separated by just a few miles from a melanic patterns in C. intermedius may reflect processes operating
population of mice (Dice and Blossom, 1937), which on a spatial or temporal scale much smaller than that
occurs on the extremely dark basalt of Black Mountain, reflected in broad-scale geographic patterns, resulting in
which we were unable to sample. Owing to their close a poor correlation between color patterns and geography
geographic proximity, it is likely that gene flow occurs and phylogeny. This result suggests that adaptation to
between these two populations. In fact, Blossom (1931) local environments (ie natural selection) is a key force
described a dark-colored race ‘nigrimontis’ from Black driving morphological diversity in this system and that
Mountain and Tumamoc Hill. Therefore, Black Mountain historical contingency plays a relatively small role.
may be the source of the Tumamoc Hill mice and may
explain why the Tumamoc mice are significantly darker
than their habitat. Using genetic data, we also explored Local adaptation
migration of melanic mice among the lava-dwelling New To determine the amount of gene flow occurring between
Mexico populations (see below). populations that differed in their substrate color (and
phenotype), we calculated gene flow in four comparisons
Genetic variation between geographically proximate populations inhabit-
In the northern range of C. intermedius, we found a high ing lava and light rock. We found that in each case
level of genetic partitioning between the five geographic substantial gene flow occurs over short distances (Nm
regions (FCT ¼ 0.64). This pattern is also reflected in the ranged from 60.5 to 2.5 using FST-based statistics and
phylogenetic tree, which shows a general grouping of 8.12–1.46 using MCMC estimates), suggesting that
individuals by geographical region, although some of natural selection must be strong in order to maintain
these groups are not monophyletic (Figure 4). The MDS habitat-specific color patterns. In the Pinacate region,
plot, based on pairwise FST estimates, also reveals strong Hoekstra et al (2004) estimated selection coefficients
clustering by geographic region (Figure 5). As rock against light mice inhabiting lava as high as 0.39. It is
pocket mice exclusively inhabit rocky areas and are important to recognize, however, that our estimates of
replaced by sand-dwelling species (eg C. penicillatus) in gene flow from FST are based on a model of migration–
nonrocky areas, their habitat is largely discontinuous drift balance at equilibrium. Inherent in these calcula-
throughout its range, which may underlie the population tions are a number of assumptions which may not be
structure observed at this geographic scale. At a smaller biologically realistic (Whitlock and McCauley, 1999).
scale, however, substantial gene flow sometimes occurs
between neighboring populations (Table 4). Evolution of melanism in New Mexico popula-
It is important to note that the genetic analyses we tions: Molecular analysis of pigmentation genes shows
performed are based on a single molecular marker, that melanic pocket mouse populations have evolved
mtDNA. In addition, mtDNA tracks only female migra- independently in the Pinacate population of Arizona and
tion. However, because there is no evidence for sex- the lava populations of New Mexico (Hoekstra and
biased dispersal in heteromyid rodents, mtDNA may be Nachman, 2003). However, it remains unclear if melanic
an accurate predictor of the average gene flow for both populations have evolved repeatedly within New
sexes (Jones, 1993). Mexico, that is, independently in Armendaris (ARM),
Carrizozo (CAR) and Kenzin (KNZ) populations. Our
Biogeographical history of C. intermedius: Phylogenetic analysis suggests that substantial gene flow may occur
analysis reveals a strong split between the two northern among these lava-dwelling pocket mouse populations
populations (BLK and WHT) and the rest of the in New Mexico, raising the possibility that melanic
populations. This northern clade is basal in the mice evolved once in New Mexico and migrated to other
phylogeny. AMOVA analysis suggests that 34% of lava populations. Again, however, we caution that
genetic variation was explained by differences between the assumption of migration–drift equilibrium is

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

227
unlikely to be met, especially among the New Mexico Together, our results suggest that strong selection is
populations. These populations are relatively young likely to be the driving force in promoting morphological
(see below) and at least some of the shared variation diversity in this system. We observe high levels of gene
among populations may reflect ancestral polymorphism flow in several neighboring populations that differ in
for which lineage sorting is incomplete, thus leading to substrate color, and the observation that mice closely
an overestimate of the true level of gene flow. match their substrate despite this gene flow indicates
The MCMC method of Nielsen and Wakely (2001), that strong selection is maintaining habitat-specific
however, assesses the relative roles of migration and phenotypes. Migration is commonly thought to be a
isolation as causes of the observed differentiation homogenizing force impeding local adaptation; however,
between populations. here we raise the possibility that migration between lava
We estimated the approximate time when pocket mice flows may also promote local adaptation by introducing
invaded New Mexico from southern Arizona in order to beneficial alleles into neighboring populations experien-
identify the maximum time for adaptation to local cing similar selective regimes.
conditions to occur in New Mexico. The average
pairwise difference between the New Mexico popula-
tions and the central Arizona populations was 0.0095. Acknowledgements
Assuming a molecular clock of 2% mtDNA divergence
per million years, these mice first appeared in New We thank K Drumm, B Haeck, J Kim, V Klein, A
Mexico roughly 500 000 years ago. The ages of both the Kurosaki, A Litt and J Storz for assistance in the field.
Kenzin (B500 000 years old) and Armendaris (B750 000 Comments from the editor and two anonymous re-
years old) lava flows are similar to or older than this viewers greatly improved the manuscript. A Redd
estimate. Thus, the maximum time for local adaptation provided valuable assistance with data analysis. We also
on these lava flows may be limited by the immigration of thank Vergial Harp of the Cabeza Prieta National
C. intermedius to New Mexico approximately 500 000 Wildlife Refuge and Thomas Waddell of the Armendaris
years ago. Ranch, Turner Enterprises Inc. for access to field sites.
One lava flow in New Mexico, Carrizozo, is less than Thanks to K Drumm and J Kim for generating some of
1000 years old, yet harbors a population of mice which the molecular data. This work was supported by an NIH
have uniformly melanic dorsal pelage (Dice and Blos- NRSA Postdoctoral Fellowship (HEH) and an NSF grant
som, 1937). Given the young age of the lava flow, this is (MWN).
somewhat surprising and leads us to consider the
relative probabilities that these melanic mice arose from
a new mutation or from migration of a new allele. We can References
get a rough estimate of these probabilities as follows. Benson SB (1933). Concealing coloration among some desert
Mutation rates are likely to be approximately 105 per rodents of the southwestern United States. Univ Calif Publ
locus per generation (eg Schlager and Dickie, 1971), Zool 40: 1–69.
corresponding to about 108 per nucleotide site per Blossom PM (1931). Relation between color of desert rodents
generation. If there are 100 sites in the genome at which and of the soil. Carnegie Inst Washington Yearb 30: 266.
mutations can produce dark color, the overall mutation Castellano S, Balletto E (2002). Is the partial Mantel test
rate to dark color is about 106 per genome per inadequate? Evolution 56: 1871–1873.
generation. If we assume that the population size is Dice L, Blossom PM (1937). Studies of mammalian ecology in
southwestern North America, with special attention to the
approximately 104 and that 103 generations have elapsed
colors of desert mammals. Publ Carnegie Inst Washington 485:
since the appearance of the lava flow (one generation per 1–25.
year), then there has been sufficient time for approxi- Excoffier L, Smouse PE, Quattro JM (1992). Analysis of
mately 10 darkening mutations to appear. In contrast, if molecular variance inferred from metric distances among
the number of migrants (Nm) from the nearest dark DNA haplotypes: application to human mitochondrial DNA
population (Armendaris) is approximately one per restriction data. Genetics 131: 479–491.
generation (FST-based estimate; Table 4) or one indivi- Haldane JBS (1948). The theory of a cline. J Genet 48: 277–284.
dual every two generations (MCMC-based estimate; Hoekstra HE, Drumm KE, Nachman MW (2004). Ecological
Table 5), then there has been sufficient time for genetics of coat color variation in pocket mice: geographic
approximately 500–1000 dark alleles to be introduced variation in selected and neutral genes. Evolution 58: 1329–
1341.
by migration. These rough calculations suggest that Hoekstra HE, Nachman MW (2003). Different genes underlie
migration is a more likely source of dark alleles in the adaptive melanism in different populations of rock pocket
Carrizozo population than is mutation. mice. Mol Ecol 12: 1185–1194.
Migration rates are likely to be high between Armen- Jones T (1993). The social systems of heteromyid rodents. In:
daris and Carrizozo because there is suitable rocky Genoways HH, Brown JH (eds) Biology of the Heteromyidae.
habitat (the San Andres Mts), which spans much of the American Society of Mammalogists, pp 575–595.
area between these lava flows. Additionally, melanic Lenormand T (2002). Gene flow and the limits to natural
mice may also have migrated from the Kenzin lava selection. Trends Ecol Evol 17: 183–189.
(although at a lower rate), increasing the likelihood of Manly BFJ (1986). Randomization and regression methods for
testing association with geographical, environmental and
introducing melanic alleles via migration. While it is
biological distances between populations. Res Popul Ecol 28:
more likely that melanic migrants are responsible for the 201–218.
melanic population at Carrizozo versus new mutations Nachman MW, Boyer SN, Searle JB, Aquadro CF (1994).
in the Carrizozo population, identifying the genetic basis Mitochondrial DNA variation and the evolution of Robert-
of melanism in these populations will directly address sonian chromosomal races of house mice, Mus domesticus.
this question. Genetics 136: 1105–1120.

Heredity
 Local adaptation in pocket mice
HE Hoekstra et al

228
Nachman MW, Hoekstra HE, D’Agostino S (2003). The genetic Slatkin M (1993). Isolation by distance in equilibrium and non-
basis of adaptive melanism in pocket mice. Proc Natl Acad Sci equilibrium populations. Evolution 43: 264–279.
USA 100: 5268–5273. Swofford DL (1999). Phylogenetic Analysis using Parsimony
Nei M, Li WH (1979). Mathematical model for studying genetic (PAUP). Sinauer: Sunderland, MA.
variation in terms of restriction endonucleases. Proc Natl Tajima F (1989). Statistical method for testing the neutral
Acad Sci USA 76: 5269–5273. mutation hypothesis by DNA polymorphism. Genetics 123:
Nielsen R, Wakeley J (2001). Distinguishing migration from 585–595.
isolation: a Markov chain Monte Carlo approach. Genetics Travisano M, Mongold J, Bennett A, Lenski R (1995). Experi-
158: 885–896. mental tests of the roles of adaptation, chance, and history in
Palsboll PJ, Berube M, Aguilar A, Notarbartolo-Di-Sciara G, evolution. Science 267: 87–90.
Nielsen R (2004). Discerning between recurrent gene flow Waterson GA (1975). On the number of segregating sites in
and recent divergence under a finite-site mutation model genetical models without recombination. Theoret Popul Biol 7:
applied to North Atlantic and Mediterranean sea fin whale 256–276.
(Balaenoptera physalus) populations. Evolution 58: 670–675. Weckerly WF (1983). Geographic Variation of the Rock Pocket
Posada D, Crandall KA (1998). MODELTEST: testing the model Mouse, Perognathus intermedius Merriam, on the Perdro Armen-
of DNA substitution. Bioinformatics 14: 817–818. dariz Lava Flow of South-Central New Mexico. Master’s Thesis.
Raufaste N, Rousset F (2001). Are partial Mantel tests adequate? Eastern New Mexico University, Portales, NM.
Evolution 55: 1703–1705. Whitlock M, McCauley D (1999). Indirect measures of gene flow
Rousset F (1997). Genetic differentiation and estimation of gene and migration: F-ST not equal 1/(4Nm+1). Heredity 82: 117–
flow from F-statistics under isolation by distance. Genetics 125.
145: 1219–1228. Wilson AC, Cann RL, Carr SM, George M, Gyllensten UB,
Rousset F (2002). Partial Mantel tests: reply to Castellano and Helmbychowski KM et al (1985). Mitochondrial DNA and
Balletto. Evolution 56: 1874–1875. two perspectives on evolutionary genetics. Biol J Linn Soc 26:
Rozas J, Rozas R (1999). DnaSP version 3: an integrated 375–400.
program for molecular population genetics and molecular
evolution analysis. Bioinformatics 15: 174–175.
Schlager G, Dickie M (1971). Natural mutation rates in the
house mouse. Estimates for five specific loci and dominant
mutations. Mutat Res 11: 89–96. Appendix A1
Schneider S, Roessli D, Excoffier L (2000). Arlequin: A Software
Program for Population Genetics Data Analysis. Genetics and Previously published mtDNA sequences from three sites
Biometry Lab, Department of Anthropology, University of (Hoekstra et al, 2004) are as follows: PIN – AY648411,
Geneva. AY648511, AY648514-9, AY648525-6; TUL – AY648417-8,
Slatkin M (1985). Gene flow in natural populations. Annu Rev AY648420-1, AY648527-8, AY648533, AY648537,
Ecol Syst 16: 393–430. AY648544, AY648552; MEX – AY648479-80, AY648503.

Heredity