DISSECTING THE GEOGRAPHICAL EXPANSION OF PEROMYSCUS LEUCOPUS IN THE
NORTHERN GREAT LAKES: INSIGHTS FROM GENETICS, MORPHOMETRICS AND
ECOLOGY

By

Rosa Anna Moscarella

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Zoology and Ecology, Evolutionary Biology and Behavior

2011

ABSTRACT

DISSECTING THE GEOGRAPHICAL EXPANSION OF PEROMYSCUS LEUCOPUS IN THE
NORTHERN GREAT LAKES: INSIGHTS FROM GENETICS, MORPHOMETRICS AND
ECOLOGY

By

Rosa Anna Moscarella

Geographical and ecological distributions of biological species are primarily limited by
prevailing climate conditions, and thus it is hardly surprising that recent rapid climate change has
triggered evolutionary processes that are deeply affecting the biology of many organisms. One
consequence of rapid climate change has been the occurrence of shifts in geographical ranges of
many species, such as the white-footed mouse (Peromyscus leucopus), which has experienced a
northern range expansion in Michigan and Wisconsin over the last few decades and as a
consequence is now common in localities where it was absent until recently. In this study, I
sought to understand some of the historical causes and consequences of the expansion process
undergone by this rodent species. I analyzed complete D-loop sequences from 595 mice
collected throughout the Great Lakes region to investigate the origin of the newly established
populations. Phylogeographic analyses strongly support two lineages, spanning populations from
Wisconsin through the western Upper Peninsula (UP), and from the Lower Peninsula (LP) to the
eastern Upper Peninsula (UP) of Michigan, respectively. With few exceptions, western
populations in the UP originate from Wisconsin while the eastern UP populations originate from
LP. Both lineages display a genetic signature of spatial expansion, whereas demographical
expansion is apparent only for the western UP-Wisconsin lineage. I also investigated whether the

rapid geographic expansion of populations of this species in the northern Great Lakes region has
been accompanied by rapid morphological divergence, and whether the pattern of divergence
suggests a neutral process or reflects the influence of non-random factors. Using morphometric
analyses of mandible shape, I found significant differences between recently expanded and long
established populations not only in the magnitude of shape variation but also in the
dimensionality of variation, which suggest that the expansion process has been accompanied by
rapid morphological change along directions of the phenotype space that are absent in ancestral
populations. Furthermore, phenotypic and phylogeographic patterns are incongruent in these
populations, suggesting that changes in morphology have not been caused by random phenotypic
drift. To further explore possible causes for the phenotypic divergence observed between old and
recent populations, I examined the association between geographical and environmental
variables and patterns of genetic and morphometric variation. Results from these analyses
suggest that even though most factors analyzed contribute to the global variation in the mandible,
cold temperatures and snowfall magnitudes appear to be the best predictors for the differentiation
between old and new populations.

Copyright by
ROSA ANNA MOSCARELLA
2011

To my parents, Emilia y Gaetano, because all I am is because of them, literally
To my best friend, my soul mate, and the love of my life: Eladio
To my daughter Veronica

v

ACKNOWLEDGMENTS

Because getting Ph.D. is not only to complete a Dissertation, I want to start my
acknowledgments with my dear husband Eladio, who has been my life coach and has given me
his love and support unconditionally. My daughter Veronica, who is the engine of my life. My
family, for their support and encouragement. I have met amazing people and have great friends,
and thanks to them I was able to keep my sanity over these years: thanks Lili, Amy, Romari,
Juliana, Martha, Karina, Nena, Michelle, Sarah, Jessica, Sandra, Ananias, Dana, Laura, Laurita,
Tina, Taty, Terri, Lauri, Emily, Carola, Fabiola, Adriana, Dr. Susan Hill, April, Mark UrbanLurain, Miriam, Don, the ZOL328 gang, David Houle and the Houligangs, my therapist … To
my advisor, Barbara Lundrigan, and committee members, Bryan Epperson, Philip Myers,
Douglas Schemske, Kim Scribner, for all their advices and suggestions that helped to improve
my Dissertation. To the following institutions for providing samples or giving a place where to
work: Michigan State Museum, Wisconsin Department of Natural Resources, University of
Wisconsin Museum of Vertebrates, Hoffman’s Lab at the Miami University in Ohio, the
University of Michigan Museum of Zoology, Duda’s Lab and the Genomic Diversity Lab at
UMMZ. To all the people that helped me in one way or another: Sam Márquez, Bryan Carstens,
Tom Duda, Smruti Damania, Kate Johnson, Lucía Luna, David Mindell, Allison Poor, Pamela
Roy, Mary Catherine Steinwald, Priscilla Tucker, Zac Taylor, Chris Yahnke, Loren Ayers, Susan
Hoffman, the MSU Zoology Department Staff, Biological Science Program Staff, DSME staff,
Museum staff, and the UMMZ staff. Thanks to MSN Messenger, Skype, facebook, Google
Voice, for allowing me to stay in touch with my loved ones. This Dissertation was partially
funded by Michigan State University Department of Zoology, EEBB Graduate Program,

vi

Graduate School, College of Natural Science, and Office of the Vice President for Research and
Graduate Studies; University of Michigan Museum of Zoology; Sigma Xi; American Society of
Mammalogists; and the John R. Shaver Research Fellowship.

vii

TABLE OF CONTENTS

LIST OF TABLES

x

LIST OF FIGURES

xiv

CHAPTER I: INTRODUCTION
Overview
Geographical expansion of Peromyscus leucopus in the northern Great Lakes
region
Northern Great Lakes Ecoregions
Objectives and Outline
Figures
Bibliography

3
5
6
9
11

CHAPTER II: PHYLOGEOGRAPHY, GENETIC STRUCTURE AND HISTORICAL
DEMOGRAPHY OF EXPANDING POPULATIONS OF PEROMYSCUS
LEUCOPUS IN THE NORTHERN GREAT LAKES
Abstract
Introduction
Materials and Methods
Results
Discussion
Tables
Figures
Bibliography

15
15
15
19
27
31
41
49
54

CHAPTER III: MORPHOMETRIC DIFFERENTIATION OF EXPANDING
POPULATIONS OF PEROMYSCUS LEUCOPUS IN THE NORTHERN GREAT
LAKES: HISTORICAL AND GENETIC FACTORS
Abstract
Introduction
Materials and Methods
Results
Discussion
Tables
Figures
Appendix
Bibliography

63
63
63
69
80
84
95
99
109
130

viii

1
1

CHAPTER IV: MORPHOMETRIC DIFFERENTIATION OF EXPANDING
POPULATIONS OF PEROMYSCUS LEUCOPUS IN THE NORTHERN GREAT
LAKES: ENVIRONMENTAL AND GEOGRAPHICAL FACTORS
Abstract
Introduction
Materials and Methods
Results
Discussion
Conclusions
Tables
Figures
Appendix
Bibliography

ix

140
140
140
143
149
155
163
166
174
195
202

LIST OF TABLES

Table II.1. Origin of the samples analyzed in this study, where n is the sample size, h is
the number of haplotypes observed, H is the haplotype diversity, S is the number of
segregating sites over the complete sequence, Π is the nucleotide diversity, and Lat and
Long are the decimal geographic coordinates of the geographic midpoint of each locality.
NA indicates that these parameters were not calculated for those populations because of
small sample size

41

Table II.2. Spatial population structure of Pln in the Northern Great Lakes, as suggested
by spatial analysis of molecular variation (SAMOVA). This analysis was run assuming K
= 2, 3 and 4 groups

45

Table II.3. Demographic parameters estimated using the isolation with migration model
implemented in IMa. 1, 2, and a are the mutation rates of population 1 (LP-EUP),
population 2 (WUP-WI), and the ancestral population, respectively; N1, N2, and N3 are
the effective population sizes of females (number of individuals) of LP-EUP, WUP-WI,
and the ancestral population, respectively; m1 and m2 are the migration rates per
generation into population 1 from population 2, and into population 2 from population 1,
respectively; and t is the time since the two groups split, in years. Conversions of model
parameters to demographic parameters were done assuming 0.5 years/generation and a
mutation rate of 1.2202*10-4 substitutions/ locus/year. HiPt is the value with the highest
frequency (i.e. mode), and HPD90Lo and HPD90Hi are the lower and upper bounds of the
estimated 90% Highest Posterior Density Interval

46

Table II.4. Comparison of nested models to the full-parameter model of isolation with
migration, using IMa. 1, 2, and a are the mutation rates for population 1 (LP-EUP),
population 2 (WUP-WI), and the ancestral population, respectively; m1 and m2 are the
migration rates per generation into population 1 from population 2, and into population 2
from population 1, respectively

47

Table II.5. Estimated parameters of the Spatial and Sudden Expansion models, and results
of tests of goodness of fit of observed data to each model for each group suggested by
SAMOVA, and for the whole sample. M is the number of migrants per generation, τ is the
time to expansion, in generations; 0 and 1 are the population mutation rates before and after
expansion, respectively; r is Harpending’s raggedness index, and P is the goodness of fit to
each model

48

x

Table III.1. Results from a nested MANOVA design measuring the effects of lineage (i.e.,
LP-EUP and WUP-WI), expansion history (i.e., old, established, vs. new, recently invaded
populations) nested within lineage, and populations, nested within old and new localities,
on shape of the mandible in Peromyscus leucopus noveboracensis. Note that Effect Size is
not additive

95

Table III.2. Results from a nested ANOVA design measuring the effects of lineage,
expansion history (i.e., old, established, vs. new, recently invaded populations) nested
within lineage, and populations, nested within old and new localities, on centroid size of
the mandible in Peromyscus leucopus noveboracensis. Note that Effect Size is not
additive

96

Table III.3. Contribution of mandible developmental blocks to the difference between old
and new populations in each of the two lineages sampled in this study. Contributions are
given in squared Procrustes units and as percentages of the total squared Procrustes
distance between the mean shape of old and new populations (in parenthesis). Distances
have been multiplied by 105

97

Table III.4. Mantel test to assess the association between morphometric and genetic
distances. Significance level corresponding to a 95% confidence interval after a
Bonferroni’s correction (α = 0.0125)

98

Table A.III.1. List of specimens used in this study, including collection information
(catalog number, museum), locality of origin (county and region), lineage, timing of
population establishment, sex (0=males, 1= females), and coordinates of geographical
midpoint of the locality of origin. LP=Michigan’s Lower Peninsula, UP=Michigan’s
Upper Peninsula, WI=Wisconsin

110

Table A.III.2. Matrix of pairwise genetic distances (Nei’s D) between populations
included in this study. Above diagonal: Average number of pairwise differences between
populations (PiXY). Diagonal elements (in bold): Average number of pairwise
differences within population (PiX). Below diagonal: Corrected average pairwise
difference (PiXY-(PiX+PiY)/2)
124

Table A.III.3. Matrix of pairwise Fst values between populations included in this study

Table A.III.4. Matrix of pairwise morphometric distances between populations included in
xi

126

this study. Distances were calculated as pairwise Euclidean distances between population
means computed from the shape variables derived from projecting Procrustes residuals
onto the first 46 eigenvectors of their covariance matrix. Values in each cell were multiply
by 103

128

Table IV.1. Relationship between mandible shape distances and geographical (Geo),
ecological (Eco), and genetic (Fst) distances based on partial Mantel tests. See text for
details about how distances were calculated. The correlation (r) is between shape distance
and the corresponding distance variable after controlling for the other two factors. Tests of
significance were based on 1,000 permutations. The Bonferroni-adjusted alpha level is
0.006. Significant correlations are indicated in boldface

166

Table IV.2. Relationship between mandible size distances and geographical (Geo),
ecological (Eco), and genetic (Fst) distances base on partial Mantel tests. See text for
details about how distances were calculated. The correlation (r) is between size distance
and the corresponding distance variable after controlling for the other two factors.
Significance tests were based on 1,000 permutations. The Bonferroni-adjusted alpha level
is 0.006. Significant correlations are indicated in boldface

167

Table IV.3. Relationship between genetic distances based on pairwise Fst and
geographical (Geo) and ecological (Eco) distances base on partial Mantel tests. . See text
for details about how distances were calculated. The correlation (r) is between genetic
distance and the corresponding distance variable after controlling for the other factor.
Significance tests were based on 1,000 permutations. The Bonferroni-adjusted alpha level
is 0.0083. Significant correlations are indicated in boldface.

168

Table IV.4. Results of a partial least-square analysis of the covariation between mean
mandible shape and temperature for populations of Peromyscus leucopus in the LP-EUP
and WUP-WI

169

Table IV.5. Results of a partial least-square analysis of the covariation between mean
mandible shape and precipitation for populations of Peromyscus leucopus in the LP-EUP
and WUP-WI

170

Table IV.6. Results of a partial least-square analysis of the covariation between mean
mandible shape and geographical coordinates for populations of Peromyscus leucopus in
the LP-EUP and WUP-WI

171

Table IV.7. Results of a partial least-square analysis of the covariation between mean
xii

mandible shape and land cover/use for populations of Peromyscus leucopus in the LPEUP and WUP-WI

172

Table A.IV.1. Temperature and precipitation variables and coordinates of the geographical
midpoint for all populations in this study. Temperatures are in Fahrenheit degrees; snow
and rain precipitations in mm2. In Timing, O and N stand for old and new, respectively. l
is Simpson’s Index of Dominance

196

Table A.IV.2. Matrix of pairwise ecological distances, (1-IM), were IM is Morista’s Index
of community similarity

198

Table A.IV.3. Matrix of pairwise geographic distances, in km

200

xiii

LIST OF FIGURES

Figure I.1. Distribution of Peromyscus leucopus, according to Lackey et al. (1985). The
numbers indicate the distribution of the 17 sub-species proposed for this mouse in North
America. The grey area with the number 14 shows the range of the subspecies P. l.
noveboracensis. Note that in this map this subspecies is not reported for the Upper
Peninsula of Michigan.

9

Figure I.2. Distribution of the white-footed mouse, Peromyscus leucopus, in Michigan
during a) 1883-1980 and b) 1981-2006 (Myers et al. 2009)

10

Figure II.1. Geographic location of Peromyscus leucopus populations analyzed in this
study. The code for each population is as in Table II.1. LI = Livingstone, Michigan; IL =
Johnson, Illinois; OH = Butler, Ohio. Black circles indicate populations in the former
northern limit of the species in the region. Grey circles indicate populations in newly
invaded localities

49

Figure II.2. Neighbor Joining tree showing the haplotype genealogy of Peromyscus
leucopus in Wisconsin (WI), and Michigan’s Lower Peninsula (LP) and Upper Peninsula
(UP). Peromyscus maniculatus is the outgroup. Some internal clades have been magnified
for better visualization. ―+‖ indicate haplotypes from Chippewa County (CW) that cluster
with the WUP-WI clade and ―*‖ haplotypes that are not in the expected clade given
collection locality. Bootstrap support values are given for the main nodes. Acronyms as in
Table II.1; LI = Livingstone, Michigan; IL = Johnson, Illinois; OH = Butler, Ohio

50

Figure II.3. Results from Mantel’s tests evaluating the correlation between genetic
(Φst/1-Φst) and geographical (km) distances between populations. a) Population
comparisons within the LP-EUP group. c) Population comparisons within the WUP-WI.
Correlation coefficients and P values are indicated in the text

51

Figure II.4. Interconnectivity among populations based on pairwise genetic distances.
The color of the line indicates genetic distance between population pairs, based on
pairwise Φst/1-Φst values, with dark blue indicating the most similar genetically and dark
red the most different. On the left are represented populations with genetic differences up
to 35%, and on the right populations with genetic differences > 35%. For interpretation of
the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation

52

xiv

Figure II.5. Distribution of pairwise differences of the observed data (dark grey) and the
simulated data (light grey) under a model of stepwise (spatial) expansion (left) and sudden
(demographic) expansion (right) for a and b (LP-EUP), c and d (WUP-WI), and e and f
(all data). Results of tests of goodness of fit of the data to each model are shown in Table
II.5

53

Figure III.1. Configuration of landmarks (white circles) and semi-landmarks (black dots)
sampled in this study. Shadowed areas indicate regions with developmental or functional
significance. With the exception of the masseteric ridge, these regions are defined to
roughly correspond to developmental units of the mandible (Atchley and Hall 1991)

99

Figure III.2. Geographic origin of the 24 populations examined in this research. The twoletter acronyms are as in Table A.III.1. Black circles represent old populations, while grey
ones are new populations, according to the criterion explained in the text

100

Figure III.3. Plots of scores of PC1 vs. PC2 (top) and PC3 vs. PC4 (bottom). Black
circles are old populations in LP-EUP and white circle is the new population in this
lineage. Black and white triangles are old and new populations, respectively, in WUP-WI.
The two-letter acronyms are as in Table A.III.1

101

Figure III.4. Shape deformations implied by first four PCs of population means: a) PC1
(26.32%), b) PC2 (23.57%), c) PC3 (11.17%), d) PC4 (7.78%). Colors represent
continuous change in surface area of infinitesimally local regions of the mandible, as
inferred via TPS interpolation, where colder colors (toward blue) are interpreted as local
contraction and warmer ones (toward red) as expansion. Numbers in the scale represent
percentage of local change (e.g., -0.1 means a local contraction of 10% with respect to the
reference). Deformation vectors and grids (but not color patterns) have been magnified x5
to improve visualization

102

Figure III.5. Box plots of canonical scores for all populations, sorted by lineage (i.e., LPEUP, WUP-WI). Old populations are indicated by a shadowed background and new
populations are on white (see text for details). The number of axes chosen for
interpretation was based on Bartlett’s test results, which estimated four significant
dimensions for the differences among population means. Proportion of the variation
explained by these canonical axes was CV1: 16.40%, CV2: 11.17%, CV3: 10.19%, and
CV4: 8.67%

103

Figure III.6. Differences between mean shapes of old and new populations within (a) LPxv

EUP and (b) WUP-WI lineages. Landmarks in mean shape reference of old samples are
connected by a black line; landmarks in mean shape of new population are connected by a
grey line. Note the difference in scale, showing more pronounced changes in (a). Colors
represent proportional change throughout the mandible, where colder colors (toward blue)
are interpreted as local contraction and warmer ones (toward red) as expansion. Numbers
in the scale are percentage of local change (e.g., -0.10 means a local contraction of 10%
with respect to the reference). To improve visibility of transformation grids, deformations
have been multiplied by 4 in (a) and by 8 in (b)

104

Figure III.7. Trend of morphological variation within old populations and from old to
new populations. The most extreme specimens are shown in each case to facilitate
visualization of the changes already suggested by deformation grids (Fig. III.6)

105

Figure III.8. UPGMA dendrogram based on morphometric distances among populations.
Distances were calculated as pairwise Euclidean distances between population means
computed from the shape variables derived from projecting Procrustes residuals onto the
first 46 eigenvectors of their covariance matrix

106

Figure III.9. UPGMA dendrogram based on pairwise Nei’s genetic distances among
populations

107

Figure III.10. UPGMA dendrogram based on pairwise Fst values among populations

108

Figure IV.1. Land Cover/Use of populations sampled in this study. Each circle covers a
radius of 454 m from the geographical midpoint of each locality. See text for additional
information

174

Figure IV.2. UPGMA dendrograms based on morphometric, genetic, and ecological
distances in LP-EUP and WUP-WI lineages

175

Figure IV.3a. Correlation between temperature and shape on PLS1 for LP-EUP, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

181

Figure IV.3b. Correlation between temperature and shape on PLS1 for WUP-WI, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

182

Figure IV.3c. Correlation between temperature and shape on PLS2 for WUP-WI, and the

183

xvi

corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

Figure IV.4a. Correlation between precipitation and shape on PLS1 for LP-EUP, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

184

Figure IV.4b. Correlation between precipitation and shape on PLS1 for WUP-WI, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

185

Figure IV.5a. Correlation between latitude and shape on PLS1 for LP-EUP, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

186

Figure IV.5b. Correlation between longitude and shape on PLS1 for WUP-WI, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

187

Figure IV.5c. Correlation between latitude and shape on PLS2 for WUP-WI, and the
corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

188

Figure IV.6a. Correlation between land cover/use and shape on PLS1 for LP-EUP, and
the corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

189

Figure IV.6b. Correlation between land cover/use and shape on PLS2 for LP-EUP, and
the corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

190

Figure IV.6c. Correlation between land cover/use and shape on PLS1 for WUP-WI, and
the corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

191

Figure IV.6d. Correlation between land cover/use and shape on PLS2 for WUP-WI, and
the corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

192

Figure IV.6e. Correlation between land cover/use and shape on PLS3 for WUP-WI, and
the corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

193

Figure IV.6f. Correlation between land cover/use and shape on PLS4 for WUP-WI, and
the corresponding deformation grid. Two-letter acronyms are as in Table A.IV.1

194

xvii

CHAPTER I

INTRODUCTION

Overview
Even though most living species have been exposed to climatic shifts throughout their history,
the current scenario of change has been deemed unique inasmuch as temperature is rising
unusually quickly. Only during the last century it has been registered an increase of 0.6º C on
average temperatures, from which 0.2º C have been gained during the last 30 years (Stott et al.
2000, Walther et al. 2002). The environmental changes associated to this phenomenon are
expected to trigger responses in organisms, including migration, adaptation, speciation, and
extinction (Lynch and Lande 1993). Indeed, the effects of this accelerated climate change are
already noticeable in a vast number of organisms, affecting their physiology, phenology, and
distribution, which has led to local extinctions in some cases and adaptation to novel local
conditions in others, presumably affecting species interactions, and hence community structure
and composition (Hughes 2000, McCarty 2001, Walther et al. 2002, Parmesan and Yohe 2003,
Root et al. 2003, Parmesan 2006, Thomas et al. 2006, Salamin et al. 2010, and references
therein).
This study focuses on some of the consequences of rapid climatic change, specifically
changes in geographical distribution of land populations undergoing range expansion. Despite
the fact that changes in distribution do not necessarily imply evolutionary changes (Thomas et al.
2001), geographic expansion might create conditions for differentiation and divergence to occur.
Moreover, such expansions can have important effects on genetic diversity and spatial structure

1

of populations, leaving characteristic footprints on the genome (Stone et al. 2003, Zheng et al.
2003). Only in the last decades, we have witnessed the invasion of boreal forests in Northern
Michigan and Wisconsin by temperate species that were previously absent in those areas, often
accompanied by the declining of the ecologically similar boreal populations (Long, 1996; Myers
et al 2009). A clear example of this phenomenon is the effect that the expansion of Peromyscus
leucopus in the northern Great Lakes seems to have on populations of the dwarf deer mouse P.
maniculatus gracilis (Pmg), where the numbers of the last species are declining noticeably while
the former are flourishing (Myers et al. 2009). These two species, however, have coexisted in
northern Michigan despite their similarity in sympatry (Baker 1983, Myers et al. 2005) with no
evidence of hybridization (Dice 1933). Their coexistence has been attributed to their dissimilar
adaptation to winter, which favors survivorship in P. maniculatus, balanced by a higher
reproductive rate: in P. leucopus (Wolff 1996). Therefore, it is not surprising that the onset of
warmer climates gives P. leucopus a competitive advantage over Pmg.
One of the greatest concerns about rapid global warming is to understand how fast
species can respond to these changes and how entire communities may be affected. The
remarkable range expansion that Pln has experienced in northern Michigan and Wisconsin makes
it an ideal model to study the short term consequences of global warming and the dynamics of
range expansion. Although the effects of climate change are currently not very dramatic in this
area, it is predicted that they will have a major impact on the Great Lakes region, through drastic
temperature increase and changes in the water level of the lakes in the next decades (Dempsey
2004). These predicted changes will likely have major effects on the biological composition of
entire communities.

2

Because the processes associated to current climate change are expected to cause not only
changes in the genetic structure of populations but also on their morphological attributes (Root et
al. 2003), this study examines the genetic and morphological consequences of the geographical
expansion of Pln in the northern Great Lakes region, as well as the ecogeographical factors
related to this phenomenon. Similar studies found in the literature usually analyze either the
genetic (e.g., Excoffier et al. 2009) or morphological (e.g., Goheen et al. 2003) consequences of
range expansion, but few offer a comprehensive evaluation of the effects of range expansion
taking into account not only genetics and morphometrics, but also ecological factors. In this
study I have focused on the newly established populations of this mouse in Michigan’s Upper
Peninsula (UP) and their putative progenitor populations in northern Michigan’s Lower
Peninsula (LP) and Wisconsin.

Geographical expansion of Peromyscus leucopus in the northern Great Lakes region
The white footed mouse Peromyscus leucopus is one of the most common mammals in the
northeastern United States. It is a generalist rodent, inhabiting habitats as diverse as woodlands,
brushlands, and farmlands (Baker 1983). It is an unusually tractable species because of its
abundance across its range, and relative ease with which they to catch. Within its range, 17
subspecies have been proposed by different authors (Fig. I.1; Hall 1981, Lackey et al. 1985).
In the Great Lakes region, P. leucopus is represented by the subspecies P. l.
noveboracensis (Pln), which is widely distributed in deciduous forests throughout the central and
northeastern United States (Lackey et al. 1985). It is of small size, similar to a house mouse
(weighing less than 30 g); it has dark dorsal pelage, usually brown or caramel, and white venter;

3

its tail is shorter than its head-body length, and its dental formula is (Baker 1983). It is
omnivorous, incorporating seeds and nuts, fruits, grasses, and insects in its diet (Holman 2001).
In Michigan, Pln is at the northern limit of its range, and in this area it seems to prefer
woodlands and brushlands, although it is also found in open fields (Baker 1983). The population
abundance of this mouse can show dramatic seasonal fluctuations. Individuals normally have a
short life span of less than one year (Baker 1983). Studies of natural populations in southern
Michigan carried out in the 1940s suggest that the reproductive season spans early March to late
October, with females producing on average 4.3 young per litter and four litters per year (Burt
1940). Nests are made of mostly organic material and are built above ground, usually under
rocks or logs or in holes at the base of trees (Holman 2001).
In a study that combined 100 years of museum records with 14 years of intensive
collecting, Myers and co-workers (2005) showed the remarkable expansion of Pln in Michigan.
Historical records indicate that the northern limit of Pln in the early 1900s barely reached 45° N
(Osgood 1909), which corresponds to approximately three fourth of the LP in latitude; museum
records, however, suggest that this mouse might have been already present in the northern LP at
this time (Myers et al. 2009). In the late 1930s, an isolated population was reported from
Menominee County in the UP (Baker 1983). No further populations in the UP outside of
Menominee were found until the early 1990s and since then, Pln has expanded eastward about
225 Km across the UP. It has become the commonest small mammal at some localities in the UP
from which it was absent only 30 years ago (Fig. I.2; Myers et al. 2005, Myers et al. 2009). Its
range has extended into the central and eastern areas of the UP (Myers et al. 2009) and
northward in Wisconsin (Long 1996) and Minnesota (Jannett et al. 2007). The geographic

4

expansion of Pln is correlated with a trend toward milder winters, suggesting that climate change
may play a causal role in this process (Myers et al. 2005).
Range expansions of populations of Pln in the Great Lakes region, however, are not a
recent phenomenon; the margins of the range of this species have extended northwards since the
late Pleistocene, following the retreat of ice after the Last Glacial Maximum (LGM). During the
LGM, when the Laurentide Ice Sheet covered the Great Lakes region (18,000 YBP; Pielou
1991), many temperate species are believed to have survived in southern refugia that remained
ice-free (Pielou 1991, Rowe et al. 2004). As the ice receded, species like Pln expanded
northward, occupying the newly available habitat (Pielou 1991). In Michigan, the ice sheet began
to retreat about 14,800 YBP. Wisconsin and northern Michigan, including the UP, were ice-free
by about 10,000 YBP (Pielou 1991). By that time, southern Michigan was already covered by
diverse forests similar to those found today (Holman 2001). Presumably, Pln would have
inhabited those forests, and subsequently would have colonized the northern LP and WI as
conditions became increasingly more favorable. Although northward expansion after LGM is a
common process for most northern species who survived the last ice age, those expansions did
not take place at the pace they have been occurring over the last few decades.

Northern Great Lakes Ecoregions
The Great Lakes region is the final product of glacial erosion during the Pleistocene, although its
actual formation started in the Paleozoic era (Holman 2001). The Laurentide Ice Sheet covered
the Great Lakes regions during the Last Glacial Maximum, ca. 18,000 YBP; in Michigan, the ice
sheet began to retreat about 14,800 YBP. By 10,000 YBP, southern Michigan was covered by

5

diverse forests similar to those in the present, while northern Michigan was still too unstable to
sustain viable communities (Holman 2001).
The northern Great Lakes is the point of encounter of two ecological regions: the
northern forests and the eastern temperate forests (CEC 1997). The transition between these two
ecoregions represents the northern or southern distributional limit of many North American
species (e.g. Hall 1981). Current climate change seems to have diffused this barrier, and some
temperate species that had their northern limit at those margins are invading the northern forest
(Myers et al 2009). The northern forests are typified by extensive boreal forests of conifers, with
mixed hardwoods toward the south; the climate of this ecoregion is characterized by long, cold
winters and short, warm summers (CEC 1997). These forests are found in northern Wisconsin,
the UP, and some areas of the northern LP, Pennsylvania, and New York. The eastern temperate
forests occupy a vast area of the eastern United States and consist mostly of beech-maple and
maple-basswood forests, although mixed-oak associations are more frequent in the northern
Midwest; climate varies over a latitudinal gradient, from cool to subtropical; this area has a very
rich biodiversity (CEC 1997).

Objectives and Outline
The main objectives of this study are to investigate the phylogeographic and genetic structure,
and morphological differentiation of populations of the white-footed mouse in the northern Great
Lakes region, where the species is undergoing geographical expansion. Additionally,
ecogeographical factors are analyzed to investigate their association with genetic and
morphological changes observed in expanding populations of Pln. Among of the unknowns
addressed in this Dissertation are the origin of the newly founded populations in the UP, the

6

phenotypic changes that have accompanied the rapid geographical expansion of these mice, and
the geographical and environmental factors that account for morphological variation and the
differentiation between old and new populations. These questions are developed in three research
chapters:


Chapter II: this chapter focuses on the genetic structure and phylogeographic pattern of
Pln in northern Michigan and Wisconsin, which were inferred based on mtDNA sequence
polymorphism. The analysis of 595 sequences from 21 populations in Wisconsin, 9
populations in the UP, and 9 populations in the LP shows two well-differentiated
lineages. One of the lineages includes all populations from western UP and Wisconsin
(WUP-WI), and the other consists of all populations in the LP and eastern UP (LP-EUP).
A clear signature of spatial expansion is evident for both lineages, although demographic
expansion is only supported for WUP-WI. Contemporary expansion is evidenced by the
geographic proximity of divergent haplotypes. These results show that the UP is a zone
of secondary contact of these two very differentiated lineages, whose divergence
occurred over 30,000 YBP, which might have important consequences for the evolution
of the resulting admixed populations.



Chapter III: In this chapter, the morphological consequences of range expansion are
investigated based on the analysis of the shape of the mandible of Pln. The main
objective was to assess whether geographic expansion has been accompanied by
morphological divergence, and whether such divergence is attributable to random factors.
Shape analysis was based on geometric morphometric techniques. Results suggest that
new populations are significantly different from the old established ones, not only in
magnitude but also in the dimensionality of variation. Because the patterns of genetic and

7

morphometric distances are incongruent, morphological differentiation of new
populations is likely caused by non-random factors. Regardless of whether phenotypic
changes have been caused by plasticity and/or adaptation, it is clear that the rapid
geographic expansion that Pln has experienced in the northern Great Lakes has already
noticeable phenotypic consequences.


Chapter IV: This last chapter seeks to find the geographical and environmental factors
associated to patterns of genetic and morphometric variation observed in the expanding
populations of Pln in the northern Great lakes. Those factors are expected to explain not
only general variation in shape among populations, but also the divergence between old
and new populations in the two lineages. To this end, several statistical analyses were
carried out to study the association between genetic and morphometric variables with
several geographical and environmental variables factors, including ecological and
geographical distances, climate, geographical coordinates, and land cover/land use. For
LP-EUP, all factors analyzed show a similar effect on the shape of the mandible, i.e.,
contraction of mandibular processes. For WUP-WI, the different factors affect different
regions of the mandible. The examination of these geographical and environmental
factors shows that climate variables, in particular those associated to cold climate (e.g.,
low temperatures and snowfall) are the best predictors that explain the differentiation
between old and new populations in both lineages.

8

Figures

Figure I.1. Distribution of Peromyscus leucopus, according to Lackey et al. (1985). The
numbers indicate the distribution of the 17 sub-species proposed for this mouse in North
America. The grey area with the number 14 shows the range of the subspecies P. l.
noveboracensis. Note that in this map this subspecies is not reported for the Upper Peninsula of
Michigan. Reproduced with permission of the authors and the American Society of
Mammalogists, Allen Press, Inc.

9

Figure I.2. Distribution of the white-footed mouse, Peromyscus leucopus, in Michigan during a) 1883-1980 and b) 1981-2006
(Myers et al. 2009).

10

BIBLIOGRAPHY

11

BIBLIOGRAPHY

Baker RH (1983) Michigan Mammals. Michigan State University Press. East Lansing, Michigan.

Burt WH (1940) Territorial behavior and populations of some small mammals in Southern
Michigan. Miscellaneous Publications Museum of Zoology, University of Michigan, 45,
1-58.

Commission for Environmental Cooperation. 1997. Ecological Regions of North America:
Toward a common perspective. Biblioteque Nationale du Quebec, Montreal, Quebec.

Dempsey D (2004) On the brink. The Great Lakes in the 21st Century. Michigan State University
press, East Lansing, Michigan.

Dice LR (1933) Fertility relationships between some of the species and subspecies of mice in the
genus Peromyscus. Journal of Mammalogy, 14, 298-305.

Excoffier L, Foll M, Petit RJ (2009) Genetic Consequences of Range Expansions. Annual
Review of Ecology, Evolution and Systematics, 40, 481-501.

Goheen, JR, Swihart RK, Robins JH (2003) The anatomy of a range expansion: changes in
cranial morphology and rates of energy extraction for North American red squirrels from
different latitudes. Oikos, 102, 33-44.

Hall RE (1981) The Mammals of North America. Volume II. 2nd edition. Wiley-Interscience
Publication, New York.

Holman JA (2001) In the quest of Great Lakes Ice Age Vertebrates. Michigan State University
press, East Lansing, Michigan.

Hughes L (2000) Biological consequences of global warming: is the signal already. Trend in
Ecology & Evolution, 15, 56-61.

Lackey JA, Huckaby DG, Ormiston BG (1985) Peromyscus leucopus. Mammalian Species, 247,
1-10.
12

Long CA (1996) Ecological replacement of the deer mouse, Peromyscus maniculatus, by the
white-footed mouse, P. leucopus, in the Great Lakes Region. Canadian Field Naturalist,
110, 271-277.

Lynch M, Lande R (1993) Evolution and extinction in response to environmental change. In:
Biotic interaction and global change (eds. Kereiva. PM, Kingsolver JG, Huey RB), pp.
234-250. Sinauer Associates Inc, Sunderland, Massachusetts.

McCarty JP (2001) Ecological consequences of recent climate change. Conservation Biology,
15, 320-331.

Myers P, Lundrigan BL, Koppel RV (2005) Climate change and the distribution of Peromyscus
in Michigan—Is global warming already having an impact? . In: Mammalian
diversification: from chromosomes to phylogeography (eds. Lacey EA, Myers P), pp.
101-126. University of California Press, Berkeley, California.

Myers P, Lundrigan BL, Hoffman SMG, Haraminac AP, Seto SH (2009) Climate-induced
changes in the small mammal communities of the Northern Great Lakes Region. Global
Change Biology, 15, 1434-1454.

Osgood WH (1909) North American Fauna. No. 28. Revision of the mice of the American genus
Peromyscus. Government Printing Office, Washington.

Parmesan C, Yohe G (2003) A globally coherent fingerprint of climate change impacts across
natural systems. Nature, 421, 37-42.

Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annual
Review of Ecology Evolution and Systematics, 37, 637-669.

Root TL, Price JT, Hall KR, Schneider SH, Rosenzweig C, Pounds JA (2003) Fingerprints of
global warming on wild animals and plants. Nature, 421, 57-60.

Stone GN, Atkinson RJ, Brown G, Rokas A, Csóka G (2003) The population genetic
consequences of range expansion: oak gallwasp as a model system. In: Genes in the

13

environment (eds. Hails RS, Beringer JE, Godfray HCJ), pp. 46-62. Blackwell Publishing
and the British Ecological Society.

Stott PA, Tett SFB, Jones GS, Allen ME, Mitchell JFB, Jenkins GJ (2000) External control of
20th Century temperature by natural and anthropogenic forcings. Science, 290, 21332137.

Thomas CD, Bodsworth EJ, Wilson RJ, Simmons AD. Davies ZG, Musche M, Conradt L (2001)
Ecological and evolutionary processes at expanding range. Nature, 411, 577-581
Walther G-R, Post E, Convey P, Menzel A,. Parmesan C, Beebee TJ, Fromentin J-M, HoeghGuldberg O, Bairlein F (2002) Ecological responses to recent climate change. Nature,
416, 389-395.

Wolff JO (1996) Coexistence of white footed mice and deer mice may be mediated by
fluctuating environmental conditions. Oecologia, 108, 529-533.

Zheng X, Arbogast BS, Kenagy GJ (2003) Historical demography and genetic structure of sister
species: deermice (Peromyscus) in the North American temperate rain forest. Molecular
Ecology, 12, 711-724.

14

CHAPTER II

PHYLOGEOGRAPHY, GENETIC STRUCTURE AND HISTORICAL DEMOGRAPHY OF
EXPANDING POPULATIONS OF PEROMYSCUS LEUCOPUS IN THE NORTHERN
GREAT LAKES

Abstract
Rapid climate change affects the biology of many organisms, which is not surprising given the
key role that climate plays in determining the limits under which an organism can live. One
consequence of rapid climate change is shifts in geographical ranges. The white-footed mouse
(Peromyscus leucopus) has experienced a northern range expansion in Michigan and Wisconsin
over the last decades and is now common in localities where it was recently absent. Complete Dloop sequences were analyzed from 595 white-footed mice collected throughout the Great Lakes
region to investigate the origin of the newly established populations. It was found that mice from
Wisconsin and the western Upper Peninsula (UP) of Michigan make up a well-differentiated
lineage, while the eastern UP and the Lower Peninsula (LP) of Michigan form another. With few
exceptions, western populations in the UP originate from Wisconsin while the eastern UP
populations originate from LP. Both lineages display a genetic signature of spatial expansion,
while demographical expansion is apparent only for the western UP-Wisconsin lineage. These
results suggest that the UP is a zone of secondary contact for these two well differentiated
lineages, which requires further investigation of the potential consequences of admixture, e.g.,
outbreeding depression.

Introduction

15

Organisms naturally move across geographical regions; they may migrate seasonally or disperse
over their lifetimes. However, their distribution is limited by the suite of environmental
conditions to which they are adapted (Gaston 1991), which basically constitutes the species’
ecological niche (Sexton et al. 2009). Understand the factors that determine a species’ range and
how its boundaries are maintained over time are among the most challenging issues in
biogeography and evolutionary biology.
Climate is a key factor in determining the suitable living environment of organisms and
the shape of their distribution (Pigot et al. 2010). Substantial evidence suggests that the present
distribution of most organisms in North America was influenced significantly by the last
glaciation and its associated climate changes (Comes & Kadereit 1998, Taberlet et al. 1998,
Hewitt 2000, Zheng et al. 2003, Rowe et al. 2004, Rowe et al. 2006, Taylor & Keller 2007,
Centeno-Cuadros et al. 2009, Kerdelhue et al. 2009). Similarly, accelerated climatic change
documented in recent decades (Stott et al. 2000) is transforming environments and triggering
shifts in distribution. The effects of climate change on the geographic range of some species are
already noticeable, with populations sometimes declining as their habitats shrink (e.g., Beever et
al. 2003, Parmesan 2006, Stirling & Parkinson 2006, Moritz et al. 2008 Myers et al. 2009) , or in
other cases expanding their range (e.g., Parmesan et al. 1999, Thomas & Lennon 1999, Thomas
et al. 2001, Parmesan 2006, Jannett et al. 2007, Moritz et al. 2008, Myers et al. 2009).
Additionally, recent climate change is affecting the physiology and phenology of a wide variety
of organisms and has caused local extinctions and promoted adaptation to novel local conditions,
presumably affecting species interactions and hence community structure and composition (e.g.,
Brown et al. 1997, Visser et al. 1998, Brown et al. 1999, Hughes 2000, Bradshaw & Holzapfel

16

2001, McCarty 2001, Parmesan & Yohe 2003, Parmesan 2006, Root et al. 2003, Thomas et al.
2006, Salamin et al. 2010).
This study examines the genetic consequences of range expansion in the white-footed
mouse, Peromyscus leucopus noveboracensis (Pln), in the Great Lakes region. This rodent
species is one of the most common mammals in the northeastern United States. In northern
Michigan, it is at the northern limit of its range, and in that region it has experienced a very rapid
geographic expansion over the last century (Myers et al. 2009). Its range has extended into the
central and eastern parts of the Upper Peninsula of Michigan (Myers et al. 2009) and northward
in Wisconsin (WI) and Minnesota (Long 1996, Jannett et al. 2007). Myers and colleagues (2005)
hypothesized that this expansion is related to a rapid increase in local temperatures and an
observed trend toward shorter winters in the northern Great Lakes region. Indeed, according to
data from the National Oceanic and Atmospheric Administration (http://www.noaa.gov), the
mean winter temperature in the midwestern United States has increased by 4.2 °C (about 7.6 °F)
since 1970. In addition, the day of ice breakup for the Grand Traverse Bay, in the northwestern
LP, now takes place almost 50 days earlier than it did at the beginning of the 20th century,
(Myers et al. 2005). This shows a remarkable advance in the 2000s, in comparison with the 11.4
days reported for this same area by Magnuson and colleagues (2000) between 1850 and 1995. In
both cases, the data show clear evidence of earlier onset of spring conditions. If the northern
limit of Pln’s range is influenced by its tolerance to winter (Wolff 1996, Myers et al. 2005),
shorter winters and earlier springs may explain the northward expansion of this species.
Little information about the phylogeography and post-glacial expansion of Pln is
available, despite it being one of the most studied rodent species in North America. Only one
study (Rowe et al. 2006) has investigated the phylogeographic patterns of this mouse, . These

17

authors used mtDNA sequence data to test the hypothesis that Pln and another species, the
eastern chipmunk, both associated with deciduous forests, share a history of expansion into
regions previously covered by the Laurentide Ice Sheet. They found support for three mutually
monophyletic clades in both species, with phylogenetic topologies suggesting a biogeographic
pattern that could not be explained by the contemporary landscape. These observations led the
authors to suggest that modern populations arose from different glacial refugia. Moreover, the
phylogeographic pattern found for these two species is congruent with that of a subset of
deciduous trees, suggesting a shared pattern of association and expansion since the Last Glacial
Maximum (LGM).
In this chapter, gene genealogies and population genetic methods based on coalescent
theory are used to investigate the phylogeography, genetic structure, and historical demography
of newly founded populations of Pln in northern WI and the UP, as well as their putative
progenitor populations in the LP and WI, with the general objective of understanding the genetic
consequences of historical and recent range expansion. Both demographic explosions and spatial
expansion leave characteristic footprints on the genome (Cann et al. 1987, Rogers & Harpending
1992, Hewitt 2000, Lessa et al. 2003, Zheng et al. 2003, Rowe et al. 2004, Rowe et al. 2006,
Taylor & Keller 2007, Centeno-Cuadros et al. 2009, Kerdelhue et al. 2009). For any population,
the genetic consequences of contemporary range expansions are strongly affected by the
evolutionary history previous to the expansion, because past population processes determine the
phylogenetic and genetic diversity of the colonizers (Taylor & Keller 2007). Additionally,
contemporary events affecting the number and diversity of sources of emigrants contributing to
the newly established populations, e.g., whether dispersers come from a single or multiple
populations (Slatkin 1977), will also impact their genetic make-up. Consequently, knowledge of

18

both the phylogeographic and genetic pattern of post-Pleistocene expansion, as well as the life
history and population dynamics of Pln are crucial for understanding the genetic consequences of
current range expansion.
We evaluate the following three hypotheses regarding the genetic diversity,
differentiation, and connectivity of Pln populations in the LP, UP and WI: (1) mouse populations
came to occupy these regions after the retreat of the Laurentide ice sheet and the concomitant
expansion of deciduous forests; therefore, a genomic signature of range expansion should be
evident in all populations, characterized by high haplotype diversity, low nucleotide diversity,
and unimodal mismatch distributions; (2) because WI and the UP are geographically continuous,
mouse populations in those regions are likely exchanging genes and should have relatively high
genetic similarity; and (3) LP lineages will be equally differentiated from WI and UP
populations, because the Great Lakes should be an effective geographical barrier.

Materials and Methods
Sampling
Five-hundred and eighty-four tissue samples were obtained from mice collected in 39 counties in
Michigan and Wisconsin, encompassing the current range of Pln including newly invaded
localities (Fig. II.1, Table II.1). Multiple trapping sites within 25 km of each other and in the
same general environment were considered to represent a single population, and the geographical
midpoint of these combined sites was determined (http://www.geomidpoint.com) and used to
calculate geographic distances among populations. In the text that follows, each population is
given the name of the county of origin, and LP (140 samples from 9 counties), UP (166 samples
from 9 counties), and WI (278 samples from 21 counties) are referred to as ―regions.‖

19

Collecting sites were in forested areas above 44° North in both Michigan and WI and
were chosen based on previous knowledge of the local mouse populations, type of vegetation,
and accessibility. Roughly 200,000 km2 was covered by this sample design. Adjacent
populations were separated by 50 km on average. Additionally, Pln tissue samples from Butler
Co., Ohio (2); Johnson Co., Illinois (5); and Livingston Co., Michigan (4) were used to represent
potential southern refuge sources. Three samples of Peromyscus maniculatus bardii (Pmb) from
the LP (1), WI (1), and OH (1), and 14 P. m. gracilis (Pmg) from the LP (4), UP (5), and WI (5)
were used as outgroups in these analyses.
The majority of samples (389) were from mice that were live-trapped between 2003 and
2007, using both small (5 x 6.4 x 16.5 cm) and large (7.6 x 8.9 x 22.9 cm) Sherman folding traps.
Captured mice were weighed to the nearest 0.1 g and the following information was recorded:
length of right ear and hind foot (in mm), age class (juvenile, sub-adult, adult), and reproductive
status (descended testicles for males; pregnancy and visible nipples for females). A piece of
tissue was clipped from the right ear and placed in a 1.6 ml centrifuge tube with SET buffer
(Brinster et al. 1985). Tissues were stored in a cooler while in the field and later moved to a -20°
C laboratory freezer. All field-captured mice were released at the site of capture. Additional
samples (195) were obtained from tissue collections of the Michigan State University Museum,
the University of Michigan Museum of Zoology, and the Museum of Natural History of the
University of Wisconsin at Stevens Point, or from frozen carcasses provided by the Wisconsin
Department of Natural Resources. Animals were handled according Michigan State University
Institutional Animal Care and Use Committee protocol # 08/04-108-00, in compliance with the
American Society of Mammalogists guidelines for the use of wild mammals in research (Gannon
et al. 2007).

20

Laboratory Procedures
DNA was extracted from mouse ear tissue using Qiagen DNeasy Tissue Kits (Qiagen Inc.),
following the manufacturer’s protocol. The complete D-loop and some nucleotides of the
adjacent tRNAs were amplified using universal primers, L15926 5’
TACACTGGTCTTGTAAACC 3’ and H00651 5’ AAGGCTAGGACCAAACCT 3’, located on
the Pro and Phe tRNAs, respectively (Kocher et al. 1989). Polymerase chain reactions (PCR)
were carried out in a volume of 25 μl, containing 1.75 mM of MgCl, 0.2 mM of each dNTP, 1
μM of each primer, 1 U of Platinum Taq (Invitrogen Co.), and 50-100 ng of DNA template, on
an Eppendorf Mastercycle thermocycler (Eppendorf Co.). After 2 min of denaturation at 95° C,
35 cycles were run under the following conditions: denaturation at 95° C for 45 s, annealing at
56° C for 60 s, and expansion at 72° C for 75 s, were run. DNA products were checked on a
1.5% agarose gel (Invitrogen Co.) and cleaned using QIAquick Gel Extraction Kits (Qiagen Inc.)
in preparation for direct automated sequencing. DNA was sequenced at the University of
Michigan DNA Sequencing Core, in an Applied Biosystems DNA Sequencer, Model 3730 XL
(Applied Biosystems Inc.). The same primers as used in PCR reactions were used for
sequencing; however, reverse sequencing was incomplete because the reaction terminated at
approximate 200 bp. This phenomenon was also reported by Rowe and collaborators (2006),
while sequencing the control region of Peromyscus leucopus and P. maniculatus. According to
these authors, sequencing from the Phe tRNA stopped after 180 bp due to a palindromic
sequence (GGGGGGAAGGGGGG) interfering with the reaction. Sequences obtained were
between 1000 and1050 bp in length. A representative of each haplotype has been submitted to
GenBank (accession numbers GU810187- GU810356).

21

Gene Genealogy
Sequences were visually checked and edited using Chromas Pro Version 1.49 (Technelysium Pty
Ltd.), aligned using ClustalX (Thompson et al. 1997), and reduced to haplotypes using Collapse
2.1 (http://darwin.uvigo.es). A gene genealogy, based on haplotypes only, was estimated through
a Neighbor Joining analysis (implemented by Mega 4,Tamura et al. 2007), with gaps and
missing data completely excluded from the analysis. Pairwise distances were estimated by the
Maximum Composite Likelihood method (Tamura et al. 2004). The model of nucleotide
substitution used by Mega 4 with this method is that of Tamura & Nei (1993). Support for each
node was evaluated with 1000 bootstraps repetitions (Felsenstein 1985).

Gene Diversity and Population Structure
Genetic diversity was measured at both nucleotide and haplotype levels. The most common
measure of nucleotide diversity is ∏ (Nei 1987), which is calculated as ∑

where xi

and xj are the frequencies of the ith and jth DNA sequences in the sample and πij is the
proportion of nucleotide differences between them. For haplotypes, H is the most appropriate
measure of genetic diversity. This parameter is defined in terms of haplotype frequencies;
equivalent to the expected heterozygosity for diploid data, it measures the probability that two
randomly chosen haplotypes in the sample are different. It is estimated as
(

∑

)

where n is the number of gene copies, k is the number of haplotypes, and pi is the frequency of

22

the ith haplotype in the sample (Nei 1987). Both ∏ and H were estimated using Arlequin 3.11
(Excoffier et al. 2005).
Spatial population structure and potential barriers to gene flow were investigated using
Spatial Analysis of Molecular Variance (SAMOVA) v.1 (Dupanloup et al. 2002). This program
performs a simulated annealing procedure that assigns populations to K groups defined a priori.
The composition of these groups takes into account the geographic proximity and genetic
homogeneity of populations, so that the proportion of total variance due to differences between
groups (Fct) is maximized. Spatial population structure was evaluated for K = 2, 3, and 4 groups,
with 100 starting conditions; statistical significance was assessed with 1000 permutations.
Population genetic structure among groups obtained with SAMOVA was inferred by
analysis of molecular variance (AMOVA, Excoffier et al. 1992). AMOVA was performed using
Arlequin 3.11 (Excoffier et al. 2005). The significance of the estimations was assessed using
1000 permutations. Because the power of statistical inferences depends on sample size, only
samples with 8 individuals or more were included in the analyses unless otherwise noted.

Migration and Divergence Time
Migration rates and the timing of population splits were estimated for the two groups obtained
from SAMOVA, using the isolation with migration model, as implemented by the IMa program
(Hey & Nielsen 2007). The full model has six parameters: the population mutation rates of the
descendent (1, 2) and ancestral populations (a); the number of migrants/1000 generations
entering population 1 from population 2 (m1), and the number of migrants/1000 generations
entering population 2 from population 1 (m2); and time since divergence (t). The program
estimates the posterior probability distributions for each parameter through MCMC simulations

23

of gene genealogies. Simulations were run assuming the HKY mutation model (Hasegawa et al.
1985). Several preliminary simulations were run to find appropriate priors for each parameter.
Once priors that produced bell-shaped curves with low dispersion of points were found for each
parameter, 1 x 106 cycles for the burn-in time and 10 x 106 cycles for the Markov chain were
run. The burn-in time is the number of initial cycles that were discarded; burn-in provides a way
to ―dememorize‖ the process and thus avoid dependence on initial conditions. To assess
consistency of the parameters estimated across runs, each simulation was repeated 3 times, each
time using a different random seed and on a different machine. The simulation that yielded the
highest effective sample size was chosen for further analysis. Additionally, nested models were
examined and compared to the full model using a likelihood ratio test (2LLR) to identify the
simplest model that fit the data.
The estimate of model parameters is given by the mode of their posterior probability
distributions, and those estimates can be converted to demographical parameters. The effective
population size of females, Nef, can be obtained from the equation Nef = /ug, where u is the
mutation rate per year per locus and g is the number of years per generation; the migration rate
between populations per generation is equal to m = mu; and the divergence time, in generations,
is t = t/u (Nielsen & Wakeley 2001, Hey & Nielsen 2004, 2007). The mutation rate was
calculated based on the genetic divergence observed between Pln and Pmg in mtDNA D-loop
sequences, assuming a molecular clock. These two species appeared in the fossil record about
500,000 years before present (Hibbard 1968). Haplotypes for 43 Pln and 37 Pm were compared;
of the 913 aligned sites (after excluding gaps and missing data), there were 82 substitutions
(~13%), which yields a mutation rate, μ, of 1.588*10-7 mutations/year/site/lineage. This estimate
is roughly the same as the one obtained by Rowe et al. (2006) and Taylor & Hoffman (2010) for
24

the same mitochondrial region, but using shorter sequences. However, any estimation based on
this rate have to be interpreted with caution because it is uncertain the actual date when Pl and
Pm diverged. The mutation rate for the whole D-loop was obtained by multiplying μ by the total
length of the sequence (934 bp), which yielded a rate of 1.2202*10-4
mutations/year/locus/lineage. Obtaining a precise estimate of the generation length for this
species is very challenging, because reproductive rates vary depending on when females are born
during the season, which is also influenced by geography, and because of overlapping
generations. In this study, 2 generations/year was assumed (Miller 1989, Rowe et al. 2006,
Centeno-Cuadros et al. 2009, Taylor & Hoffman 2010).

Isolation by Distance
To investigate whether the pattern of genetic structure observed in Pln might be explained by
limited dispersal and isolation by distance, a Mantel test was performed to assess the correlation
between genetic and geographic distances. Pairwise geographical distances between the
geographical midpoints of each pair of populations were calculated as the shortest distance over
the earth’s surface. For the matrix of genetic distances, pairwise st were obtained with Arlequin
3.11 (Excoffier et al. 2005). For the correlation, linearized st (i.e. st /(1-st) and the log of
genetic distance were used (Rousset 1997). The significance of the correlation was evaluated
using 1000 permutations. All calculations for this test, and obtaining values for the geographic
distance matrix, were done with Matlab R2006a (The MathWorks, Inc.). This analysis was
carried out for the whole sample and for each group obtained with SAMOVA.

Demographical and Spatial Expansions

25

Sudden (demographical) expansion (Rogers & Harpending 1992) was modeled using a method
implemented by Arlequin 3.01 (Excoffier et al. 2005), based on the distribution of the frequency
of mistmatches among all sequences. The demographical parameters estimated are τ, the time to
expansion in generations, and θ0 and θ1, the neutral population mutation rates before and after
the expansion. Given that τ = 2ut, where u is the neutral mutation rate per locus per year (Sherry
et al. 1994), θ0 = 2μN0 and θ1 = 2μN1, the time to the expansion (in generations) and change in
population size were calculated. To convert t to YBP, t was multiplied by the number of years in
one generation (i.e., 0.5). Confidence intervals for these parameters were calculated using a
parametric bootstrap that compares the observed mismatch distribution to a simulated
distribution under the sudden expansion model. The sum of square differences (SSD) between
the observed and the simulated mismatches was used to calculate a P-value, which served as the
test statistic. This P-value takes on large values if the observed distribution fits a model of
sudden expansion, and small values if the data correspond to a stationary population. The
Harpending’s raggedness index, r, was also calculated; this index takes on large values in
stationary populations and small values in populations that have undergone demographical
expansion.
Spatial expansion was examined using a continent-island model implemented in Arlequin
3.01 (Excoffier et al. 2005), where the island exchanges m migrants with a continent that has a
population of infinite size. Three parameters were estimated: τ, θ (= θ0 = θ1), and M (= 2Nm),
from the distribution of the frequency of mistmatches among all sequences. The statistical test
and the fit of the observed mismatch to a model of spatial expansion were performed following
the same procedure as for the model of sudden expansion described above. The mismatch

26

distributions for the whole sample, and the groups identified by SAMOVA, were obtained and
compared to models of sudden and spatial expansion.
In addition, Fu’s Fs values (Fu 1997) for each group were calculated using DnaSp 5.1
(Librado & Rozas 2009). Fs takes on negative values when there are many new mutations (i.e.
rare alleles), which is usually the case in expanding populations; large negative values of Fs
indicate departure from neutrality and therefore is an indicative of demographic expansion.

Results
Sequencing and Gene Genealogy
Aligning, editing, and fitting all sequences to the same length resulted in a 934 bp segment that
included the complete D-loop and a few nucleotides from each of the flanking tRNAs. The 595
total sequences included 170 unique haplotypes. The majority of haplotypes were found in 1-3
copies, but a handful were present in more than 10 and the maximum number of copies observed
for a haplotype was 38. The LP and WI regions shared no haplotypes, the UP and WI shared 16,
and the LP and UP shared 4. The 4 haplotypes shared between LP and UP were apparently of LP
origin: Mackinac (UP) shared one haplotype with Alpena, Cheboygan, Presque Isle, and
Crawford counties (all in the LP); Chippewa (UP) shared one haplotype with Cheboygan Co. and
one with Alpena Co. (both in the LP); and two haplotypes from Ontonagon Co. (UP) were also
found in Otsego Co. (LP).
A Neighbor Joining tree based on 187 sequences (i.e., the total number of observed Pln
haplotypes) plus 17 Pm (outgroup) sequences (GenBank accession numbers GU810357GU810373), revealed two well-supported (>90% bootstrap) clades: Pln and Pm (Fig. II.2). The
Pln clade has two main sub-branches, one clustering the haplotypes from the Western UP (WUP)

27

and WI (72% bootstrap) and the other grouping the haplotypes from the LP and Eastern UP
(EUP, 40% bootstrap). With the exception of those from two counties in the eastern UP
(Chippewa and Mackinac), mice from the UP clustered with those from WI. Chippewa and
Mackinac mice were genetically more similar to LP than WI mice, with only 2 haplotypes from
Chippewa and 1 from Mackinac found in the WUP-WI clade (―+‖ in Fig. II.2). While most other
haplotypes were in the expected clade based on geography, there were a few exceptions (―*‖ in
Fig. II.2): one haplotype from the LP (Benzie Co.) clustered with the WUP-WI haplotypes, and
one haplotype from WI (Jackson Co.) and another from the UP (Ontonagon Co.) were found in
the LP-EUP clade. All of the haplotypes from Illinois and Ohio were contained in the LP-EUP
clade. Although a significant genetic distance separates WUP-WI and LP-EUP clades, the
internal clades are shallow, with very short branches suggesting few nucleotide differences.

Gene Diversity and Differentiation
Haplotype diversity across populations was relatively high, between 0.69 and 0.96, while
nucleotide diversity was rather low, from 0.04 to 0.18 (Table II.1). Marathon Co., in central WI,
showed the highest haplotype diversity, while Mackinac Co., in Michigan’s eastern UP, had the
lowest. The highest and lowest nucleotide diversities were found in the eastern UP (Chippewa
and Delta Counties, respectively). Regional differences for both parameters were small, with the
highest haplotype diversity observed in WI and the lowest haplotype diversity in Michigan’s UP.
The SAMOVA separates the Pln populations into 2 groups (Table II.2), in agreement
with the results of the gene genealogy. This analysis suggests that the populations from WUP
and WI make up one group and LP and EUP populations define a second. Increasing the value of
K creates groups each comprising only one population, while the variance due to inter-group

28

differences (Фct) does not change, thus not adding information. Differences among all
populations (Фst) explain over 70% of the genetic variance observed, while differences between
the two groups defined by SAMOVA account for over 65% of this variance; differentiation
among populations within groups is moderate (Фsc = 14%).

Migration and Divergence Time
The three simulations that were run after appropriate priors had been identified produced very
similar results. For each parameter of the model, the posterior probability distributions were
smooth and bell-shaped. The results suggest that the population that gave rise to those currently
in the LP-EUP and WUP-WI split in the late Pleistocene, about 34,000 YBP (Table II.3). The
descendent populations probably experienced demographical expansion, since their female
effective population sizes are larger than the one estimated for the ancestral population. Current
migration between the two descendent populations is near zero and historical migration between
them has been asymmetrical. The comparison of nested models to the full-parameter model was
run twice. In one of the runs, only one model was accepted (Table II.4), the model in which the
female effective population sizes of the ancestral population and LP-EUP were the same.
However, the value of 2LLR was very close to the X2 critical value for 1 degree of freedom and
 = 0.05. In an independent run, none of the nested models was accepted, suggesting that the
full-parameter model best explains the data.

Isolation by Distance
The Mantel test indicated a high and significant positive correlation between genetic and
geographic distance among all populations (r = 0.41, P < 0.001). Because this test can be

29

sensitive to extreme values and does not take into account geographic barriers like the Great
Lakes, the test was performed separately for each partition obtained with SAMOVA (Fig. II.3, a
and b). The positive correlation between genetic and geographical distances remained significant
in each of the groups (LP-EUP: r = 0.58, P < 0.01; WUP-WI: r = 0.33, P < 0.01), even after
adjusting for multiple comparisons using a Bonferroni’s correction. To better visualize the
relationship between genetic and geographic distances, I mapped populations and connected
them with color-coded lines based on their levels of genetic similarity, i.e., colder colors for
more similar genetically and warmer colors for more differentiated. Two maps were generated,
one connecting populations within lineages with up to 35% of genetic differences, and another
connecting populations across lineages with differences over 35% (Fig. II.4). At those levels of
genetic differentiation, all WI and western UP populations are connected, and the eastern UP is
connected to the entire LP; however, UP populations are the most different in both lineages. The
greatest genetic distances are between populations from the LP-EUP and those from the WUPWI, with the exception of one locality in the eastern UP (i.e. Chippewa), which stands out as the
least different from populations in WUP-WI group.

Demographical and Spatial Expansions
The analysis of pairwise differences (mismatch distributions) showed that for the two groups of
Pln defined by SAMOVA, and for all populations pooled together, the data fit a model of spatial
expansion (Table II.5, Figs. II.5 a, c, and e). However, the distribution of LP-EUP is ragged (Fig.
II.5a), and the distribution of all populations considered at the same time is bi-modal (Fig. II.5e).
It appears that LP-EUP started to expand first, as early as 26,000 YBP, followed later by WUPWI populations, about 18,000 YBP (Table II.5). Grouping all populations suggests that the

30

spatial expansion process might have begun around 60,000 YBP. Estimated levels of gene flow
are high within each group, with an average 13 migrants per generation in LP-EUP, almost 35 in
WUP-WI, and over 50 when all populations are considered.
The comparison of the mismatch distributions to Roger and Harpending’s model of
sudden expansion (Rogers & Harpending 1992) is significant for the WUP-WI, and when all
populations are considered together, but not for the LP-EUP (Table II.5), indicating that the latter
group has not experienced recent demographical growth. The mismatch distribution of LP-EUP
is ragged (Fig. II.5b), and the distribution of all data is bi-modal (Fig. II.5f). This result is not
congruent with Fu’s Fs, which is significant for the two groups and all populations, indicating a
departure from neutrality, usually caused by changes in population size (Table II.5). The
estimated times since demographical expansion for each group are similar to those for spatial
expansion: LP-EUP about 29,000 YBP, WUP-WI roughly 19,000 YBP, and all populations
approximately 45,000 YBP.

Discussion
This study used variation in mtDNA control region sequences to investigate the phylogeography,
genetic structure, and historical demography of populations of Pln in the northern Great Lakes
region. The results showed that populations of Pln sort into two well-differentiated groups, one
comprising populations from the WUP-WI, and the other populations from the LP-EUP. This
grouping was based on genetic similarity and was concordant with geographic proximity. The
two groups appeared to be reciprocally monophyletic, and in both, a signature of spatial
expansion was identified. However, demographical expansion is apparent only for WUP-WI.
The two groups diverged from each other about 34,000 YBP, and shortly afterward (26,000

31

YBP), the LP-EUP population started to expand, followed by the WUP-WI population about
6,000 years later.

Gene genealogy
The tree topology has the characteristics typical of populations that have gone through recent
demographic/spatial expansion, with most haplotypes represented either by a single or few
copies, and only a small number of haplotypes widespread. This produces a shallow tree, which
is more evident for the WUP-WI clade than for LP-EUP one. With the exception of a few
haplotypes, these two clades appear clearly differentiated. However, because of some unresolved
nodes in each clade, the bootstrap support for LP-EUP is only 40% and for WUP-WI, 72%. A
few haplotypes were not found in the cluster expected given their collection locality; this might
be explained by ancestral polymorphism, natural dispersal (across ice or by rafting), and/or
human-mediated dispersal. Given that WI and LP have no haplotypes in common and are
reciprocally monophyletic, ancestral polymorphism is an unlikely explanation for most of these
exceptional haplotypes. Natural dispersal, although possible, is also unlikely because it would
require that mice cross Lake Michigan. Instead, human-mediated dispersal is the most plausible
explanation. The haplotype from Benzie (LP) found in the WUP-WI clade, and the haplotype
from Jackson (WI) found in the LP-EUP clade, might reflect transport of mice across the lake by
ferry. Two vehicle and passenger ferries cross Lake Michigan, one between Muskegon (MI) and
Milwaukee (WI), and another between Ludington (MI) and Manitowoc (WI). These piers are not
necessarily close to the localities were mice used in this study were collected, but they provide a
potential mechanism for these mice to cross the lake. These mice are frequently found occupying
human dwellings and even nesting in vehicles; they undoubtedly are carried by people on

32

occasion. Similarly, the mice of likely LP origin found in UP counties (Chippewa, Mackinac,
and Ontonagon) might have crossed the Straits of Mackinac with humans via the Mackinac
Bridge. The other two haplotypes found in Chippewa were most similar to haplotypes from
Illinois and not found in either the LP or WI. Because it is unlikely that the similarity between
those haplotypes would be caused by gene flow between Illinois and the UP, ancestral
polymorphism would be a better explanation.
The gene genealogy found in this study is consistent with the findings of Rowe et al.
(2006), who analyzed 575 partial sequences from the mtDNA D-loop of Peromyscus leucopus
from 50 sites in the eastern United States, comprising almost the entire distribution of the
species. Those authors found 3 clades: a western clade, grouping haplotypes from Wisconsin,
Minnesota, Illinois, Indiana, Iowa, and Missouri; a central clade, consisting of haplotypes from
southern Illinois, southern Indiana, Kentucky, Ohio, Tennessee, and Pennsylvania; and an
eastern clade, with haplotypes from Michigan, southern Illinois, Indiana, Ohio, Tennessee, West
Virginia, North Carolina, Pennsylvania, New York, Rhode Island, Massachusetts, New Jersey,
Nova Scotia, Maine, and New Hampshire. Although their sampling did not include extensive
representation of northern WI and northern Michigan, the WUP-WI lineage found in this study
fell in their western clade, and the LP-EUP clade in their eastern clade.

Gene Diversity and Differentiation
Further evidence of a recent demographic/spatial expansion is provided by gene diversity
indices. The high haplotype diversity and low nucleotide diversity are the result of a high
frequency of derived haplotypes differentiated by only a few mutations. Spatial and
demographical expansion are usually accompanied by an initial reduction in haplotype diversity,

33

and the consequent loss of information about deep coalescent events, due to a reduction in
population size associated with bottlenecks or founder effects (Excoffier et al. 2009).
The SAMOVA yielded two interesting results: (1) Lake Michigan is (not surprisingly) an
effective barrier to gene flow between the LP and WI; and (2) the newly founded mouse
populations in the UP have two different origins, in the LP and in WI. Genetic differences
between mice from these two regions explain more than 65% of the total genetic variation
observed. The partition supported by this analysis is consistent with the LP-WUP and WUP-WI
clades found in the gene genealogy.

Isolation by Distance
White-footed mice have a home range of about 1,000 m2 (Burt 1940) and may disperse up to
1km in an unfragmented landscape (Krohne et al. 1984); although Maier (2002) reported
exceptional long-distance movements of two female in New England, which were recaptured at
almost 15 and 8 km from the original capture site. Given these short dispersion distances, it is
not surprising that populations of these mice display a pattern of isolation by distance. Moreover,
it has been shown that dispersal in this species is male biased (Jacquot and Vessey 1995), and
because genetic distances were estimated from polymorphism in mtDNA, which is maternally
inherited, this pattern reflects female philopatry. Although significant for both regions, the
pattern of isolation by distance was more evident for LP-EUP. The largest genetic distance was
between LP and WI (Fig. II.4), supporting the idea that Lake Michigan has been an effective
barrier for preventing gene flow.

Spatial and Demographical Expansion

34

Depending on the levels of historical and current gene flow among subpopulations, both
demographical and spatial expansions can have similar genetic consequences, characterized by a
smooth, unimodal mismatch distribution (Harpending 1994, Excoffier 2004, Excoffier et al.
2009). In light of this, goodness of fit of models of sudden (pure demographical) expansion and
stepwise (spatial) expansion were fitted to the data. As expected, both regions showed a
significant signature of spatial expansion. Independently of where populations survived during
the glacial maximum, they moved northward after the ice sheet started to retreat. In both regions,
gene flow among populations is high; the estimated number of migrants/generation is higher for
WUP-WI than for LP-EUP, probably reflecting differences in the landscape permeability of
these two regions. Habitat fragmentation does not limit the dispersal capability of Pln (Mossman
and Waser 2001), but urban areas do, as reported by Munshi-South and Kharchenko (2010) in
New York City, where populations of this mouse were fragmented and genetically differentiated.
Range expansion has occurred repeatedly in the history of terrestrial organisms
(Excoffier et al. 2009) and is perhaps the most remarkable feature of the Quaternary Period
(Hewitt 2000). I found evidence of Pln range expansion from different sources. As mentioned
above, the shallow depth of the gene genealogy and the characteristics of the gene diversity
indices are typical of populations that have undergone recent range expansion (Hewitt 2000,
Taylor & Keller 2007). Moreover, the analysis of pairwise differences suggests a major role of
post-glacial range expansions in generating the pattern of polymorphism observed in sequences
of the mtDNA control region of mice in these populations. This genetic signature of post-glacial
expansion is a common trait of many organisms in the northern hemisphere (Lessa et al. 2003).
Its effect was so pervasive here that when all samples were analyzed together, it was not possible
to distinguish the effects of contemporary expansions, presumably because populations have not

35

yet recovered from the loss of genetic structure caused by post-Pleistocene range expansion
(Pinho et al. 2007, Taylor & Keller 2007, Kerdelhue et al. 2009). Qualitative evidence of
contemporary range expansion is only apparent when the pattern of isolation by distance is
examined. Several populations in WUP-WI (Fig. II.4), share genetically similar haplotypes
while they are geographically distant from one another, which suggests movement of those
haplotypes as populations expanded.
During spatial expansion, population size is usually reduced due to founder effects; if the
population is successfully established, population size should increase and a demographical
expansion footprint should be evident in the genome. Contrary to this expectation, a model of
demographical expansion did not fit the LP-EUP data, which displayed a ragged mismatch
distribution (Fig. 5a). Additionally, the only nested model tested with IMa that fit the data
significantly (in only one of the two simulations carried out) was the model in which a
(ancestral population mutation rate) and 1 (LP-EUP population mutation rate) were equal
(Table II.4), also suggesting a stationary population. A plausible explanation for this incongruent
pattern could be that the signature of geographical expansion has been obscured by admixture.
The pattern of mismatch distribution found here shows a number of pairwise differences > 20
and in low frequency (Fig. II.5b), indicative of an old expansion (Rogers & Harpending 1992).
Rowe et al. (2004) found a similar pattern for the eastern chipmunk, Tamias striatus, and argued
that this is the signature of two independent expansion episodes. Evidence that admixture can
erase (or overwhelm) the signature of demographical expansion is found in the North American
populations of the invasive plant Silene latifolia (Taylor & Keller 2007), which does not show a
pattern of demographical expansion despite the fact that its ancestors did experience post-glacial
expansion. This weed is original from Europe and was introduced into North America around

36

200 YBP. While European populations show a clear footprint of post-glacial demographical
expansion, North America populations display very low genetic diversity and evidence of
admixture caused by introductions from two differentiated lineages, a history that is reflected in
the bimodal mismatch distribution of haplotypes from this region.
Similarly, LP-EUP mice should have experienced a demographical expansion, as
suggested by the significantly negative Fu’s Fs, but the genetic signature of that expansion might
have been erased by at least two waves of expansion in this region. The agreement of the
estimated time ranges of expansion for the study populations (26,000 to 14,000 YBP for LPEUP, and 18,000 to 9,000 YBP for WUP-WI) with the actual date that the Laurentide ice sheet
began to retreat after the LGM is striking. Southern Michigan was ice free about 14,800 YPB,
while WI and northern Michigan were ice free around 10,000 YBP (Pielou 1991). Moreover,
these expansion times fall within the range obtained by Rowe et al. (2006) for the expansion of
Pln and Tamias striatus, and by Taylor & Hoffman (2010) for Pmg, from this same region.
Although these results should be interpreted with caution, because they are affected by the
estimated mutation rate and the assumed generation time, such consistency among three different
species strongly and coherence with known geological events suggest a common demographical
history for small terrestrial mammals in the region.

Differentiation between LP and WI, and Origin of UP Populations
These data suggest that the LP and WI populations began to diverge several thousand years
before the formation of Lake Michigan, about 45,000 YBP (Table II.3). At that time, they may
have occupied a common refugium or several different refugia south of the glaciated region
(Rowe et al. 2006); subsequently, as they expanded to newly open northern lands, populations

37

would have become isolated by cooling cycles. Later, the formation of Lake Michigan would
have reinforced this isolation. The estimation of demographical parameters from a divergence
with migration model implemented by IMa suggests that gene flow between the LP-EUP and
WUP-WI is extremely low (Table II.3). Mice from these two regions are now genetically
distinct, and thus the hypothesis that Lake Michigan has been an effective barrier to gene flow
between these two regions is well supported. However, this Lake has been a permeable barrier at
its northern end, where there is evidence of gene flow between populations of Pmg in the UP,
northern LP and Beaver islands (an island group in Lake Michigan; Taylor & Hoffman 2010).
Taylor & Hoffman (2010) found a similar date (59,000 YBP) for the divergence between eastern
and western clades of Pmg in this same region, suggesting a common history of expansion.
I hypothesized that recently established populations in the UP originated in WI, but
results showed that there have been two expansion pathways: one west to east, from WI to the
WUP, and another south to north from the LP to the EUP. A similar pair of pathways was found
for UP populations of Pmg, except that for that (boreal) species, the source of EUP populations is
to the north in Canada (Taylor & Hoffman 2010). In addition, Pmg is an old resident of the UP
that has inhabited this region for thousands of years, while Pln is a recent invader that has been
found over most of this area only over the last two decades (Burt 1957; Myers et al. 2009). This
is an intriguing pattern of faunal sorting, and I wonder if other species show a similar pattern of
distribution. Michigan’s Upper Peninsula is known to have a very complex climate, with
conditions changing dramatically from shore to inland and from the mainland to the peninsula.
The WUP has a strongly continental climate, with very cold winters and vegetation characterized
by mostly northern hardwoods (Albert 1995). The EUP climate is greatly influenced by the Great
Lakes, with warmer winters and cooler summers relative to other areas at the same latitude in the

38

WUP and WI, and with vegetation characterized by northern hardwoods, pine forests, hardwoodconifer and conifer swamps (Anderson 1986). Such differences in climate and vegetation surely
should play a key role as selective agents on the distribution of organisms in the UP.

Evolutionary Implications and Future Directions
This study suggests different expansion histories for Pln in the LP-EUP vs. WUP-WI. It is likely
that the populations that ultimately gave rise to these two clades not only survived in different
refugia during the last episode of glaciation, but also that the WUP-WI populations originated
from a single source (as suggested by the unimodal mismatch distribution), while the LP-EUP
populations were likely founded by mice from more than one source (as shown by the ragged
mismatch distribution). Because the sampling here did not include populations from the southern
limit of the range of Pln, and because the expansion process has likely resulted in the loss of
deep coalescence events, it is not clear where the LP populations originated. The WI population,
however, may have arisen from animals that survived the last glaciation in a more northern
refugium, as suggested by Rowe et al. (2004, 2006) for T. striatus and Pln currently found in WI
and northern Illinois. Alternatively, if WI populations originated in some southern refugium, the
observed pattern could have been caused by the loss of all ancestral haplotypes, either because of
a selective sweep on the mtDNA or as a consequence of the expansion process. The fact that two
different studies show a similar genetic pattern for two different species (T. striatus and Pln) in
the same region, suggests a common demographical history that affected their genomes in
comparable ways. This alternative hypothesis can be assessed by broadly sampling Pln
populations from unglaciated locations to identify possible refugial origins. Additionally, the

39

examination of other molecular markers would allow a test of whether a selective sweep might
account for the lack of deep coalescent events in our samples.
Spatial expansion of Pln has had an impact on the genetic diversity of northern
populations of this mouse. Michigan’s Upper Peninsula is a potential zone of secondary contact
for two differentiated lineages of Pln. Further research, with nuclear DNA, needs to be carried
out to evaluate the extent of gene flow between eastern and western UP populations, and the
possible phenotypic effects of admixture. Additionally, the expansion of this mouse has affected
the composition of northern Great Lakes communities, where it is outcompeting Pmg (Myers et
al. 2005). As discussed by Myers et al. (2009), although the consequences of this change are
unpredictable, it can be anticipated that it will have a major impact given the ecological role of
small mammals in these communities. Also of great concern is the fact that Peromyscus
leucopus is a natural reservoir of Borrelia burgdoferi, the causative agent of Lyme disease in this
region, and the expansion of Pln may influence the spread of that disease.

40

Tables

Table II.1. Origin of the samples analyzed in this study, where n is the sample size, h is the
number of haplotypes observed, H is the haplotype diversity, S is the number of segregating sites
over the complete sequence, Π is the nucleotide diversity, and Lat and Long are the decimal
geographic coordinates of the geographic midpoint of each locality. NA indicates that these
parameters were not calculated for those populations because of small sample size.
Michigan’s Lower Peninsula (LP)
County

Code

n

h

Alpena

AP

17

7

Benzie

BZ

15

9

Charlevoix

CX

17

7

Cheboygan

CH

19

7

Crawford

CR

20

8

Emmet

EM

9

4

Grand
Traverse

GD

8

4

Otsego

OT

17

12

Presque Isle

PI

18

7

140

65

Total LP

H

S

0.8676
(+/- 0.0503)
0.9333
(+/- 0.0397)
0.8603
(+/- 0.0475)
0.8480
(+/- 0.0515)
0.8684
(+/- 0.0475)
0.7500
(+/- 0.1121)
0.7500
(+/- 0.1391)
0.9485
(+/- 0.0368)
0.8562
(+/- 0.0466)
0.9856

24

41

22

35

27

19

19

29

21
56

41

Π
0.009069
(+/- 0.004946)
0.013877
(+/- 0.007434)
0.008063
(+/- 0.004438)
0.010896
(+/- 0.005825)
0.008566
(+/- 0.004645)
0.007540
(+/- 0.004442)
0.007988
(+/- 0.004765)
0.009213
(+/- 0.005018)
0.008857
(+/- 0.004820)
0.010087

Lat
(north)

Long
(west)

44.9997 83.5471

44.6067 85.1494

45.2614 84.8037

45.5634 84.6637

44.7888 84.7340

45.5431 85.0478

44.6448 85.4962

45.0549 84.5688

45.4372 84.0489

Table II.1. Continued
(+/- 0.0023)

( +/- 0.005162)

Michigan’s Upper Peninsula (UP)
County

n

h

Alger

AL

12

8

Chippewa

CW

17

8

Delta

DE

20

7

Gogebic

GO

17

11

Mackinac

MW

20

7

Menominee

MN

20

6

Marquette

MA

20

7

Ontonagon

ON

22

10

Schoolcraft

SC

18

6

166

71

n

h

20

13

Total UP

H

S

0.9242
(+/- 0.0575)
0.8529
(+/- 0.0663)
0.7684
(+/- 0.0689)
0.9412
(+/- 0.0364)
0.6895
(+/- 0.1053)
0.7579
(+/- 0.0645)
0.8105
(+/- 0.0689)
0.8918
(+/- 0.0419)
0.7582
(+/- 0.0701)
0.9797
(+/- 0.0032)

18

48

15

24

37

13

18

41

11

75

Π
0.005281
(+/- 0.003118)
0.018187
(+/- 0.009534)
0.004057
(+/- 0.002386)
0.006490
(+/- 0.003644)
0.010401
(+/- 0.005560)
0.004875
(+/- 0.002798)
0.005543
(+/- 0.003134)
0.010757
(+/- 0.005709)
0.004353
(+/- 0.002552)

Lat
(north)

Long
(west)

46.3412 86.7884

46.4320 84.7266

46.0788 86.8345

46.3437 89.3678

46.0025 84.8929

45.1580 87.7000

46.8736 87.8958

46.7364 89.7394

46.0886 86.2278

0.013992
(+/- 0.007011)

Wisconsin (WI)
County
Adams

AD

H

S

0.9421
(+/- 0.0340)
42

32

Π
0.007680
(+/- 0.004203)

Lat
(north)

Long
(west)

43.9012 89.8676

Table II.1. Continued
Ashland

AS

10

8

Barron

BN

18

11

Bayfield

BY

20

13

Buffalo

BF

20

12

Burnett

BU

15

9

Douglas

DG

11

6

Forest

FO

1

NA

Iron

IR

13

7

Jackson

JK

20

11

Marathon

MH

20

14

Marinette

MR

18

9

Oneida

ON

11

6

Outagamie

OU

16

9

Ozaukee

OZ

10

6

0.9556
(+/- 0.0594)
0.9412
(+/- 0.0328)
0.9421
(+/- 0.0340)
0.9421
(+/- 0.0318)
0.8762
(+/- 0.0697)
0.8364
(+/- 0.0887)
NA

22

27

26

26

29

19
NA

0.8718
(+/- 0.0670)
0.8684 (+/0.0640)
0.9632
(+/- 0.0255)
0.8889
(+/- 0.0488)
0.8000
(+/- 0.1138)
0.8583 (+/0.0772)
0.8444
(+/- 0.1029)

43

19

44

29

19

16

26

12

0.007552
(+/- 0.004388)
0.007005
(+/- 0.003890)
0.006675
(+/- 0.003700)
0.007313
(+/- 0.004020)
0.007001
(+/- 0.003937)
0.005567
(+/- 0.003297)
NA
0.004738
(+/- 0.002814)
0.008678
(+/- 0.004702)
0.007333
(+/- 0.004030)
0.006571
(+/- 0.003672)
0.005830
(+/- 0.003436)
0.006456
(+/- 0.003643)
0.004589
(+/- 0.002812)

46.1568 90.4755

45.6839 92.0936

46.3705 91.3376

44.5892 91.8022

46.0418 92.1906

46.2722 91.6982
45.9640 88.9674
46.3579 90.3774

44.2055 90.4484

44.8175 89.7509

45.5900 87.9464

45.8709 89.6527

44.3961 88.6663

43.4039 87.9906

Table II.1. Continued
0.8039

PO

18

9

Price

PR

4

NA

NA

NA

NA

45.6751 90.1751

Rusk

RU

3

NA

NA

NA

NA

45.3686 90.5151

Shawano

SH

2

NA

NA

NA

NA

44.7173 88.5953

Taylor

TA

12

8

Wood

WO

16

8

Total WI

278/
268*

159

(+/- 0.0907)

0.8939
(+/- 0.0777)
0.8833
(+/- 0.0522)
0.9941
(+/- 0.0010)

29

0.007666

Portage

29

18

110

(+/- 0.004222)

0.007995
(+/- 0.004531)
0.005679
(+/- 0.003247)

44.5412 89.5654

45.3127 90.6268

44.3785 90.1090

0.007639
(+/- 0.003982)

*The number above indicates the total number of samples from the region; the number below,
which omits small samples, was used for the calculation of h, H, S, and Π.

44

Table II.2. Spatial population structure of Pln in the Northern Great Lakes, as suggested by
spatial analysis of molecular variation (SAMOVA). This analysis was run assuming K = 2, 3 and
4 groups.
Source of Variation

Fixation Index† Percentage of Variation§

2 groups
Among all populations

st = 0.70653

70.66

Among groups

c t= 0.65897

65.90

Among populations within groups

sc = 0.13945

4.76

Within populations

29.35

3 groups
Among all populations

st = 0.70271

70.27

Among groups

ct = 0.66452

66.45

Among populations within groups

sc = 0.11383

3.82

Within populations

29.73

4 groups
Among all populations

st = 0.69927

69.93

Among groups

ct = 0.66127

66.13

Among populations within groups

sc = 0.11219

3.80

Within populations

30.07

†All p-values were < 0.001
§Note that the variance explained by differences among all populations is equal to the variance
explained by differences among groups plus the variance explained by differences among
populations within groups.

45

Table II.3. Demographic parameters estimated using the isolation with migration model implemented in IMa. 1, 2, and a are the
mutation rates of population 1 (LP-EUP), population 2 (WUP-WI), and the ancestral population, respectively; N1, N2, and N3 are the
effective population sizes of females (number of individuals) of LP-EUP, WUP-WI, and the ancestral population, respectively; m1 and
m2 are the migration rates per generation into population 1 from population 2, and into population 2 from population 1, respectively;
and t is the time since the two groups split, in years. Conversions of model parameters to demographic parameters were done
assuming 0.5 years/generation and a mutation rate of 1.2202*10-4 substitutions/ locus/year. HiPt is the value with the highest
frequency (i.e. mode), and HPD90Lo and HPD90Hi are the lower and upper bounds of the estimated 90% Highest Posterior Density
Interval.

Value

Population mutation rate

Females effective population size

migration rates/generation

Time

1

2

a

N1

N2

Na

m1

m2

t

HiPt

87.1341

237.2285

39.508

1428194

3888354

647566

4.54525*10-06

6.7111*10-07

34469.76

HPD90Lo

66.5683

199.7049 20.3854

1091105

3273314

334132.1

1.7998*10-06

1.28121*10-07 28274.05

282.329

1859905

4627586

1309916

9.60908*10-06 1.92792*10-06 45090.97

HPD90Hi 113.4728

79.918

46

Table II.4. Comparison of nested models to the full-parameter model of isolation with
migration, using IMa. 1, 2, and a are the mutation rates for population 1 (LP-EUP),
population 2 (WUP-WI), and the ancestral population, respectively; m1 and m2 are the migration
rates per generation into population 1 from population 2, and into population 2 from population
1, respectively.
Model

Log (P)

2LLR

d.f.

1 2 a m1 m2

-5.2036

‒

‒

1 2 a m1 = m2

-7.9641

5.521

1

1 2 a m1 m2 = 0

-9.2249

8.0425

1*

1 2 a m1 = 0 m2

-29.5687

48.7303

1*

1 2 a m1 = m2 = 0

-38.1349

65.8627

2*

1 = 2 a m1 m2

-16.9191

23.4311

1

1 = 2 = a m1 m2

-23.4555

36.5038

2

1 = 2 a m1 = m2

-18.0838

25.7604

2

1 = 2 a m1 = m2 = 0

-56.8541

103.301

3*

1 = 2 = a m1 = m2

-24.9274

39.4476

3

1 = 2 = a m1 = m2 = 0

-64.6224

118.8376

4*

1= a 2 m1 m2

-7.0983

3.7894 §

1

1= a 2 m1 = m2

-9.9595

9.5118

2

1= a 2 m1 = m2 = 0

-40.3996

70.3921

3*

1 2 = a m1 m2

-11.722

13.0369

1

1 2 = a m1 = m2

-14.0896

17.772

2

1 2 = a m1 = m2 = 0

-45.6183

80.8294

3*

§The only model that was not rejected; marginally significant (X21, 0.5 = 3.841). In an
independent run, this model was not accepted (2LLR = 3.9149)
*When one or more parameters are set to zero, the expected distribution of 2LLR is a mixture. In
the case when only one parameter is fixed, 2LLR will distribute as a random variable that would
take both the value of zero and the value of X21 with probability 0.05 (see Hey and Nielsen
2007)

47

Table II.5. Estimated parameters of the Spatial and Sudden Expansion models, and results of tests of goodness of fit of observed data
to each model for each group suggested by SAMOVA, and for the whole sample. M is the number of migrants per generation, τ is the
time to expansion, in generations; 0 and 1 are the population mutation rates before and after expansion, respectively; r is Harpending’s
raggedness index, and P is the goodness of fit to each model.
Spatial Expansion
Regions

LPEUP
WUPWI

All

N

M
(95%
CI)

τ
(95%
CI)

Estimated
expansion
time (YBP)

177

12.77
(6.4527.9)

10.31
(6.8712.93)

14,08026,500

397

34.44
(21.7182.69)

6.27
(4.478.91)

9,16318,262

574

52.97
(0.49137.31)

3.734
(1.8831.43)

3,85264,396

Sudden Expansion
P

0
(95%
CI)

1
(95%
CI)

τ
(95%
CI)

Estimated
expansion
time (YBP)

r

P

Fu’s
Fs

0.10*

0
(02.96)

37.02
(23.43419.76)

10.94
(6.7314.15)

13,79729,000

0.02*

0.01

-4.40*

0.77*

0.446
(02.41)

57.93
(32.551184.81)

7.09
(4.659.57)

9,53319,614

0.51

0.43*

97.32*

0.002 0.49*

13.869
(031.80)

173.57
(31.102824.82)

3.84
(1.8622.20)

3,80945,493

0.002

0.24*

76.69*

r

0.02

0.01

*p-values  0.01

48

Figures

Figure II.1. Geographic location of Peromyscus leucopus populations analyzed in this study.
The code for each population is as in Table II.1. LI = Livingstone, Michigan; IL = Johnson,
Illinois; OH = Butler, Ohio. Black circles indicate populations in the former northern limit of the
species in the region. Grey circles indicate populations in newly invaded localities.

49

Figure II.2. Neighbor Joining tree showing the haplotype genealogy of Peromyscus leucopus in Wisconsin (WI), and Michigan’s
Lower Peninsula (LP) and Upper Peninsula (UP). Peromyscus maniculatus (Pm) is the outgroup. Some internal clades have been
magnified for better visualization. ―+‖ indicate haplotypes from Chippewa County (CW) that cluster with the WUP-WI clade and ―*‖
haplotypes that are not in the expected clade given collection locality. Bootstrap support values are given for the main nodes.
Acronyms as in Table II.1; LI = Livingstone, Michigan; IL = Johnson, Illinois; OH = Butler, Ohio.

50

a
Figure II.3. Results from Mantel’s tests evaluating the correlation between genetic (Φst/1-Φst) and geographical (km) distances
between populations. a) Population comparisons within the LP-EUP group. c) Population comparisons within the WUP-WI.
Correlation coefficients and P values are indicated in the text.

51

b

Figure II.4. Interconnectivity among populations based on pairwise genetic distances. The color of the line indicates genetic distance
between population pairs, based on pairwise Φst/1-Φst values, with dark blue indicating the most similar genetically and dark red the
most different. On the left are represented populations with genetic differences up to 35%, and on the right populations with genetic
differences > 35%. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

52

Spatial Expansion

Sudden Expansion

b

a
0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0

0
1

6

11

16

21

26

31

36

1

6

11

c

21

26

31

36

d

0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0

Frequency

16

0
1

6

11

16

21

26

1

31

6

11

e

16

21

26

31

f

0.1

0.1

0.05

0.05

0

0
1

6

11

16

21

26

31

36

1

6

11

16

21

26

31

Number of pairwise differences
Figure II.5. Distribution of pairwise differences of the observed data (dark grey) and the
simulated data (light grey) under a model of stepwise (spatial) expansion (left) and sudden
(demographic) expansion (right) for a and b (LP-EUP), c and d (WUP-WI), and e and f (all
data). Results of tests of goodness of fit of the data to each model are shown in Table II.5.

53

36

BIBLIOGRAPHY

54

BIBLIOGRAPHY

Albert DA (1995) Regional landscape ecosystems of Michigan, Minnesota, and Wisconsin: a
working map and classification. Gen. Tech. Rep. NC-178. St. Paul, MN: U.S.
Department of Agriculture, Forest Service, North Central Forest Experiment Station.
Northern Prairie Wildlife Research Center Home Page.
http://www.npwrc.usgs.gov/resource/1998/rlandscp/rlandscp.htm (Version 03JUN98).

Beever EA, Brussard PE, Berger J (2003) Patterns of apparent extirpation among isolated
populations of pikas (Ochotona princeps) in the Great Basin. Journal of Mammalogy, 84,
37-54.

Bradshaw WE, Holzapfel CM (2001) Genetic shift in photoperiodic response correlated with
global warming. Proceedings of the National Academy of Sciences of the United States of
America, 98, 14509-14511.
Brinster RL, Chen HY, Trumbauer ME, Yagle MK, Palmiter RD (1985) Factors affecting the
efficiency of introducing foreign DNA into mice by microinjecting eggs. Proceedings of
the National Academy of Sciences of the United States of America, 82, 4438-4442.
Brown JH, Valone TJ, Curtin CG (1997) Reorganization of an arid ecosystem in response to
recent climate change. Proceedings of the National Academy of Sciences of the United
States of America, 94, 9729-9733.

Brown JL, Li SH, Bhagabati N (1999) Long-term trend toward earlier breeding in an American
bird: A response to global warming? Proceedings of the National Academy of Sciences of
the United States of America, 96, 5565-5569.
Burt WH (1940) Territorial behavior and populations of some small mammals in southern
Michigan. University of Michigan Press, Ann Arbor, Michigan.

Burt WH (1957) Mammals of the Great Lakes. University of Michigan Press, Ann Arbor,
Michigan.

Cann RL, Stoneking M, Wilson AC (1987) Mitochondrial DNA and human evolution. Nature,
325, 31-36.

Centeno-Cuadros A, Delibes M, Godoy JA (2009) Phylogeography of southern water vole
(Arvicola sapidus): evidence for refugia within the Iberian glacial refugium? Molecular
Ecology, 18, 3652-3667.
55

Comes HP, Kadereit JW (1998) The effect of quaternary climatic changes on plant distribution
and evolution. Trends in Plant Science, 3, 432-438.

Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing approach to define the
genetic structure of populations. Molecular Ecology, 11, 2571-2581.

Excoffier L (2004) Patterns of DNA sequence diversity and genetic structure after a range
expansion: lessons from the infinite-island model. Molecular Ecology, 13, 853-864.

Excoffier L, Foll M, Petit RJ (2009) Genetic Consequences of Range Expansions. Annual
Review of Ecology, Evolution and Systematics, 40, 481-501.

Excoffier L, Laval G, Schneider S (2005) Arlequin (version 3.0): An integrated software package
for population genetics data analysis. Evolutionary Bioinformatics, 1, 47-50.

Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric
distances among dna haplotypes - application to human mitochondrial-DNA restriction
data. Genetics, 131, 479-491.

Felsenstein J (1985) Confidence limits on phylogenies - an approach using the bootstrap.
Evolution, 39, 783-791.

Gannon WL, Sikes RS, Comm ACU (2007) Guidelines of the American Society of
Mammalogists for the use of wild mammals in research. Journal of Mammalogy, 88, 809823.

Gaston KJ (1991) How large is a species geographic range. Oikos, 61, 434-438.

Harpending H (1994) Gene frequencies, DNA sequences, and human origins. Perspectives in
Biology and Medicine, 37, 384-394.

Hasegawa M, Kishino H, Yano TA (1985) Dating of the human ape splitting by a molecular
clock of mitochondrial-DNA. Journal of Molecular Evolution, 22, 160-174.

56

Hewitt G (2000) The genetic legacy of the Quaternary ice ages. Nature, 405, 907-913.
Hey J, Nielsen R (2004) Multilocus methods for estimating population sizes, migration rates and
divergence time, with applications to the divergence of Drosophila pseudoobscura and
D. persimilis. Genetics, 167, 747-760.

Hey J, Nielsen R (2007) Integration within the Felsenstein equation for improved Markov Chain
Monte Carlo methods in population genetics. Proceedings of the National Academy of
Sciences of the United States of America, 104, 2785-2790.

Hibbard CW (1968) Paleontology. In: Biology of Peromyscus (Rodentia) (ed. King JA), pp. 6-24.
Special Publication No. 2 of the American Society of Mammalogists.

Hughes L (2000) Biological consequences of global warming: is the signal already apparent?
Trends in Ecology & Evolution, 15, 56-61.

Jacquot JJ, Vessey SH (1995) Influence of the natal environment on dispersal of white-footed
mice. Behavioral Ecology and Sociobiology, 37, 407-412.

Jannett FJ, Broschart MR, Grim LH, Schaberl JP (2007) Northerly range extensions of
mammalian species in Minnesota. American Midland Naturalist, 158, 168-176.

Kerdelhue C, Zane L, Simonato M, et al. (2009) Quaternary history and contemporary patterns
in a currently expanding species. BMC Evolutionary Biology, 9, 220.

Kocher TD, Thomas WK, Meyer A, et al. (1989) Dynamics of mitochondrial DNA evolution in
animals - Amplification and sequencing with conserved primers. Proceedings of the
National Academy of Sciences of the United States of America, 86, 6196-6200.

Krohne DT, Dubbs BA, Baccus R (1984) An analysis of dispersal in an unmanipulated
population of Peromyscus leucopus. American Midland Naturalist, 112, 146-156.

Lackey JA, Huckaby DG, Ormiston BG (1985) Peromyscus leucopus. Mammalian Species, 247,
1-10.

57

Lessa EP, Cook JA, Patton JL (2003) Genetic footprints of demographical expansion in North
America, but not Amazonia, during the Late Quaternary. Proceedings of the National
Academy of Sciences of the United States of America, 100, 10331-10334.

Librado P, Rozas J (2009) DnaSP v5: a software for comprehensive analysis of DNA
polymorphism data. Bioinformatics, 25, 1451-1452.

Long CA (1996) Ecological replacement of the Deer Mouse, Peromyscus maniculatus, by the
White-footed Mouse, P. leucopus, in the Great Lakes region. Canadian Field Naturalist,
110, 271-277.

Magnuson JJ, Robertson DM, Benson BJ, et al. (2000) Historical trends in lake and river ice
cover in the Northern Hemisphere. Science, 289, 1743-1746.

Maier TJ (2002) Long-distance movements by female White-footed Mice, Peromyscus leucopus,
in extensive mixed-wood forest. Canadian Field-Naturalist, 116, 108-111.

McCarty JP (2001) Ecological consequences of recent climate change. Conservation Biology, 15,
320-331.

Moritz C, Patton JL, Conroy CJ, et al. (2008) Impact of a century of climate change on smallmammal communities in Yosemite National Park, USA. Science, 322, 261-264.

Mossman CA, Waser PM (2001) Effects of habitat fragmentation on population genetic structure
in the white-footed mouse (Peromyscus leucopus). Canadian Journal of Zoology, 79,
285-295.

Munshi-South J, Kharchenko K (2010) Rapid, pervasive genetic differentiation of urban whitefooted mouse (Peromyscus leucopus) populations in New York City. Molecular Ecology,
19, 4242-4254.
Myers P, Lundrigan BL, Hoffman SMG, Haraminac AP, Seto SH (2009) Climate-induced
changes in the small mammal communities of the Northern Great Lakes Region. Global
Change Biology, 15, 1434-1454.

Myers P, Lundrigan BL, Koppel RV (2005) Climate change and the distribution of Peromyscus
in Michigan—Is global warming already having an impact? . In: Mammalian

58

diversification: from chromosomes to phylogeography (eds. Lacey EA, Myers P), pp.
101-126. University of California Press, Berkeley, California.

Myers P, Lundrigan BL, Koppel RV (2005) Climate change and the distribution of Peromyscus
in Michigan—Is global warming already having an impact? . In: Mammalian
diversification: from chromosomes to phylogeography eds. Lacey EA, Myers P), pp.
101-126. University of California Press, Berkeley, California.

Nei M (1987) Molecular evolutionary genetics. Columbia University Press, New York.
Nielsen R, Wakeley J (2001) Distinguishing migration from isolation: A Markov Chain Monte
Carlo approach. Genetics, 158, 885-896.

Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annual
Review of Ecology Evolution and Systematics, 37, 637-669.

Parmesan C, Ryrholm N, Stefanescu C, et al. (1999) Poleward shifts in geographical ranges of
butterfly species associated with regional warming. Nature, 399, 579-583.

Parmesan C, Yohe G (2003) A globally coherent fingerprint of climate change impacts across
natural systems. Nature, 421, 37-42.

Pielou EC (1991) After the Ice Age: the return of life to glaciated North America. University of
Chicago Press, Chicago, Illinois.

Pigot AL, Owens IPF, Orme CDL (2010) The environmental limits to geographic range
expansion in birds. Ecology Letters, 13, 705-715.

Pinho C, Harris DJ, Ferrand N (2007) Contrasting patterns of population subdivision and
historical demography in three western Mediterranean lizard species inferred from
mitochondrial DNA variation. Molecular Ecology, 16, 1191-1205.

Rogers AR, Harpending H (1992) Population growth makes waves in the distribution of pairwise
genetic differences. Molecular Biology and Evolution, 9, 552-569.

Root TL, Price JT, Hall KR, et al. (2003) Fingerprints of global warming on wild animals and
plants. Nature, 421, 57-60.
59

Rosenzweig C, Karoly D, Vicarelli M, et al. (2008) Attributing physical and biological impacts
to anthropogenic climate change. Nature, 453, 353-357.

Rousset F (1997) Genetic differentiation and estimation of gene flow from F-statistics under
isolation by distance. Genetics, 145, 1219-1228.

Rowe KC, Heske EJ, Brown PW, Paige KN (2004) Surviving the ice: Northern refugia and
postglacial colonization. Proceedings of the National Academy of Sciences of the United
States of America, 101, 10355-10359.

Rowe KC, Heske EJ, Paige KN (2006) Comparative phylogeography of eastern chipmunks and
white-footed mice in relation to the individualistic nature of species. Molecular Ecology,
15, 4003-4020.

Salamin N, Wuest RO, Lavergne S, Thuiller W, Pearman PB (2010) Assessing rapid evolution in
a changing environment. Trends in Ecology & Evolution 25, 692-698.

Sexton JP, McIntyre PJ, Angert AL, Rice KJ (2009) Evolution and Ecology of Species Range
Limits. Annual Review of Ecology Evolution and Systematics 40, 415-436.
Slatkin M (1977) Gene flow and genetic drift in a species subject to frequent local extinctions.
Theoretical Population Biology, 12, 253-262.

Stirling I, Parkinson CL (2006) Possible effects of climate warming on selected populations of
polar bears (Ursus maritimus) in the Canadian Arctic. Arctic, 59, 261-275.

Stott PA, Tett SFB, Jones GS, et al. (2000) External control of 20th century temperature by
natural and anthropogenic forcings. Science, 290, 2133-2137.

Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF (1998) Comparative phylogeography and
postglacial colonization routes in Europe. Molecular Ecology, 7, 453-464.

Tamura K, Dudley J, Nei M, Kumar S (2007) MEGA4: Molecular evolutionary genetics analysis
(MEGA) software version 4.0. Molecular Biology and Evolution, 24, 1596-1599.

60

Tamura K, Nei M (1993) Estimation of the number of nucleotide substitutions in the control
region of mitochondrial DNA in humans and chimpanzees. Molecular Biology and
Evolution, 10, 512-526.

Tamura K, Nei M, Kumar S (2004) Prospects for inferring very large phylogenies by using the
neighbor-joining method. Proceedings of the National Academy of Sciences of the United
States of America, 101, 11030-11035.

Taylor DR, Keller SR (2007) Historical range expansion determines the phylogenetic diversity
introduced during contemporary species invasion. Evolution, 61, 334-345.

Taylor ZS, Hoffman SMG (2010) MtDNA genetic structure transcends natural boundaries in
Great Lakespopulations of deer mice (Peromyscus maniculatus gracilis). Canadian
Journal of Zoology, 88, 404-415.

Thomas CD, Bodsworth EJ, Wilson RJ, et al. (2001) Ecological and evolutionary processes at
expanding range margins. Nature, 411, 577-581.

Thomas CD, Franco AMA, Hill JK (2006) Range retractions and extinction in the face of climate
warming. Trends in Ecology & Evolution, 21, 415-416.

Thomas CD, Lennon JJ (1999) Birds extend their ranges northwards. Nature, 399, 213-213.

Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG (1997) The CLUSTAL X
windows interface: flexible strategies for multiple sequence alignment aided by quality
analysis tools. Nucleic Acids Research, 25, 4876-4882.

Visser ME, van Noordwijk AJ, Tinbergen JM, Lessells CM (1998) Warmer springs lead to
mistimed reproduction in great tits (Parus major). Proceedings of the Royal Society of
London Series B-Biological Sciences, 265, 1867-1870.

Wolff JO (1996) Coexistence of white footed mice and deer mice may be mediated by
fluctuating environmental conditions. Oecologia, 108, 529-533.

61

Zheng XG, Arbogast BS, Kenagy GJ (2003) Historical demography and genetic structure of
sister species: deermice (Peromyscus) in the North American temperate rain forest.
Molecular Ecology, 12, 711-724.

62

CHAPTER II

MORPHOMETRIC DIFFERENTIATION OF EXPANDING POPULATIONS OF
PEROMYSCUS LEUCOPUS IN THE NORTHERN GREAT LAKES: HISTORICAL AND
GENETIC FACTORS
Abstract
The white-footed mouse, Peromyscus leucopus, has experienced a remarkable northward range
expansion in the Great Lakes region over the last decades. In this study, I investigate whether the
rapid geographic expansion of populations of this species in the northern Great Lakes region has
been accompanied by rapid morphological divergence, and whether the pattern of divergence
suggests a neutral process or reflects the influence of non-random factors. Shape was measured
and analyzed using geometric morphometric techniques for statistical inference and graphical
interpretation. Significant differences were found between recently expanded and long
established populations not only in the magnitude of shape variation but also in the
dimensionality of variation, which suggest that the expansion process has been accompanied by
rapid phenotypic change. Additionally, the pattern of shape variation is incongruent with the
phylogeographic pattern observed in these populations, ruling out that changes in morphology
are caused by random phenotypic drift. These phenotypic changes might have been promoted by
phenotypic plasticity and/or adaptation, and provide clear evidence that although the geographic
expansion of Peromyscus leucopus in the northern Great Lakes is relatively recent, the
phenotypic consequences are already noticeable.

Introduction

63

Overview
A rapid geographical expansion of the white-footed mouse, Peromyscus leucopus, in the
northern Great Lakes region has been taking place over the last 30 years. This expansion is
particularly remarkable in Michigan’s Upper Peninsula (UP), where this species was almost
completely absent just a few decades ago (with the exception of populations found in the
southernmost county, Menominee), and is now common throughout most of the peninsula
(Myers et al. 2009). Phylogeographic analyses of white-footed mouse populations in the northern
Great Lakes, including northern Wisconsin (WI), the UP, and Michigan’s northern Lower
Peninsula (LP), show a clear genetic footprint of expansion and reveal an admixture of newly
founded populations in the UP with western populations originating in WI and eastern
populations from the LP (see Chapter II). The short time frame of this demographic expansion
raises questions about the evolutionary processes that might have allowed the rapid invasion of
novel environments at higher latitudes. These processes might be reflected in geographic patterns
of phenotypic variation in ecologically relevant phenotypic traits, such as foraging structures
(e.g., the mammalian mandible).
Rapid morphological divergence is not uncommon among populations of small
mammals. For example, Pergams and colleagues (Pergams & Ashley 1999; Pergams & Lacy
2008; Pergams & Lawler 2009) found rapid morphological changes (over circa 100 years) in
several species of rodents inhabiting both mainland and island environments. Morphological
changes in island mammals, especially changes in size, are usually attributed to island effect, e.g.
limited resources and reduced diversity, predation and competition (Millien 2006). In contrast,
the rapid morphological changes observed in mainland mammals are more often associated with
the impact of human activities and climate alterations.

64

Evidence of rapid morphological change has been found in expanding populations of
North American red squirrels in the central United States (Goheen et al. 2003). These squirrels
are experiencing a geographic expansion similar to that of Peromyscus leucopus. Since the late
1800s, they have expanded their range from central hardwoods to conifer-dominated forests.
Differences in diet between animals inhabiting these habitats are reflected in cranial
morphologies, which have diverged significantly in accord with changes in mandibular force.
Another example of rapid morphological change is found in house mouse populations on
the Mediterranean islands of Corsica and Sardinia and the small islet of Piana, near France
(Renaud & Auffray 2010). Populations on the two larger islands are thought to have arrived with
the first human settlers, while the population on the islet was first documented about 200 years
ago. The shape of the mandible of animals from the mainland is different from that on the
islands, and mandible shape on the small islet differs from that on both the mainland and the two
larger islands. The authors argue that because multiple colonizations of the larger islands would
counteract the effect of drift, morphological differences between those mice and mainland mice
are likely related to ecological causes, in particularly widening of the dietary niche of the island
mice. In contrast, differentiation of mice on the islet might reflect founder effect and drift, a
consequence of reduced population size and lack of gene flow with the adjacent island
populations.
As the preeminent anatomical element of masticatory function, and given its correlation
with dietary breadth (Dixon et al. 1997; Kantomaa & Rönning 1997), the mandible is an ideal
structure for studies of evolutionary change in the phenotype. It is reasonable to suspect that the
mandible would also be among the first structures to respond to novel environmental conditions,
as has been recently demonstrated in avian systems, where a combination of environmentally-

65

induced variation and modular development have favored a rapid adaptive divergence in size and
shape of the house finch beak, following expansion into novel environments (Badyaev 2009,
2010).

The mammalian mandible as an evolutionary model system
In addition to its ecological relevance, the mandible has become a popular model system for the
study of the origin of phenotypic variation in general (Atchley & Hall 1991; Hall & Miyake
2000; Klingenberg et al. 2003), which is likely a reflection of the relative simplicity of its largely
modular early development (Hall & Miyake 2000) and its amenability to theoretical and
experimental functional research (Herring 1993). To understand the link between the modular
developmental origin of the mandible and its evolutionary potential, Atchley & Hall (1991)
proposed a quantitative genetic model that describes the sources of variation in the mammalian
mandible. In that model, selection acts on variation caused by genetic and epigenetic factors, and
maternal and environmental effects refine and reshape the different regions of the mandible
during ontogeny, producing a functionally integrated structure that is not only adapted to the
organism’s immediate environment, but that is also partly shaped by it. In this context, the
mandible is seen as a complex structure, i.e., comprising an array of functionally and
developmentally integrated components. In the model’s most basic form, a total of six major
partitions have been postulated for the mandible (Atchley & Hall 1991; Atchley 1993; Hall &
Miyake 2000), which corresponds to the mesenchymal condensations that give rise to the
coronoid, condyloid, and angular processes, the ramus, and the incisor and molar alveoli in the
developed dentary bone (Fig. III.1). According to this model, these basic units of the mandible
may vary independently or become integrated by function and development (Zelditch et al.

66

2009). The mandible thus provides a theoretical framework that can take into account potential
biases introduced by function and development in analyses of patterns of variation and
diversification due to ecological pressures.
Because the development of the mammalian mandible has an important epigenetic
component, the interpretation of mandible variation can be very challenging; plasticity,
adaptation, and chance are all potential causes of phenotypic divergence. The mandible is often
regarded as highly plastic, because variation in its form has a large environmental component
(Hanken & Hall 1993). In mice and rats, for instance, a large proportion of the phenotypic
variability observed in some traits of the adult mandible (e.g., over 70% in the mouse coronoid
area, and over 50% in the rat coronoid and condyloid areas) has been shown to be caused by
residual environmental variance, likely associated with diet and other ecological factors (Atchley
1993). Also, it is possible to alter the shape of the mandible experimentally by raising animals
under different diets (Corruccini & Beecher 1982; Bouvier & Hylander 1984; Ulgen et al. 1997;
Liu et al. 1998; Kiliaridis et al. 1999; Maki et al. 2002; Renaud et al. 2010), inducing functional
alterations of masticatory muscles (Ghafari & Heeley 1982; Bresin & Kiliaridis 2002; Renaud et
al. 2010), or allowing experimental subjects to develop under colder than normal temperatures
(Steegman.At & Platner 1968).
At a macroevolutionary level, adaptation, especially to different diets, as suggested by the
diversity of mandible shapes, has been postulated as one of the main causes for the radiation of
platyrrhine primates (Anapol & Lee 1994), bats (Nogueira et al. 2009; Monteiro & Nogueira
2010), and rodents (Michaux et al. 2007; Samuels 2009). In contrast, random factors like founder
effects and subsequent genetic drift due to reduced population size have also been causally
linked to the emergence of large phenotypic differences in the shape of the mandible among

67

populations inhabiting small islands or otherwise isolated habitats (LeBoulenge et al. 1996;
Renaud & Auffray 2010).
Even though adaptation and plasticity are often treated as alternatives to explain
phenotypic differences across environments, they do not constitute mutually exclusive
explanations for the origin of such divergence. On the contrary, developmental plasticity is
currently broadly seen as a possible link between adaptation and evolutionary diversification
(West-Eberhard 2003; Angers et al. 2010; Pfennig et al. 2010; Young & Badyaev 2010). Thus,
adaptation of organisms to new habitats, especially in fluctuating environments, might be
facilitated by the presence of highly plastic responses to changing conditions (Lande 2009;
Canale & Henry 2010), much in the same way as it is facilitated by greater standing genetic
variation.
The idea that plasticity plays an important role in diversification has found some
opponents, who argue that although there is evidence for developmental plasticity, its central role
in the evolution of phenotypic novelty lacks support (de Jong & Crozier 2003). Developmental
plasticity has also been misinterpreted by some as neo-Lamarckism, as pointed out by Pfennig
and collaborators ( 2010). Currently, several epigenetic mechanisms are known to regulate gene
expression (Angers et al. 2010; Katsnelson 2010; Pfennig et al. 2010), and a role of phenotypic
plasticity in morphological divergence has been suggested for house finches (Badyaev 2009) and
shrews (Young & Badyaev 2010).
This study examines the pattern of variation in size and shape of the mandible of
Peromyscus leucopus in relation to lineage divergence and geographical range expansion. The
central aim is to address whether the rapid geographic expansion of populations of this species in
the northern Great Lakes region has been accompanied by rapid morphological divergence, and

68

whether the pattern of divergence suggests a neutral process or reflects the influence of local
eco-geographical factors. Shape was measured and analyzed using geometric morphometric
techniques. Because shape variation is expected to reflect differences in diet, a nested MANOVA
design was used to assess the potential role of non-dietary factors (i.e., history, genetics) on the
morphological differentiation of these populations. In addition, a hypothesis of phenotypic drift
was assessed by studying the strength of the association between matrices of morphometric and
genetic distances, with the latter based on a neutral marker (the mtDNA control region). Results
support the hypothesis that differences in the shape of the mandible between old and newly
established populations are related to, and plausibly caused by, the expansion process.
Additionally, the absence of evidence of phenotypic drift suggests that exposure to new
environments might contribute to the emergence of novel phenotypes. These results further
suggest that a plasticity-mediated process such as Baldwin effect and/or local adaptation might
play a causal role in the geographic diversification of Peromyscus leucopus in the northern Great
Lakes.

Materials and Methods
Samples
Mandibles were sampled from 24 of the 39 populations (Fig. III.2) included in the genetic
analyses (Chapter II). Only mandibles from adult specimens (i.e., with the third molar erupted)
and in excellent condition were included. Selection of populations was based on availability of
specimens in museum collections. Populations analyzed comprised 12 localities from WI, 5
localities from the LP, and 7 localities from the UP. With few exceptions, 10 specimens were
measured per locality, for a total of 227 specimens. Samples were obtained from the collections

69

of Michigan State University Museum, University of Michigan Museum of Zoology, and
University of Wisconsin at Stevens Point Museum of Natural History (Table A.III.1).
Geographic characteristics of samples, such as region (i.e., WI, LP, UP) and geographic
midpoint of collection sites, were defined using the same criteria as in previous genetic analyses
(Chapter II). Samples were grouped by lineage according to the partitions obtained from
SAMOVA (Chapter II), which are also the monophyletic groups observed in the haplotype
genealogy (e.g., Fig. II.4 in Chapter II). Finally, sampling localities were classified as ―old‖ or
―new‖ with respect to the timing of establishment of their corresponding populations based on
whether the population was known to be present at the site for more than 30 years, or was
discovered within the last 30 years, respectively. All specimens included in this study were
collected within the last 30 years.

Morphometric Data
Variation in morphological shape and size was analyzed using 2-D landmark coordinates placed
on homologous sites of digital images of the mandible. For each specimen, hemi-mandibles were
carefully separated prior to photography; only the lateral (buccal) surface of the right hemimandible was photographed. Hemi-mandibles were photographed with the dentary parallel to the
focal plane and against a contrasting background. In those few cases where the right side was
severely damaged, the left hemi-mandible was used instead. A reference ruler was included in all
photographs. To minimize the effect of measurement error, all photographs were taken by the
same person (B. Lundrigan), and each mandible was positioned and photographed twice;
analyses were then based on the average of the landmarks extracted from the two images.

70

Photographs were acquired using a Fuji FinePix S5 Pro digital camera with a 70mm Sigma
macro lens.
A total of 15 2-D landmarks and 5 curves (Fig. III.1) were obtained from the digital
images using TPSDig2.1 (Rohlf 2006). Curves were sub-sampled to obtain 20 regularly spaced
pseudo-landmarks (Bookstein 1997), using the software SemiThinner (Marquez 2006). Pseudolandmarks were allowed to slide along their respective curves to minimize the Procrustes
distance between each individual configuration and the mean (Zelditch et al. 2004), using the
software SemiLand6 (Sheets 2009). After sliding, pseudo-landmarks were considered semilandmarks (Bookstein 1997). Semi-landmarks differ from ordinary landmarks in that their
homology can only be established in relation to a smooth curve or edge, where characteristic
homologous features are normally absent. For that reason, variation in the spacing of semilandmarks is arbitrary and thus minimized, resulting in semi-landmarks being constrained to vary
only in the direction that is perpendicular to the margin (Zelditch et al. 2004).

Geometric Morphometric Methods
Shape is defined as the component of form that remains after removing information about
location, scale, and orientation (Kendall 1977; Bookstein 1991; Zelditch et al. 2004). In
landmark-based geometric morphometrics, these factors are removed by means of Generalized
Least Square Procrustes superimposition (Rohlf & Slice 1990; Zelditch et al. 2004). In this
procedure, each individual landmark (and semi-landmark) configuration is (1) translated so that
it is centered at the Cartesian origin; (2) scaled so that its centroid size equals one, where
centroid size is computed as the square root of the sum of square differences between each
landmark and the centroid of the configuration (i.e., the average position of all its landmarks);

71

and (3) iteratively rotated so as to minimize the sum of squared differences between the
landmarks of that specimen and those of the sample mean, which is also iteratively computed.
Following Procrustes superimposition, each configuration is mapped as a single point in a
multidimensional space, termed Kendall’s shape space, which possesses a non-Euclidean
geometry.
Following superimposition, and prior to their use in multivariate statistical analyses,
shape-space data are customarily projected onto a Euclidean space that is tangent to the shape
space at a predetermined reference point, which is usually the mean shape (Rohlf & Slice 1990;
Zelditch et al. 2004). This projection can be accomplished by transforming each individual
configuration of landmarks into a collection of residual values relative to the reference (i.e.,
mean) configuration after superimposition, so that for k landmarks, 2k of these Procrustes
residuals are obtained (Dryden & Mardia 1998), which for most practical purposes can be treated
as Euclidean shape variables. Procrustes residuals were calculated in Matlab® (The Mathworks
2006); implementations of the algorithms and equations are provided in Dryden & Mardia
(1998).
The extraction of shape information from configuration data results in the loss of four
degrees of freedom for 2-D landmarks: one from scaling by a scalar factor, two from translation
along x and y axes, and one from rotation by a scalar angle. Additionally, each semi-landmark
contributes only a single degree of freedom, since the other is lost during optimization of its
position along a curve (Bookstein 1997). Consequently, the number of degrees of freedom of the
set of Procrustes residuals is 2k + l – 4, where k and l denote the number of landmarks and semilandmarks, respectively. This creates a discrepancy between the number of variables (2k + 2l)
and the degrees of freedom of a data set, and results in covariance matrices of shape variables

72

that are not of full rank. These matrices are not usable in statistical analysis seeking to quantify
the magnitude of group mean differences relative to within-group variation, because such
analyses (e.g., MANOVA and CVA) require the inversion of the pooled within-group covariance
(Krzanowski 2000). In morphometric applications, this discrepancy in degrees of freedom is
dealt with by projecting landmark data onto a space of the correct dimensionality, such as the
principal warps of the bending energy matrix (Bookstein 1991), or by using generalized inverse
matrices in statistical analyses (Dryden & Mardia 1998). For this study, data were transformed
by projecting them onto a space defined by the first 2k + l - 4 eigenvectors of the covariance
matrix of the Procrustes residuals obtained from a superimposition of the combined set of all
samples. Thus, the final number of degrees of freedom of samples used in this study, applying
the formula above, was (15×2) + 20 – 4 = 46.
To obtain a summary metric for the overall difference between the two lineages
examined in this study, I used partial Procrustes distance, computed as the squared root of the
minimized sum of square differences between homologous landmarks (Dryden & Mardia 1998).
This distance represents a close approximation to the actual distance between each specimen and
the sample mean in shape space (Zelditch et al. 2004).
Landmark data were also used to measure mandible size, estimated, as it is customary in
geometric morphometric applications, as centroid size (CS, Dryden & Mardia 1998), defined
above. Given that CS is a one-dimensional measure, statistical analyses of size variation were
carried out using the univariate equivalents of the methods used to analyze shape.

Statistical Analyses

73

Statistical analyses were conceived as a hierarchy of geographical and genetic factors potentially
able to explain the morphological differences observed between the mandibles of ―old‖ and
―new‖ populations. This hierarchical approach seeks to address three questions, two of which are
primarily concerned with broad patterns of geographical and historical variation in mandible
shape, and the other with more specific associations between genetic and morphological
differentiation:
1. Are mandible phenotypes distinguishable between lineages, between old and new
localities within lineages, and/or among populations within lineages?
2. Are differences in phenotypic means between old and new populations characterized by
the emergence of novel directions of variation relative to the directions observed among
old populations alone?
3. Are morphological differences among populations congruent with, and therefore
explained by, their phylogeographic pattern based on a neutral genetic marker?
Prior to carrying out tests to address these questions, sexual dimorphism in shape and size
was ruled out by testing for the effects of population, sex, and their interaction using MANOVA
and ANOVA, respectively. Because the latter two effects were found to be non-significant (P >
0.05), males and females were pooled for subsequent statistical analyses.

Genealogical and Historical Effects
I used a nested MANOVA design to test whether there are significant differences in mean shape
(multivariate dependent variable) between the LP-EUP and WUP-WI lineages, between old and
new localities within lineages, and among populations within old and new localities (independent
variables), corresponding to the model:
74

where L denotes the effect of lineage (i.e., LP-EUP vs. WUP-WI); H denotes the historical effect
(i.e., old vs. new populations), nested within lineage; P denotes the population effect, nested
within lineage and historical groups; and μ and ɛ represent a mean vector and error matrix,
respectively (Sokal & Rohlf 1995).
In order to apply MANOVA to test differences in shape among groups, two requirements
need to be met: the measurement space must be Euclidean, and the intra- and inter-group
covariance matrices must be non-singular (i.e., the number of variables included in analyses must
equal the number of degrees of freedom of the variation). As described above, Procrustes
superimposition and tangent space projection ensure that both of these conditions are met when
using landmark data. To test the null hypothesis of no differences in mean shape among groups, I
used the Hotelling-Lawley trace criterion and its associated F- and P- values, which is one of the
multiple testing criteria often used in MANOVA. This statistic is computed as the sum of the
eigenvalues of BW-1, with W and B denoting the pooled within- and factor-specific betweengroup sums of squares and cross-products (SSCP) matrices, respectively. The Hotelling-Lawley
trace thus estimates the total variance captured by a test factor after correcting for intra-group
variation. The F-value is also a function of W and B, providing a measure of the magnitude of
the variance of the dependent variables that is explained by the independent variables (Zelditch
et al. 2004). The effects of lineage, expansion history, and population on centroid size (CS)
variation were analyzed using a univariate equivalent to the nested MANOVA model used for
shape data, i.e., nested ANOVA.
In addition to testing for statistical significance, the effect size of each factor of the model
was estimated. Effect size for CS was estimated as the value of R2 for each effect in the ANOVA
75

model. For shape, an effect size measurement analogous to R2 was calculated using the
following procedure (E. Marquez, personal communication). The least squares means
corresponding to the predicted values for each treatment and factor in the nested MANOVA
model were pooled, and the covariance matrix of these means was used to extract all
eigenvectors with a corresponding non-zero eigenvalue. For the lineage effect, this
decomposition yields one eigenvector, equivalent to the difference between the LP-EUP and
WUP-WI mean vectors, whereas for the historical effect within lineage and the effect of
population nested within lineage and historical groups, there are 3 and 23 non-zero eigenvectors,
respectively. Next, the original shape data were projected onto these eigenvectors and the
variance in each case was calculated. The estimated multivariate effect size is equal to
∑

, where the numerator is the accumulated variance explained by the eigenvectors vi

of the values predicted for the jth effect, and the denominator is the total shape variance,
computed as the trace of their covariance matrix. The effect sizes in this analysis are not additive
(i.e., they are not intended to add to 100).
Principal Component Analysis (PCA) and Canonical Variate Analysis (CVA) were used
for visual inspection of patterns of geographical and historical variation among populations.
PCA finds linear combinations of the original variables that maximize the amount of information
(i.e., variation) present in a homogeneous sample. In this application, PCA was used on
population means to produce shape deformation grids (Bookstein 1991) depicting the main
directions of mandibular diversity observed in the study. CVA uses group information to find
linear combinations that optimally explain variation across means, maximizing the ratio of
among- to pooled within-group variation. This is accomplished by computing the eigenvectors of

76

the matrix W-1B, where B and W are as defined above (Krzanowski 2000). Specimen scores on
canonical axes were obtained by projecting Procrustes residuals onto these eigenvectors. CVA
on shape data was used to visualize the distribution of differences among the 24 populations
included in the study. These canonical scores were then grouped according to historical and
geographical categories to visually assess the extent to which differences among populations are
explicable in terms of lineage membership or history of population expansion.
To further investigate morphological differences between old and new populations, I
addressed the question of whether variation in the dimensionality of population mean differences
within each lineage is related to the timing of geographic expansion. This analysis complements
results from MANOVA, which indicates whether at least one of the sampled groups differs from
the rest in terms of magnitude, without considering the direction of such differences. There are
two possible outcomes: (1) phenotypic differences appearing during or after geographic
expansion might be accompanied by an increase in dimensionality; or (2) phenotypic differences
might accumulate along the same dimensions of variation as already present among old
populations. To assess these possibilities, I estimated and contrasted the dimensionality of mean
population differences between two data sets, one containing only old populations and another
containing both old and new populations combined. By treating old vs. (old + new) as two
different blocks, I can address the question of whether new populations introduce morphological
variation consistent with already existing variability, in which case no change in dimensionality
would be expected, or whether the addition of new populations results in variation in new
features, as indicated by the addition of new dimensions. The rationale for this comparison is that
differences in dimensionality between the two blocks of populations would have to be accounted
for by the new populations, wherein additional dimensions represent novel directions of

77

mandible variation. To estimate dimensionality, I used Bartlett’s statistic (Krzanowski 2000),
which tests for equality of the q smallest CVA eigenvalues, hence defining dimensionality as the
maximum number of canonical axes that capture a distinct (i.e., non-error) proportion of the total
shape variation (Anderson 1963). Bartlett’s statistic is based on a likelihood ratio test and
approximately follows a χ2 distribution (Krzanowski 2000).
Nested MANOVA statistics were computed using SAS procedures. GML and IML (SAS
Institute 2004), and principal components, canonical axes, and Bartlett’s statistics were
computed using Matlab® functions PCACOV and MANOVA1 (The Mathworks 2006).

Genetic Effects
The tests described above address geographic and historical variation in phenotypic means,
asking whether differences found among populations can be explained by geographic and/or
historical effects. However, even if a hypothesis invoking those effects is supported, a thorough
examination of underlying patterns of genealogical associations is still necessary to fully account
for possible sources of morphological variation. Because geographic, historical, and genealogical
effects are not mutually exclusive, morphological similarities might merely mirror genealogical
associations among populations, supporting a hypothesis of neutral divergence (Monteiro et al.
2005; Nogueira et al. 2005; Perez et al. 2009).
To assess the hypothesis that patterns of inter-population differentiation based on neutral
genetic markers and on morphometric data are similar enough to be considered causally linked, I
computed correlation coefficients between genetic and morphological distance matrices, and
carried out Mantel tests to estimate the statistical significance of those correlations (Sokal &
Rohlf 1995). In these tests, observed Pearson’s correlation coefficients between the elements of

78

the two matrices are compared to a distribution of equivalent correlations computed after the
rows and columns of one of the matrices have been randomly permuted. This procedure
simulates the distribution of correlation values expected under the null hypothesis of no
association between matrices; P-values for significance tests are computed as the proportion of
correlations between randomized matrices that equals or exceeds the observed correlation (Dietz
1983; Manly 2007). I computed a total of 1,000 random permutations for each Mantel test; these
were carried out using the CORRCOEF and RANDPERM Matlab® functions (The Mathworks
2006).
Mantel tests were based on two kinds of genetic metrics, pairwise Nei’s distances (Nei's
D, Nei 1987) and Fst estimates (Reynolds et al. 1983), which represent distinct but similarly
informative genetic criteria. Nei’s D assumes that a population is in mutation-drift equilibrium; it
is calculated as DNei = -ln I, where I is a measure of the identity of genes from two populations.
Nei (1973) pointed out that this measure of genetic distance has some interesting properties,
including that it is (a) related to the kinship coefficient, (b) linearly correlated with evolutionary
time if the substitution rate is constant, and (c) linearly related to geographic distance under the
assumptions of some migration models. In contrast, Fst can be used as an indicator of short-term
genetic distance, where divergence between populations is assumed to have been caused by drift
only. Assuming short divergence times, Fst-based pairwise distances were computed as DFst = log(1-Fst) ~ t/N, where t is the number of generations since divergence and N is the size of
haploid populations (Reynolds et al. 1983). Estimates of both distances were computed using
Arlequin 3.11 (Excoffier et al. 2005), with the output formatted as distance matrices.
To estimate pairwise morphological distances between populations, I computed
Euclidean distances between population means obtained from the shape variables projected onto

79

a 46-dimensional planar tangent space as described above (Table A.III.4). These distances
approximate the corresponding distances between mean shapes in shape space. A Euclidean
metric was also used to compute distances among population means for centroid size.
Graphical visualization of the information in distance matrices was summarized using
dendrograms created with the UPGMA method, which uses group averages to link clusters
(Rencher 2002). These were produced using the Matlab® functions DENDROGRAM,
LINKAGE, and COPHENET with the genetic and morphometric distance matrices described
above.

Results
Effects of Genealogy and Expansion History on Morphological Variation

Mean Shape Variation Between Lineages and Among Populations
Results from the nested MANOVA on shape variables indicate that all tested factors are highly
significant and have moderate to large effects (Table III.1). Effect size estimates indicate that the
effect of lineage explains almost 9% of the overall variation in shape, whereas historical effect,
which is nested within lineage, explains over 24% of the total shape variation. Differences
among populations, which include the cumulative effect of lineage and history, explain over 85%
of the differences in the shape of the mandible of Peromyscus leucopus in the geographical area
covered by this study. The remaining 15% of the variation is not explained by any of the factors
included in this nested model, and might be related to variation within populations. In contrast,
the analysis based on CS indicates that none of the factors tested in the nested design have a
significant effect on the pattern of mandible size variation (Table III.2).

80

The first four principal components computed from mean population shapes jointly
explain 69% of the total variation among those means, suggesting a relatively high
dimensionality of population variation, with no major trends. Plots of PC scores (Fig. III.3)
suggest that the major axes of variation are unrelated to any of the factors of interest in this
study; i.e. none of these axes provides unambiguous support for differentiation between the two
lineages or between old and new populations.
Shape deformations implied by PCs 1 through 4 are shown in Figures 4a-d. PC1 (which
accounts for 26% of the total variation) appears to be associated with contraction /expansion of
the anterior margin of the coronoid process and caudoventral angular process, and indicates that
these structures are highly variable with respect to mean shape. These changes are coupled with
expansion /contraction of the caudal end of the masseteric ridge, which corresponds to the base
of the incisor alveolus. Two UP populations belonging to the WUP-WI lineage, Marquette (MA)
and Delta (DE), exhibit the most extreme shapes along this axis (Fig. III.3, top). PC2 (24%)
captures a pattern of variation that simultaneously involves elevation /reduction of the condyloid
process, expansion /contraction near the caudal margins of the coronoid and angular processes,
displacement of the masseteric ridge, and expansion /contraction of the anterior incisor alveolus.
Although the most contrasting populations along PC2 correspond to different lineages (Crawford
–CR– and Schoolcraft –SC–), there is substantial overlap of the two lineages along this axis (Fig.
III.3, top). Inspection of PC3 (11%) and PC4 (8%) suggest covariation between the coronoid and
angular processes, and between these processes and the alveolar region, respectively (Figs. III.4
c and d), in neither case clearly separating the groups of interest (Fig. III.3, bottom).
Canonical axes are directions of maximal differentiation between population means
relative to within-population variation and thus show different patterns of population distribution

81

than PCA axes. According to Bartlett’s test, only four CVs capture a significant proportion of the
shape variation among populations, and thus no additional axes were computed. These four
canonical axes account for only 46% of the standardized variation among groups. In Figure III.5,
CV scores are arranged by population (i.e., by county of collection) and region, making it
possible to visually assess whether morphological differences among populations are structured
by these biogeographical factors. The box plot of population scores on the first canonical axis
(16.4% of the variation) suggests a contrast between old and new populations, with the most
extreme values corresponding to newly invaded localities (e.g., Schoolcraft, Delta, Bayfield, and
Chippewa). Subsequent canonical axes do not show clear patterns that can be related to the
difference between old and new localities.
Visual inspection of deformation grids contrasting old and new populations shows that
the patterns of morphological variation in the mandible differ between the two lineages (Fig.
III.6). In the LP-EUP lineage, change from old to new involves moderate expansion of the
ascending ramus, accompanied by contraction near the caudal end of the masseteric ridge and in
the region of the angular process (Figs. III.6 a and III.7). In the WUP-WI lineage, change from
old to new is characterized by contraction of the 1st molar alveolus and base of the coronoid
process, increased curvature of the caudal margin connecting condyloid and angular processes,
and expansion of the incisor alveolus (Fig. III.6b). Distances between the mean shapes of old and
new populations equal 2.81×10-4 and 3.29×10-5 squared Procrustes units for the LP-EUP and
WUP-WI lineages, respectively. Table III.3 shows proportional contributions of the major
developmental blocks of the mandible to these distances, computed using the landmarks
contained within each region (Fig. III.1). These contributions range from 10% (coronoid process)
to 29% (masseteric ridge) in the LP-EUP and from 7% (incisor alveolus) to 28% (condyloid

82

process) in the WUP-WI, and suggest differences in the pattern of morphological divergence
between these lineages. Shape change in the LP-EUP is focused primarily on the masseteric
ridge and angular process, whereas in the WUP-WI, the condyloid process and 1st molar alveolus
undergo the most change.

Changes in Shape Variation Dimensionality during Population Expansion
Comparison of estimates of dimensionality of shape variation in samples of means from old vs.
(old + new) populations show a significant increase in dimensionality of the whole mandible
only for WUP-WI. When the anterior and posterior regions of the mandible are analyzed
separately, both lineages show an increase in dimensionality associated with geographical
expansion (α = 0.0167, after Bonferroni’s correction; P < 0.01, for all comparisons).
In the WUP-WI lineage, Bartlett’s test estimated an increase from 2 to 3 dimensions on
addition of the new populations to the old ones, when the whole mandible was analyzed.
Repeating this test by removing new populations one at a time does not change this result,
suggesting that the gain in dimensionality is not caused by any single outlier population. In
contrast, in the LP-EUP lineage Bartlett’s test did not show a change in dimensionality for the
whole mandible (3 for old and (old + new) populations). However, there is an evident shape
change between old and new populations in the regions of the angular and coronoid processes
and caudal masseteric ridge (Fig. III.6a). Separate Bartlett’s tests of the posterior (angular,
condyloid, and coronoid processes, and masseteric ridge) and anterior (incisor and molar alveoli)
regions of the mandible found an increase of from 1 to 2 dimensions in the posterior mandible
and unchanged dimensionality (=1) in the anterior mandible. Similarly, an analysis of
dimensionality of anterior and posterior regions of the mandible in the WUP-WI lineage revealed

83

an increase of dimensionality in both partitions, with the posterior mandible displaying a more
striking change (from 1 to 3 dimension) than the anterior mandible (from 1 to 2 dimensions).
Together, these results support the hypothesis that geographic expansion in both lineages was
accompanied by the emergence of novel directions of morphological change, a pattern involving
the entire mandible in the WUP-WI lineage, but restricted to the posterior region in the LP-EUP
lineage.

Genetic Effects
Correlations Between Genetic and Morphological Patterns
Mantel tests did not find a statistically significant association between genetic and morphological
distances (Table III.4). Neither pairwise Nei’s D nor pairwise Fst’s were significantly correlated
with pairwise Euclidean distances between shapes (P > 0.0125), although for both genetic
distance metrics, the correlation coefficients are of moderate magnitude. The incongruence
between genetic and morphometric patterns is further evinced in the UPGMA dendograms based
on these distance matrices (Figs. III.8-III.10). Together, these data suggest that the differences in
mandible shape detected by MANOVA cannot be explained by genealogical relationships.
Similarly, differentiation in mandible size (CS) is not significantly correlated with genetic
distance (full dataset r = 0.0423, P =0.5285; LP-EUP-only r = 0.0677, P = 0.4945; WUP-WIonly r = 0.1844, P = 0.2178).

Discussion
The questions addressed in this study have been organized as a hierarchy of potential
explanations for observed patterns of geographic and historical variation in mandible

84

morphology. Thus, at a large scale, analyses focused on overall comparison of shape means
among populations, whereas at a smaller scale, comparisons took into account patterns of
differentiation along historical and genetic dimensions of variation, with the ultimate aim of
explaining phenotypic differences between old and newly established populations. This
hierarchical organization offers an innovative approach to search for possible causes of
underlying patterns of morphological differentiation associated with the process of geographic
expansion, by dissecting layers of information from broader to more specific potential causal
factors.
The broadest level in the hierarchy addressed whether population means are different,
within the framework of a nested MANOVA design. Results in this case suggest that not only
are populations different overall, but there are also significant effects associated with lineage and
expansion history. The effect of lineage accounts for only 9% of the shape variation observed
among individuals; the effect of history, which takes into account the effect of lineage, explains
over 24% of shape variation. In other words, over 24% of the shape variation among individuals
is explained by membership in old or new populations that belong to the LP-EUP or to the WUPWI lineage. Most of the shape variation explained by these factors involves the masseteric region
(in the LP-EUP lineage), and condyloid process (in WUP-WI). While not the focus of this study,
the largest effect on shape variation explained by this nested model was for differences among
populations, which accounts for 85% of the total variation in shape among individuals. A large
value is expected here, as this factor is the most inclusive source of variation in this nested
MANOVA; it includes the cumulative effects of lineage and history among populations, in
addition to other sources of geographical variation not examined in this study, such as local
ecology.

85

Whereas a large population effect is expected in the absence of detailed information
about local environments, the relative magnitude of the effect of lineage and expansion history
within lineages suggests that the process of geographical expansion undergone by these
populations in recent decades has resulted in morphological divergence comparable to that
between the two lineages, which are estimated to have split circa 35,000 YBP (Chapter II).
These results are consistent with a scenario in which abundant morphological diversification
accumulated independently within each lineage, with variation accruing along multiple
phenotypic dimensions until intra-lineage variation far exceeded inter-lineage divergence, a
pattern that is clearly illustrated in the PCA and CVA results.
The next hierarchical level of explanation refined the previous question (whether old and
new population means are different) by asking whether diversification has accrued along a single
direction of divergence (as expected of a geographic or ecological cline, or where developmental
biases restrict the direction of variation), or alternatively, whether the expansion process has
produced novel directions of variation. Results from Bartlett’s test comparing the dimensionality
of shape variation between old vs. (old + new) populations indicate that the occupation of new
environments during expansion has been accompanied by divergence along novel directions of
variation, suggesting the possibility of adaptive processes brought about by exposure to novel
selection pressures.
Graphical comparisons of mean shapes between old and new populations (Figs. III.6 and
II.7) suggest that the expansion process has affected the mandible differently in the two lineages.
In the LP-EUP, the angular process and masseteric region contribute the most to the divergence
between old and new populations; both are points of insertion of masticatory muscles, and they
jointly account for about 54% of the difference between group means. In contrast, in the WUP-

86

WI, the largest contribution to the divergence between old and new populations comes from the
condyloid process, in the posterior mandible, and 1st molar alveolus, an anterior structure.
Jointly, these regions account for about 51% of the difference between old and new populations
(Table III.3, Fig. III.6).
Despite this difference in the pattern of variation between old and new populations, a
common feature of both lineages is the tendency for variation to be focused on the region of
insertion of masticatory muscles, which includes the posterior end of the masseteric ridge, and
the coronoid, condyloid, and angular processes (Atchley et al. 1992) (Table III.3). With few
exceptions, most studies have also found that mandibular processes are the most variable region
of the mandible. Some studies have found that a single mandibular process bears most of the
variation (Cheverud et al. 1991); in other cases, a combination of processes shows most of the
variation (Klingenberg et al. 2001a; Klingenberg et al. 2001b; Renaud & Michaux 2004; Perez et
al. 2009); and in other instances, all threes processes combined concentrate most of the variation
(Atchley et al. 1985; Duarte et al. 2000; Monteiro et al. 2005). However, other studies have also
found a pattern of variation involving the anterior region of the mandible as well. For instance,
Atchley et al. (1992) found that posterior structures like the coronoid and condyloid processes
and the lower edge of the angular process show relatively large additive and phenotypic variance
in random-bred house mouse, but also that the incisor alveolus (an anterior structures ) shows
similarly large variance. A contrasting pattern was reported by Cheverud et al. (1991), who
found in a comparison of 10 inbred mouse strains that most of the variation in the mandible was
restricted to the anterior mandible, and to a lesser extent, the coronoid process.
Several other studies have also reported that variation of the mandibular processes
contributes significantly to shape variation associated with geographic differentiation. Duarte et

87

al. (2000) found that differences between spiny rats (Thrichomys apereoides) from Palmeiras
and other populations are largely restricted to the coronoid and angular processes, and the
alveoli. Monteiro et al. (2005) expanded the number of localities, skull views, and ecological
variables in an analysis of the same group, and described similar changes in the mandible, which
they attributed to a latitudinal ecological gradient. Similarly, Renaud (2005) found that the main
axes of variation among populations of the European wood mouse (Apodemus sylvaticus) in
Germany involve variation in the coronoid and angular processes. In a broader geographic study
of this same species of mouse, which included mainland and insular populations over a
latitudinal gradient, Renaud & Michaux (2007) found that differences among insular populations
were mostly due to variation in the alveolar region in addition to the coronoid and angular
processes. Because the pattern of variation of insular and mainland mice was different, these
authors concluded that island ecology plays a key role in the morphological variation of the
insular mice.
The prominence of the mandibular processes and masseteric ridge in overall shape
variation is likely related to their functional role as insertion sites and loading regions for
important masticatory muscles. These include the masseter superficialis and lateralis (attached to
the angular process), masseter medialis and temporalis (coronoid process), masseter lateralis
(masseteric ridge), and pterygoideus externus (condyloid process; Rinker 1954). The close
interdependency between the development of bone and muscle tissue and the impact on
morphology of mechanical forces generated during musculoskeletal (e.g., masticatory) function
are well established (Herring 1993). Thus, the observation that most of the variation among
populations is concentrated in the region of insertion of masticatory muscles is plausibly a
consequence of dietary differences, in particular the physical properties of food, such as

88

hardness, size, and/or manageability. Several studies have measured the effects of dietary
consistency on masticatory function and development of the mandible in rodents (Watt &
Williams 1951; Bouvier & Hylander 1984; Myers et al. 1996; Ulgen et al. 1997; Kiliaridis et al.
1999; Maki et al. 2002), as well as humans and non-human primates (Bouvier & Hylander 1981;
Taylor 2002; Varrela 2006). In most of these studies, differences in diet had a significant effect
on mandible shape. However, in the few studies that actually measured that effect, the magnitude
represented a small portion of the total variation. For example, Myers and colleagues (1996)
found a significant effect of diet (i.e. hard versus soft dietary regimens) on cranial shape in
prairie mice (Peromyscus maniculatus bairdii), but the magnitude of the effect was extremely
small, suggesting that the effect of diet on mandible shape was almost negligible. Bouvier &
Hylander (1984) did not find statistically significant differences in mandibular shape between
rats fed soft versus hard diets, although animals on a soft diet had maxillary and mandibular
measurements that were smaller on average than those on a hard diet regime. In a similar study
by Watt & William (1951), differences between soft diet and hard diet groups were also small,
accounting for only 1 to 3% of the observed differences in mandible morphology. Maki et al.
(2002) reported a reduction of 6% of the posterior jaw height of laboratory mice fed a soft diet,
which was likely manifest as a reduction in both coronoid and angular processes. In another
study with laboratory mice, Renaud et al. (2010) investigated the effects of masticatory function
on the shape of mandible. Although they did not measure the effect of diet on mandible shape
variation directly, they reported that PC2, which accounted for most of the differences associated
with food consistency, captured around 14% of the total variation.
Although a comparison across studies is made difficult by differences in methodological
approaches, there is support for the hypothesis that diet affects the shape of the mandible, most

89

likely through masticatory activity, and that this effect is predominantly observed in the
processes of the mandible. In my comparison of P. leucopus populations, change in the angular
(LP-EUP) and condyloid (WUP-WI) processes was characterized by contraction (relative to
surrounding morphological regions) in the newly established populations, suggesting a shift to a
softer diet (Maki et al. 2002; Renaud et al. 2010). Such a shift would presumable require lower
mastication forces, and therefore, smaller muscles. One possibility is that insects have become
more prominent in the diet of UP mice than in mice from WI or LP, which probably rely more on
hard-shell seeds. Peromyscus leucopus is an omnivorous rodent whose diet varies seasonally and
geographically. Populations inhabiting deciduous forest in the Eastern United States feed on
arthropods and insects, which may represent over 45% of their diet. However, the relative
proportion of nuts and seeds increases when those become available, reaching about 30% in
winter (Wolff et al. 1985). Another possibility is that the physical properties (e.g., hardness) of
seeds and nuts differ between northern and southern localities.
None of these studies has explicitly measured changes in the dimensionality of mandible
variation. However, results from the present study demonstrate that the dimensionality of shape
variation can offer valuable insight in the study of phenotypic divergence, and should be
included in comparison of populations, since regions showing a similar magnitude of variation
can in fact diverge along different directions. This is a simple consequence of the
multidimensional nature of morphological variation, and it illustrates the value of focusing on
directions of variation, rather than the variation of anatomical regions of the mandible alone. As
indicated by the nested MANOVA model used here, variation in mandible shape among
populations within lineages is substantial, accounting for 85% of the total variation in shape,
which is sufficient to obscure the differences between two long-isolated lineages, and old and

90

newly established populations within these, which account for 9% and 24% of the variation,
respectively. Yet a pattern was elucidated, indicating the emergence of novel dimensions of
divergence in an expansion process that is only a few decades old. Coupled with the fact that this
invasion has resulted in the occupancy of novel habitats, i.e., boreal forests, with their own set of
ecological challenges, it is plausible to hypothesize that ecological factors, more specifically diet,
have had a significant role in the shape diversification of the mandible of expanding populations
of Peromyscus leucopus, as they had on the radiation of old world rats and mice (Michaux et al.
2007; Samuels 2009).
The third level in the hierarchy of explanations explored in this study asks whether mean
shape differences, which were previously shown to be significant in magnitude and direction,
could be explained by genetic distances computed using neutral markers, in which case no role
for local ecological processes has to be invoked. Comparisons between phylogeographic
differentiation, based on mtDNA sequence data from the control region, and morphological
disparity, based on geometric shape measurements, found a lack of congruence. This lack of
association between genetic and morphometric distances is consistent with disparate causes
underlying these patterns, and implies that morphological differentiation in these populations is
not simply due to phenotypic drift.
Because patterns of morphometric variation observed in this study cannot be attributed to
phenotypic drift, phenotypic plasticity and selection remain as plausible explanations. The only
new population that was sampled from the LP-EUP lineage, namely Chippewa, is markedly
different from other populations in the same lineage not only morphologically, as shown above,
but also genetically (Chapter II). The specimens analyzed in this study came from an isolated
patch of red oak forest in the Hiawatha Forest in the eastern UP, which has an unusual

91

abundance of oaks compared to most other localities in the UP (P. Myers, personal
communication). It is thus possible that the divergence of this population has resulted from the
combined effects of unique ecological conditions and genetic isolation. In contrast, new
populations in the WUP-WI lineage do not show a sharp genetic differentiation with respect to
other populations within that lineage (Chapter II), thus suggesting that morphological
differentiation between old and new populations might predominantly reflect functional
differences associated with local ecology, presumably due to of phenotypic plasticity or local
adaptation.
Previous studies of wild populations of P. maniculatus (Holbrook 1982) and P. leucopus
(Elrod & Kennedy 1995) have also found morphological differentiation at a microgeographic
scale among populations inhabiting different habitats. Holbrook (1982) reported that mandibles
of mice from woodlands and grasslands had significantly different morphologies, which the
author attributed to diet, suggesting that high morphological plasticity in these generalist rodents
might have enabled them to exploit a broad variety of microhabitats. Elrod & Kennedy (1995)
compared cranial and mandibular measurements of mice from eight localities and discovered a
pattern of morphological variation that was seemingly unrelated to geographic proximity, habitat
type, or broad climatic variables. While authors in this case favored selection as the mechanism
underlying differentiation, plasticity or a joint effect of plasticity and adaptation cannot be ruled
out as explanations.
A large proportion of the studies of geographical variation in the mouse mandible
attribute differences in morphology to dietary factors. The predominance of diet as an
explanation reflects both the role of the mandible as the primary foraging structure, and the
developmental plasticity of the mandible, a structure that continues to be influenced by

92

masticatory function during postnatal life (Atchley 1993; Herring 1993; Renaud et al. 2010b).
Thus, the developmental dynamics of the mandible favor the integration of environmental
sources of variation on shape, which is a basis for plasticity (Badyaev & Foresman 2000; Young
& Badyaev 2010) and could provide favorable conditions for rapid morphological divergence.
Mammalian mandibles have been postulated to evolve rapidly, presumably because they have a
single mechanical function, and are thus not generally subject to conflicting selective pressures
(Caumul & Polly 2005). However, even though causal links between mandible morphology and
masticatory function have often been hypothesized and assumed in studies of geographical
variation, the necessary information to rigorously test these ideas is generally lacking, while
critical research is still required in order to elucidate the complexity of the interaction between
features of food items (e.g., hardness, nutritional content, mobility, toxicity) and their combined
effect on mandible form.
In summary, results of this study show a pattern of rapid and abundant differentiation in
the mandible, some of which appears to be the result of an expansion process into novel habitats.
Rapid evolutionary rates in this structure are expected for a number of reasons, among the most
important being that masticatory function, which is central for the morphogenesis of the
mandible, continues to play a dynamic role in bone remodeling throughout ontogeny. The
increased sensitivity to environmental variation that results might be translated into reaction
norms with larger gradients, which could facilitate the invasion into new habitats as average
plastic responses in the population would effectively shift approximately instantaneously into the
adaptive zone characterizing the new environment (Pfennig et al. 2010).
An equivalent adaptive process relying on the spread of an advantageous mutation would
take a greater number of generations, and might require the population to cross through adaptive

93

valleys before reaching a new peak (Pfennig et al. 2010). However, a lack of functional
constraints in the mandible (Caumul & Polly 2005) would presumably accelerate the rate of
adaptation in this structure, suggesting that both selection and plasticity might simultaneously
contribute to its rapid rate of evolution. These two processes are suitable explanations for the
differences in the shape of the mandible found between old and new populations, within each
lineage, but by no means are they mutually exclusive, as argued above.
Two distinct mechanisms have been proposed to explain phenotypic evolution mediated
by plasticity, namely Baldwin’s effect and genetic assimilation (Crispo 2007; Pfennig et al.
2010); in neither case is plasticity put forward as the direct cause of phenotypic evolution.
Instead, it is the developmental processes underlying phenotypic responses to environmental
variation that are considered potential targets of selection (Schlichting & Pigliucci 1998;
Schlichting 2004; Crispo 2007; Angers et al. 2010; Pfennig et al. 2010). Both of these
mechanisms assume the presence of environmentally induced phenotypic changes as a basis for
the action of selection.
Results of this study suggest that even though the geographical expansion of P. leucopus
in the northern Great Lakes is relatively recent, the phenotypic consequences are already
noticeable. Therefore, it appears that in this case range expansion should not be regarded as an
instance of habitat tracking. Even if diet is treated as the sole ecological factor playing a role in
mandible form variation, it remains unclear which aspects of the diet are associated with
particular directions of variation, or how distinct food properties might cancel or enhance each
other’s effect on the mandible, for all of which additional investigation is required in order to
improve our understanding of the evolutionary and environmental forces shaping the mandible of
these populations.

94

Tables

Table III.1. Results from a nested MANOVA design measuring the effects of lineage (i.e., LPEUP and WUP-WI), expansion history (i.e., old, established, vs. new, recently invaded
populations) nested within lineage, and populations, nested within old and new localities, on
shape of the mandible in Peromyscus leucopus noveboracensis. Note that Effect Size is not
additive.
Effect

Hotelling-Lawley trace

F-value

P-value

Effect Size

Lineage, Li

0.499

1.71

0.008

8.9%

Expansion history, Hij

1.479

2.52

< 0.0001

24.29%

Population, Pijk

8.166

1.39

< 0.0001

85.1%

95

Table III.2. Results from a nested ANOVA design measuring the effects of lineage, expansion
history (i.e., old, established, vs. new, recently invaded populations) nested within lineage, and
populations, nested within old and new localities, on centroid size of the mandible in Peromyscus
leucopus noveboracensis. Note that Effect Size is not additive.

Effect

Mean Square

F-value

P-value

Effect Size

Lineage, li

441.89

3.02

0.084

1.3%

Expansion history, hij

211.72

1.45

0.238

1.3%

Population, pijk

146.57

1.00

0.462

8.8%

96

Table III.3. Contribution of mandible developmental blocks to the difference between old and
new populations in each of the two lineages sampled in this study. Contributions are given in
squared Procrustes units and as percentages of the total squared Procrustes distance between the
mean shape of old and new populations (in parenthesis). Distances have been multiplied by 105.

Lineage

Angular
Process

Condyloid
Process

Coronoid
Process

Masseteric
Region

Incisor
Alveolus
(anterior
end)

LPEUP

0.709
(25.22%)

0.318
(11.30%)

0.275
(9.77%)

0.818
(29.09%)

0.372
(13.24%)

0.320
(11.38%)

WUPWI

0.065
(19.78%)

0.092
(27.87%)

0.026
(7.96%)

0.046
(14.09%)

0.023
(7.00%)

0.076
(23.30%)

97

Molar
Alveolus

Table III.4. Mantel test to assess the association between morphometric and genetic distances.
Significance level corresponding to a 95% confidence interval after a Bonferroni’s correction (α
= 0.0125).
Nei’s D
r

Fst
P

r

P

Full dataset

0.2645

0.05

0.3302

0.0260

LP-EUP

0.4989

0.1140

0.5106

0.1520

WUP-WI

0.3671

0.0930

0.4047

0.0730

98

Figures

Figure III.1. Configuration of landmarks (white circles) and semi-landmarks (black dots) sampled in this study. Shadowed areas
indicate regions with developmental or functional significance. With the exception of the masseteric ridge, these regions are defined to
roughly correspond to developmental units of the mandible (Atchley and Hall 1991).

99

Figure III.2. Geographic origin of the 24 populations examined in this research. The two-letter
acronyms are as in Table A.III.1. Black circles represent old populations, while grey ones are
new populations, according to the criterion explained in the text.

100

Figure III.3. Plots of scores of PC1 vs. PC2 (top) and PC3 vs. PC4 (bottom). Black circles are
old populations in LP-EUP and white circle is the new population in this lineage. Black and
white triangles are old and new populations, respectively, in WUP-WI. The two-letter acronyms
are as in Table A.III.1.

101

a)

b)

c)

d)

Figure III.4. Shape deformations implied by first four PCs of population means: a) PC1 (26.32%), b) PC2 (23.57%), c) PC3
(11.17%), d) PC4 (7.78%). Colors represent continuous change in surface area of infinitesimally local regions of the mandible, as
inferred via TPS interpolation, where colder colors (toward blue) are interpreted as local contraction and warmer ones (toward red) as
expansion. Numbers in the scale represent percentage of local change (e.g., -0.1 means a local contraction of 10% with respect to the
reference). Deformation vectors and grids (but not color patterns) have been magnified x5 to improve visualization.

102

Figure III.5. Box plots of canonical scores for all populations, sorted by lineage (i.e., LP-EUP,
WUP-WI). Old populations are indicated by a shadowed background and new populations are on
white (see text for details). The number of axes chosen for interpretation was based on Bartlett’s
test results, which estimated four significant dimensions for the differences among population
means. Proportion of the variation explained by these canonical axes was CV1: 16.40%, CV2:
11.17%, CV3: 10.19%, and CV4: 8.67%.

103

a)

b)

Figure III.6. Differences between mean shapes of old and new populations within (a) LP-EUP
and (b) WUP-WI lineages. Landmarks in mean shape reference of old samples are connected by
a black line; landmarks in mean shape of new population are connected by a grey line. Note the
difference in scale, showing more pronounced changes in (a). Colors represent proportional
change throughout the mandible, where colder colors (toward blue) are interpreted as local
contraction and warmer ones (toward red) as expansion. Numbers in the scale are percentage of
local change (e.g., -0.10 means a local contraction of 10% with respect to the reference). To
improve visibility of transformation grids, deformations have been multiplied by 4 in (a) and by
8 in (b).

104

Figure III.7. Trend of morphological variation within old populations and from old to new
populations. The most extreme specimens are shown in each case to facilitate visualization of the
changes already suggested by deformation grids (Fig. III.6).

105

Figure III.8. UPGMA dendrogram based on morphometric distances among populations.
Distances were calculated as pairwise Euclidean distances between population means computed
from the shape variables derived from projecting Procrustes residuals onto the first 46
eigenvectors of their covariance matrix.

106

Figure III.9. UPGMA dendrogram based on pairwise Nei’s genetic distances among
populations.

107

Figure III.10. UPGMA dendrogram based on pairwise Fst values among populations.

108

APPENDIX

109

Table A.III.1. List of specimens used in this study, including collection information (catalog number, museum), locality of origin
(county and region), lineage, timing of population establishment, sex (0=males, 1= females), and coordinates of geographical midpoint
of the locality of origin. LP=Michigan’s Lower Peninsula, UP=Michigan’s Upper Peninsula, WI=Wisconsin.
Catalog number

Museum

County

ID

Region

Lineage

Timing

Sex

Lat

Log

37246

MSU

Adams

AD

WI

WUP-WI

Old

0

43.9012

-89.8676

37247

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37248

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37249

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37250

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37251

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37252

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37435

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37441

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37444

MSU

Adams

AD

WI

WUP-WI

Old

1

43.9012

-89.8676

37321

MSU

Ashland

AS

WI

WUP-WI

New

1

46.1568

-90.4755

37322

MSU

Ashland

AS

WI

WUP-WI

New

1

46.1568

-90.4755

37323

MSU

Ashland

AS

WI

WUP-WI

New

1

46.1568

-90.4755

37324

MSU

Ashland

AS

WI

WUP-WI

New

0

46.1568

-90.4755

110

Table A.III.1. Continued.
37325

MSU

Ashland

AS

WI

WUP-WI

New

1

46.1568

-90.4755

37327

MSU

Ashland

AS

WI

WUP-WI

New

1

46.1568

-90.4755

37328

MSU

Ashland

AS

WI

WUP-WI

New

0

46.1568

-90.4755

167367

UMMZ

Bayfield

BY

WI

WUP-WI

New

0

46.3705

-91.3376

167369

UMMZ

Bayfield

BY

WI

WUP-WI

New

1

46.3705

-91.3376

167372

UMMZ

Bayfield

BY

WI

WUP-WI

New

1

46.3705

-91.3376

167373

UMMZ

Bayfield

BY

WI

WUP-WI

New

1

46.3705

-91.3376

167375

UMMZ

Bayfield

BY

WI

WUP-WI

New

0

46.3705

-91.3376

167376

UMMZ

Bayfield

BY

WI

WUP-WI

New

0

46.3705

-91.3376

37319

MSU

Bayfield

BY

WI

WUP-WI

New

1

46.3705

-91.3376

37320

MSU

Bayfield

BY

WI

WUP-WI

New

0

46.3705

-91.3376

8204

UWSP

Bayfield

BY

WI

WUP-WI

New

1

46.3705

-91.3376

8205

UWSP

Bayfield

BY

WI

WUP-WI

New

1

46.3705

-91.3376

37255

MSU

Buffalo

BF

WI

WUP-WI

Old

1

44.5892

-91.8022

37260

MSU

Buffalo

BF

WI

WUP-WI

Old

1

44.5892

-91.8022

37262

MSU

Buffalo

BF

WI

WUP-WI

Old

0

44.5892

-91.8022

37263

MSU

Buffalo

BF

WI

WUP-WI

Old

0

44.5892

-91.8022

111

Table A.III.1. Continued.
37268

MSU

Buffalo

BF

WI

WUP-WI

Old

1

44.5892

-91.8022

37270

MSU

Buffalo

BF

WI

WUP-WI

Old

0

44.5892

-91.8022

37276

MSU

Buffalo

BF

WI

WUP-WI

Old

0

44.5892

-91.8022

37279

MSU

Buffalo

BF

WI

WUP-WI

Old

1

44.5892

-91.8022

37284

MSU

Buffalo

BF

WI

WUP-WI

Old

1

44.5892

-91.8022

37287

MSU

Buffalo

BF

WI

WUP-WI

Old

0

44.5892

-91.8022

37203

MSU

Charlevoix

CX

LP

LP-EUP

Old

1

45.2615

-84.8037

37204

MSU

Charlevoix

CX

LP

LP-EUP

Old

1

45.2615

-84.8037

37208

MSU

Charlevoix

CX

LP

LP-EUP

Old

1

45.2615

-84.8037

37211

MSU

Charlevoix

CX

LP

LP-EUP

Old

1

45.2615

-84.8037

37216

MSU

Charlevoix

CX

LP

LP-EUP

Old

0

45.2615

-84.8037

37220

MSU

Charlevoix

CX

LP

LP-EUP

Old

1

45.2615

-84.8037

37221

MSU

Charlevoix

CX

LP

LP-EUP

Old

0

45.2615

-84.8037

37223

MSU

Charlevoix

CX

LP

LP-EUP

Old

0

45.2615

-84.8037

37225

MSU

Charlevoix

CX

LP

LP-EUP

Old

0

45.2615

-84.8037

37228

MSU

Charlevoix

CX

LP

LP-EUP

Old

0

45.2615

-84.8037

175151

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

1

45.5634

-84.6637

112

Table A.III.1. Continued.
175214

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

1

45.5634

-84.6637

175216

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

1

45.5634

-84.6637

175217

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

1

45.5634

-84.6637

175218

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

0

45.5634

-84.6637

175220

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

1

45.5634

-84.6637

175221

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

0

45.5634

-84.6637

175226

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

1

45.5634

-84.6637

175228

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

0

45.5634

-84.6637

175230

UMMZ

Cheboygan

CH

LP

LP-EUP

Old

0

45.5634

-84.6637

177248

UMMZ

Chippewa

CW

UP

LP-EUP

New

1

46.4320

-84.7266

177249

UMMZ

Chippewa

CW

UP

LP-EUP

New

0

46.4320

-84.7266

177250

UMMZ

Chippewa

CW

UP

LP-EUP

New

0

46.4320

-84.7266

177253

UMMZ

Chippewa

CW

UP

LP-EUP

New

1

46.4320

-84.7266

177257

UMMZ

Chippewa

CW

UP

LP-EUP

New

1

46.4320

-84.7266

177259

UMMZ

Chippewa

CW

UP

LP-EUP

New

1

46.4320

-84.7266

177260

UMMZ

Chippewa

CW

UP

LP-EUP

New

0

46.4320

-84.7266

177261

UMMZ

Chippewa

CW

UP

LP-EUP

New

1

46.4320

-84.7266

113

Table A.III.1. Continued.
177267

UMMZ

Chippewa

CW

UP

LP-EUP

New

1

46.4320

-84.7266

177269

UMMZ

Chippewa

CW

UP

LP-EUP

New

0

46.4320

-84.7266

132046

UMMZ

Crawford

CR

LP

LP-EUP

Old

1

44.7888

-84.7340

132047

UMMZ

Crawford

CR

LP

LP-EUP

Old

0

44.7888

-84.7340

132049

UMMZ

Crawford

CR

LP

LP-EUP

Old

0

44.7888

-84.7340

132050

UMMZ

Crawford

CR

LP

LP-EUP

Old

1

44.7888

-84.7340

156424

UMMZ

Crawford

CR

LP

LP-EUP

Old

1

44.7888

-84.7340

156425

UMMZ

Crawford

CR

LP

LP-EUP

Old

0

44.7888

-84.7340

167287

UMMZ

Crawford

CR

LP

LP-EUP

Old

0

44.7888

-84.7340

167329

UMMZ

Crawford

CR

LP

LP-EUP

Old

1

44.7888

-84.7340

167330

UMMZ

Crawford

CR

LP

LP-EUP

Old

1

44.7888

-84.7340

167331

UMMZ

Crawford

CR

LP

LP-EUP

Old

1

44.7888

-84.7340

175993

UMMZ

Delta

DE

UP

WUP-WI

New

0

46.0788

-86.8345

175996

UMMZ

Delta

DE

UP

WUP-WI

New

0

46.0788

-86.8345

176004

UMMZ

Delta

DE

UP

WUP-WI

New

1

46.0788

-86.8345

176005

UMMZ

Delta

DE

UP

WUP-WI

New

0

46.0788

-86.8345

176006

UMMZ

Delta

DE

UP

WUP-WI

New

0

46.0788

-86.8345

114

Table A.III.1. Continued.
176007

UMMZ

Delta

DE

UP

WUP-WI

New

0

46.0788

-86.8345

176014

UMMZ

Delta

DE

UP

WUP-WI

New

1

46.0788

-86.8345

176015

UMMZ

Delta

DE

UP

WUP-WI

New

1

46.0788

-86.8345

176016

UMMZ

Delta

DE

UP

WUP-WI

New

1

46.0788

-86.8345

176046

UMMZ

Delta

DE

UP

WUP-WI

New

1

46.0788

-86.8345

165217

UMMZ

Emmet

EM

LP

LP-EUP

Old

1

45.5431

-85.0478

165265

UMMZ

Emmet

EM

LP

LP-EUP

Old

1

45.5431

-85.0478

165267

UMMZ

Emmet

EM

LP

LP-EUP

Old

0

45.5431

-85.0478

165270

UMMZ

Emmet

EM

LP

LP-EUP

Old

1

45.5431

-85.0478

165277

UMMZ

Emmet

EM

LP

LP-EUP

Old

1

45.5431

-85.0478

165278

UMMZ

Emmet

EM

LP

LP-EUP

Old

0

45.5431

-85.0478

165285

UMMZ

Emmet

EM

LP

LP-EUP

Old

1

45.5431

-85.0478

167388

UMMZ

Emmet

EM

LP

LP-EUP

Old

0

45.5431

-85.0478

167390

UMMZ

Emmet

EM

LP

LP-EUP

Old

0

45.5431

-85.0478

167394

UMMZ

Emmet

EM

LP

LP-EUP

Old

0

45.5431

-85.0478

177159

UMMZ

Gogebic

GO

UP

WUP-WI

New

1

46.3437

-89.3678

177175

UMMZ

Gogebic

GO

UP

WUP-WI

New

0

46.3437

-89.3678

115

Table A.III.1. Continued.
177388

UMMZ

Gogebic

GO

UP

WUP-WI

New

1

46.3437

-89.3678

177389

UMMZ

Gogebic

GO

UP

WUP-WI

New

0

46.3437

-89.3678

177390

UMMZ

Gogebic

GO

UP

WUP-WI

New

1

46.3437

-89.3678

177391

UMMZ

Gogebic

GO

UP

WUP-WI

New

1

46.3437

-89.3678

177392

UMMZ

Gogebic

GO

UP

WUP-WI

New

0

46.3437

-89.3678

177393

UMMZ

Gogebic

GO

UP

WUP-WI

New

0

46.3437

-89.3678

37331

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37333

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37334

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37336

MSU

Iron

IR

WI

WUP-WI

New

0

46.3579

-90.3774

37337

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37339

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37342

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37343

MSU

Iron

IR

WI

WUP-WI

New

0

46.3579

-90.3774

37344

MSU

Iron

IR

WI

WUP-WI

New

1

46.3579

-90.3774

37346

MSU

Iron

IR

WI

WUP-WI

New

0

46.3579

-90.3774

37373

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

116

Table A.III.1. Continued.
37375

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

37376

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

37377

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

37378

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

37380

MSU

Marathon

MH

WI

WUP-WI

Old

1

44.8175

-89.7509

37381

MSU

Marathon

MH

WI

WUP-WI

Old

1

44.8175

-89.7509

37382

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

37384

MSU

Marathon

MH

WI

WUP-WI

Old

0

44.8175

-89.7509

37391

MSU

Marathon

MH

WI

WUP-WI

Old

1

44.8175

-89.7509

37400

MSU

Marinette

MR

WI

WUP-WI

Old

0

45.5900

-87.9464

37401

MSU

Marinette

MR

WI

WUP-WI

Old

1

45.5900

-87.9464

37402

MSU

Marinette

MR

WI

WUP-WI

Old

1

45.5900

-87.9464

37406

MSU

Marinette

MR

WI

WUP-WI

Old

0

45.5900

-87.9464

37407

MSU

Marinette

MR

WI

WUP-WI

Old

0

45.5900

-87.9464

37409

MSU

Marinette

MR

WI

WUP-WI

Old

1

45.5900

-87.9464

37412

MSU

Marinette

MR

WI

WUP-WI

Old

0

45.5900

-87.9464

37417

MSU

Marinette

MR

WI

WUP-WI

Old

0

45.5900

-87.9464

117

Table A.III.1. Continued.
37427

MSU

Marinette

MR

WI

WUP-WI

Old

1

45.5900

-87.9464

37430

MSU

Marinette

MR

WI

WUP-WI

Old

1

45.5900

-87.9464

177394

UMMZ

Marquette

MA

UP

WUP-WI

New

1

46.8736

-87.8958

177395

UMMZ

Marquette

MA

UP

WUP-WI

New

1

46.8736

-87.8958

177396

UMMZ

Marquette

MA

UP

WUP-WI

New

1

46.8736

-87.8958

177397

UMMZ

Marquette

MA

UP

WUP-WI

New

0

46.8736

-87.8958

177398

UMMZ

Marquette

MA

UP

WUP-WI

New

0

46.8736

-87.8958

177399

UMMZ

Marquette

MA

UP

WUP-WI

New

0

46.8736

-87.8958

177400

UMMZ

Marquette

MA

UP

WUP-WI

New

1

46.8736

-87.8958

177401

UMMZ

Marquette

MA

UP

WUP-WI

New

0

46.8736

-87.8958

175318

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

176017

UMMZ

Menominee

MN

UP

WUP-WI

Old

0

45.1580

-87.7000

176022

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

176026

UMMZ

Menominee

MN

UP

WUP-WI

Old

0

45.1580

-87.7000

176027

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

176029

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

176047

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

118

Table A.III.1. Continued.
176181

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

176182

UMMZ

Menominee

MN

UP

WUP-WI

Old

1

45.1580

-87.7000

176191

UMMZ

Menominee

MN

UP

WUP-WI

Old

0

45.1580

-87.7000

37310

MSU

Oneida

OD

WI

WUP-WI

New

0

45.8709

-89.6527

37311

MSU

Oneida

OD

WI

WUP-WI

New

1

45.8709

-89.6527

37312

MSU

Oneida

OD

WI

WUP-WI

New

1

45.8709

-89.6527

37313

MSU

Oneida

OD

WI

WUP-WI

New

0

45.8709

-89.6527

37314

MSU

Oneida

OD

WI

WUP-WI

New

1

45.8709

-89.6527

37315

MSU

Oneida

OD

WI

WUP-WI

New

0

45.8709

-89.6527

37316

MSU

Oneida

OD

WI

WUP-WI

New

0

45.8709

-89.6527

37317

MSU

Oneida

OD

WI

WUP-WI

New

1

45.8709

-89.6527

37318

MSU

Oneida

OD

WI

WUP-WI

New

1

45.8709

-89.6527

177181

UMMZ

Ontonagon

ON

UP

WUP-WI

New

0

46.7364

-89.7394

177184

UMMZ

Ontonagon

ON

UP

WUP-WI

New

0

46.7364

-89.7394

177186

UMMZ

Ontonagon

ON

UP

WUP-WI

New

0

46.7364

-89.7394

177194

UMMZ

Ontonagon

ON

UP

WUP-WI

New

1

46.7364

-89.7394

177201

UMMZ

Ontonagon

ON

UP

WUP-WI

New

1

46.7364

-89.7394

119

Table A.III.1. Continued.
177211

UMMZ

Ontonagon

ON

UP

WUP-WI

New

0

46.7364

-89.7394

177212

UMMZ

Ontonagon

ON

UP

WUP-WI

New

0

46.7364

-89.7394

177214

UMMZ

Ontonagon

ON

UP

WUP-WI

New

1

46.7364

-89.7394

177216

UMMZ

Ontonagon

ON

UP

WUP-WI

New

1

46.7364

-89.7394

277222

UMMZ

Ontonagon

ON

UP

WUP-WI

New

1

46.7364

-89.7394

170135

UMMZ

Otsego

OT

LP

LP-EUP

Old

0

45.0549

-84.5688

170137

UMMZ

Otsego

OT

LP

LP-EUP

Old

0

45.0549

-84.5688

170316

UMMZ

Otsego

OT

LP

LP-EUP

Old

1

45.0549

-84.5688

170323

UMMZ

Otsego

OT

LP

LP-EUP

Old

1

45.0549

-84.5688

170454

UMMZ

Otsego

OT

LP

LP-EUP

Old

0

45.0549

-84.5688

175470

UMMZ

Otsego

OT

LP

LP-EUP

Old

1

45.0549

-84.5688

175471

UMMZ

Otsego

OT

LP

LP-EUP

Old

1

45.0549

-84.5688

175472

UMMZ

Otsego

OT

LP

LP-EUP

Old

0

45.0549

-84.5688

175473

UMMZ

Otsego

OT

LP

LP-EUP

Old

0

45.0549

-84.5688

175474

UMMZ

Otsego

OT

LP

LP-EUP

Old

1

45.0549

-84.5688

37230

MSU

Outagamie

OU

WI

WUP-WI

Old

1

44.3961

-88.6663

37233

MSU

Outagamie

OU

WI

WUP-WI

Old

1

44.3961

-88.6663

120

Table A.III.1. Continued.
37234

MSU

Outagamie

OU

WI

WUP-WI

Old

0

44.3961

-88.6663

37235

MSU

Outagamie

OU

WI

WUP-WI

Old

1

44.3961

-88.6663

37236

MSU

Outagamie

OU

WI

WUP-WI

Old

0

44.3961

-88.6663

37238

MSU

Outagamie

OU

WI

WUP-WI

Old

1

44.3961

-88.6663

37239

MSU

Outagamie

OU

WI

WUP-WI

Old

1

44.3961

-88.6663

37240

MSU

Outagamie

OU

WI

WUP-WI

Old

0

44.3961

-88.6663

37244

MSU

Outagamie

OU

WI

WUP-WI

Old

1

44.3961

-88.6663

37245

MSU

Outagamie

OU

WI

WUP-WI

Old

0

44.3961

-88.6663

8162

UWSP

Portage

PO

WI

WUP-WI

Old

0

44.5412

-89.5654

5385

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

5467

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

5468

UWSP

Portage

PO

WI

WUP-WI

Old

0

44.5412

-89.5654

5469

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

5535

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

5539

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

5671

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

6055

UWSP

Portage

PO

WI

WUP-WI

Old

1

44.5412

-89.5654

121

Table A.III.1. Continued.
8163

UWSP

Portage

PO

WI

WUP-WI

Old

0

44.5412

-89.5654

175526

UMMZ

Schoolcraft

SC

UP

WUP-WI

New

0

46.0886

-86.2278

175527

UMMZ

Schoolcraft

SC

UP

WUP-WI

New

1

46.0886

-86.2278

175528

UMMZ

Schoolcraft

SC

UP

WUP-WI

New

0

46.0886

-86.2278

175541

UMMZ

Schoolcraft

SC

UP

WUP-WI

New

0

46.0886

-86.2278

175542

UMMZ

Schoolcraft

SC

UP

WUP-WI

New

0

46.0886

-86.2278

37289

MSU

Taylor

TA

WI

WUP-WI

Old

1

45.3127

-90.6268

37291

MSU

Taylor

TA

WI

WUP-WI

Old

1

45.3127

-90.6268

37294

MSU

Taylor

TA

WI

WUP-WI

Old

1

45.3127

-90.6268

37295

MSU

Taylor

TA

WI

WUP-WI

Old

0

45.3127

-90.6268

37296

MSU

Taylor

TA

WI

WUP-WI

Old

0

45.3127

-90.6268

37298

MSU

Taylor

TA

WI

WUP-WI

Old

1

45.3127

-90.6268

37299

MSU

Taylor

TA

WI

WUP-WI

Old

1

45.3127

-90.6268

37304

MSU

Taylor

TA

WI

WUP-WI

Old

0

45.3127

-90.6268

37305

MSU

Taylor

TA

WI

WUP-WI

Old

0

45.3127

-90.6268

37308

MSU

Taylor

TA

WI

WUP-WI

Old

1

45.3127

-90.6268

37347

MSU

Wood

WO

WI

WUP-WI

Old

0

44.3785

-90.1090

122

Table A.III.1. Continued.
37348

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

37349

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

37351

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

37352

MSU

Wood

WO

WI

WUP-WI

Old

0

44.3785

-90.1090

37353

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

37356

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

37359

MSU

Wood

WO

WI

WUP-WI

Old

0

44.3785

-90.1090

37361

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

37362

MSU

Wood

WO

WI

WUP-WI

Old

1

44.3785

-90.1090

123

Table A.III.2. Matrix of pairwise genetic distances (Nei’s D) between populations included in this study. Above diagonal: Average
number of pairwise differences between populations (PiXY). Diagonal elements (in bold): Average number of pairwise differences
within population (PiX). Below diagonal: Corrected average pairwise difference (PiXY-(PiX+PiY)/2).

AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

AD
7.1
0.2
0.4
0.4
19.0
19.4
10.3
17.8
1.0
20.8
0.7
0.6
0.8
0.4
0.6
2.1
0.2
0.7
19.1
1.0
1.4
1.9
0.5
0.8

AS
7.2
7.0
-0.2
0.3
19.0
19.6
10.4
17.7
0.9
20.8
0.2
0.2
0.5
0.4
0.4
2.5
0.1
0.4
19.1
1.0
0.7
2.2
0.2
0.7

BY
7.0
6.3
6.1
0.6
19.3
19.8
10.6
18.1
1.0
21.2
0.8
0.2
1.0
0.6
0.8
2.4
0.8
0.8
19.5
1.0
1.2
2.1
0.4
1.1

BF
7.5
7.4
7.2
7.1
17.7
18.0
9.1
16.4
0.4
19.4
0.4
0.6
0.3
0.1
0.9
1.6
0.7
0.5
17.8
0.2
0.9
1.2
0.3
0.4

CX
26.3
26.1
26.1
24.9
7.4
1.2
5.0
-0.1
18.5
-0.2
19.4
18.5
19.4
18.6
19.0
18.8
19.6
15.6
0.3
19.1
17.7
18.9
18.0
19.9

CH
28.0
28.1
27.9
26.6
9.9
10.0
4.3
0.9
18.6
1.7
20.0
18.9
19.8
19.1
19.4
18.7
20.1
15.6
0.1
19.3
18.2
18.8
18.6
20.5

124

CW
22.2
22.3
22.1
21.0
17.1
17.7
16.8
4.4
9.6
5.4
10.8
10.2
10.6
10.0
10.3
9.9
10.9
7.9
5.1
10.2
9.6
10.0
9.5
11.2

CR
25.3
25.1
25.1
23.9
7.5
9.9
16.7
7.9
17.1
0.5
18.1
17.3
18.1
17.3
17.8
17.3
18.4
14.4
0.0
17.7
16.4
17.4
16.8
18.6

DE
6.4
6.2
5.9
5.8
24.1
25.4
19.8
22.9
3.7
20.5
1.0
0.9
1.1
0.5
1.2
1.4
1.3
0.7
18.3
0.4
1.6
0.9
0.8
1.6

EM
27.9
27.8
27.7
26.4
7.0
10.2
17.3
7.9
25.8
6.9
21.2
20.4
21.1
20.5
20.5
20.6
21.4
17.2
1.1
21.1
19.5
20.7
19.6
21.6

GO
7.2
6.7
6.9
7.0
26.2
28.0
22.2
25.1
5.9
27.7
6.1
0.6
-0.1
0.6
0.7
2.7
0.6
0.3
19.5
0.9
0.8
2.3
0.4
0.5

IR
6.4
6.0
5.6
6.4
24.5
26.2
20.8
23.5
5.0
26.1
5.9
4.5
0.9
0.9
0.8
3.0
0.9
0.5
18.6
1.1
1.7
2.5
0.6
0.9

Table A.III.2. Continued.

AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

MH
8.1
7.6
7.8
7.6
26.8
28.5
22.6
25.7
6.6
28.2
6.7
6.8
7.4
0.6
1.2
2.7
0.8
0.6
19.4
0.9
0.9
2.2
0.7
0.3

MR
7.0
6.9
6.8
6.7
25.4
27.2
21.4
24.3
5.4
27.0
6.7
6.2
7.4
6.2
0.9
0.8
0.5
0.9
18.6
0.2
1.0
0.7
0.4
0.6

MA
6.8
6.6
6.5
7.1
25.4
27.1
21.4
24.4
5.7
26.7
6.5
5.7
7.6
6.7
5.4
3.1
0.5
0.6
19.2
1.5
2.3
2.8
0.1
1.1

MN
7.9
8.2
7.7
7.4
24.7
26.0
20.5
23.4
5.5
26.3
8.0
7.5
8.6
6.2
8.1
4.5
2.9
2.8
18.4
1.5
2.0
-0.1
2.3
3.1

OD
6.6
6.4
6.7
7.1
26.2
28.0
22.1
25.2
6.0
27.7
6.4
6.0
7.4
6.4
6.1
8.0
5.7
0.6
19.7
1.4
1.6
2.6
0.8
0.9

ON
9.3
8.9
8.9
9.1
24.3
25.7
21.3
23.4
7.6
25.8
8.4
7.8
9.3
9.0
8.3
10.1
8.5
10.1
15.4
1.2
1.3
2.4
0.6
0.9

125

OT
26.9
26.9
26.8
25.6
8.3
9.3
17.7
8.2
24.4
8.8
26.8
25.1
27.3
26.0
26.1
24.9
26.8
24.8
8.5
19.0
17.8
18.5
18.2
20.0

OU
7.6
7.6
7.2
6.8
25.9
27.4
21.6
24.7
5.4
27.7
7.1
6.4
7.6
6.4
7.3
6.8
7.4
9.3
26.3
6.2
1.5
1.0
0.6
0.9

PO
8.5
7.8
7.9
8.0
25.0
26.8
21.5
23.9
7.0
26.6
7.5
7.6
8.1
7.7
8.5
7.8
8.0
10.0
25.6
8.1
7.2
1.9
1.2
1.7

SC
7.5
7.8
7.3
6.9
24.7
25.9
20.4
23.4
4.9
26.3
7.4
6.9
8.0
5.9
7.6
4.3
7.6
9.6
24.9
6.2
7.6
4.2
2.0
2.8

TA
7.7
7.4
7.1
7.5
25.4
27.3
21.6
24.4
6.4
26.7
7.1
6.5
8.0
7.2
6.5
8.3
7.3
9.3
26.1
7.4
8.5
7.8
7.4
0.7

WO
7.2
7.0
7.0
6.8
26.4
28.3
22.4
25.4
6.3
27.9
6.3
5.9
6.8
6.5
6.6
8.2
6.6
8.8
27.1
6.9
8.1
7.7
7.2
5.6

Table A.III.3. Matrix of pairwise Fst values between populations included in this study.

AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

AD
0.0
0.0
0.1
0.1
0.7
0.7
0.5
0.7
0.2
0.7
0.1
0.1
0.1
0.1
0.1
0.3
0.0
0.1
0.7
0.1
0.2
0.2
0.1
0.1

AS

BY

BF

CX

CH

CW

CR

DE

EM

GO

IR

0.0
0.0
0.0
0.7
0.7
0.4
0.7
0.2
0.7
0.0
0.0
0.1
0.1
0.1
0.3
0.0
0.0
0.7
0.1
0.1
0.3
0.0
0.1

0.0
0.1
0.7
0.7
0.5
0.7
0.2
0.8
0.1
0.0
0.1
0.1
0.1
0.3
0.1
0.1
0.7
0.1
0.2
0.3
0.1
0.2

0.0
0.7
0.7
0.4
0.7
0.1
0.7
0.1
0.1
0.0
0.0
0.1
0.2
0.1
0.0
0.7
0.0
0.1
0.2
0.0
0.1

0.0
0.1
0.3
0.0
0.8
0.0
0.7
0.7
0.7
0.7
0.8
0.8
0.7
0.6
0.0
0.7
0.7
0.8
0.7
0.8

0.0
0.2
0.1
0.7
0.1
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.6
0.0
0.7
0.7
0.7
0.7
0.7

0.0
0.3
0.5
0.3
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.4
0.3
0.5
0.4
0.5
0.4
0.5

0.0
0.7
0.1
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.6
0.0
0.7
0.7
0.7
0.7
0.7

0.0
0.8
0.2
0.2
0.2
0.1
0.2
0.3
0.2
0.1
0.8
0.1
0.2
0.2
0.2
0.3

0.0
0.8
0.8
0.7
0.8
0.8
0.8
0.8
0.6
0.1
0.8
0.7
0.8
0.7
0.8

0.0
0.1
0.0
0.1
0.1
0.3
0.1
0.0
0.7
0.1
0.1
0.3
0.1
0.1

0.0
0.1
0.1
0.1
0.4
0.2
0.0
0.7
0.2
0.2
0.4
0.1
0.1

126

Table A.III.3. Continued.
MH
AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

MR

MA

MN

OD

ON

OT

OU

PO

SC

TA

WO

0.0
0.1
0.2
0.3
0.1
0.1
0.7
0.1
0.1
0.3
0.1
0.0

0.0
0.1
0.1
0.1
0.1
0.7
0.0
0.1
0.1
0.1
0.1

0.0
0.4
0.1
0.1
0.7
0.2
0.3
0.4
0.0
0.2

0.0
0.4
0.3
0.7
0.2
0.3
0.0
0.3
0.4

0.0
0.1
0.7
0.2
0.2
0.4
0.1
0.1

0.0
0.6
0.1
0.1
0.2
0.1
0.1

0.0
0.7
0.7
0.7
0.7
0.7

0.0
0.2
0.2
0.1
0.1

0.0
0.2
0.1
0.2

0.0
0.3
0.4

0.0
0.1

0.0

127

Table A.III.4. Matrix of pairwise morphometric distances between populations included in this study. Distances were calculated as
pairwise Euclidean distances between population means computed from the shape variables derived from projecting Procrustes
residuals onto the first 46 eigenvectors of their covariance matrix. Values in each cell were multiply by 103.

AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

AD
0
13
14
9
16
20
14
19
14
17
16
11
11
15
24
13
10
14
15
12
14
18
12
11

AS

BY

BF

CX

CH

CW

CR

DE

EM

GO

IR

0
16
14
16
17
21
21
21
14
13
8
15
20
19
19
14
15
15
14
15
20
13
16

0
16
18
22
13
20
17
18
16
15
14
16
23
18
14
18
17
16
17
16
15
13

0
17
18
15
16
14
15
15
11
13
15
25
13
11
12
17
16
17
22
17
14

0
17
20
22
23
18
14
14
14
16
17
18
17
18
13
12
11
20
15
14

0
23
16
23
13
16
15
20
20
16
19
18
15
15
18
16
26
21
20

0
17
12
22
17
18
16
14
28
14
14
17
18
18
20
20
20
13

0
15
17
20
18
20
16
23
14
17
13
16
21
21
26
24
19

0
21
21
20
17
14
28
11
15
16
19
21
23
22
21
15

0
17
14
20
20
19
18
18
14
19
20
20
25
19
20

0
11
14
18
20
17
12
15
13
15
13
20
16
15

0
12
17
19
17
11
14
13
11
12
20
14
13

128

Table A.III.4. Continued.
MH
AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

MR

MA

MN

OD

ON

OT

OU

PO

SC

TA

WO

0
16
22
17
10
17
14
11
11
16
12
8

0
22
11
15
16
14
17
17
23
18
13

0
23
23
21
18
20
17
27
21
21

0
14
12
15
18
18
24
19
14

0
13
14
13
13
19
14
11

0
16
18
18
25
20
16

0
11
9
18
15
14

0
8
17
11
12

0
18
13
13

0
14
17

0
13

0

129

BIBLIOGRAPHY

130

BIBLIOGRAPHY

Anapol F, Lee, S (1994) Morphological adaptation to diet in Platyrrhine primates. American
Journal of Physical Anthropology, 94, 239-261.

Anderson TW (1963) Asymptotic Theory for Principal Component Analysis. Annals of
Mathematical Statistics, 34, 122-148.

Angers B, Castonguay E, Massicotte R (2010) Environmentally induced phenotypes and DNA
methylation: how to deal with unpredictable conditions until the next generation and
after. Molecular Ecology, 19, 1283-1295.

Atchley W R (1993) Genetic and Developmental Aspects of Variability in the Mammalian
Mandible. In: The Skull (eds. Hanken J, Hall, BK), pp. 207-247. University of Chicago
Press, Chicago, Illinois.

Atchley W R, Cowley DE, Vogl C, McLellan T (1992) Evolutionary divergence, shape change,
and genetic correlation structure in the rodent mandible. Systematic Biology , 41, 196221.

Atchley W R, Hall BK (1991) A model for development and evolution of complex
morphological structures. Biological Reviews of the Cambridge Philosophical Society 66,
101-157.

Atchley WR, Plummer AA, and Riska B (1985) Genetics of mandible form in the mouse.
Genetics, 111, 555-577.

Badyaev AV (2009) Evolutionary significance of phenotypic accommodation in novel
environments: an empirical test of the Baldwin effect. Philosophical Transactions of the
Royal Society B-Biological Sciences, 364, 1125-1141.

Badyaev AV (2010) The beak of the other finch: coevolution of genetic covariance structure and
developmental modularity during adaptive evolution. Philosophical Transactions of the
Royal Society B-Biological Sciences, 365, 1111-1126.

131

Badyaev AV, Foresman, KR (2000) Extreme environmental change and evolution: stressinduced morphological variation is strongly concordant with patterns of evolutionary
divergence in shrew mandibles. Proceedings of the Royal Society B-Biological Sciences,
267, 371-377.

Bookstein F L (1991) Morphometric tools for landmark data: geometry and biology. Cambridge
University Press, Cambridge [England]; New York.

Bookstein FL (1997) Landmark methods for forms without landmarks: morphometrics of group
differences in outline shape. Medical Image Analysis, 1, 225-243.

Bouvier M, Hylander WL (1981) The effect of dietary consistency on the morphology of the
mandibular condylar cartilage in macaques. American Journal of Physical Anthropology,
54, 203-204.

Bouvier M, Hylander WL (1984) The effect of dietary consistency on gross and histologic
morphology in the craniofacial region of young rats. American Journal of Anatomy, 170,
117-126.

Bresin A, Kiliaridis S (2002) Dento-skeletal adaptation after bite-raising in growing rats with
different masticatory muscle capacities. European Journal of Orthodontics, 24, 223-237.

Canale C I, Henry PY (2010) Adaptive phenotypic plasticity and resilience of vertebrates to
increasing climatic unpredictability. Climate Research, 43, 135-147.

Caumul R, Polly PD (2005) Phylogenetic and environmental components of morphological
variation: Skull, mandible, and molar shape in marmots (Marmota, Rodentia). Evolution,
59, 2460-2472.

Cheverud JM, Hartman SE, Richtsmeier JT, Atchley WR (1991) A Quantitative genetic analysis
of localized morphology in mandibles of inbred mice using finite element scaling
analysis. Journal of Craniofacial Genetics and Developmental Biology, 11, 122-137.

Corruccini RS, Beecher RM (1982) Occlusal variation related to soft diet in a non-human
primate. Science, 218, 74-76.

132

Crispo E (2007) The Baldwin effect and genetic assimilation: Revisiting two mechanisms of
evolutionary change mediated by phenotypic plasticity. Evolution, 61, 2469-2479.

de Jong G, Crozier RH (2003) Developmental plasticity and evolution. Nature, 424, 16-17.
Dietz EJ (1983) Permutation tests for association between 2 distance matrices. Systematic
Zoology, 32, 21-26.

Dixon AD, Hoyte DAN, Rönning O (1997) Fundamentals of craniofacial growth. CRC Press,
Boca Raton, Florida.

Dryden IL, Mardia KV (1998) Statistical shape analysis. Wiley, New York.

Duarte LC, Monteiro LR, Von Zuben FJ, Dos Reis SF (2000) Variation in mandible shape in
Thrichomys apereoides (Mammalia : Rodentia): Geometric analysis of a complex
morphological structure. Systematic Biology, 49, 563-578.

Elrod DA, Kennedy ML (1995) Microgeographic variation in morphometric characters of the
white-footed mouse, Peromyscus leucopus. Southwestern Naturalist, 40, 42-49.

Excoffier L, Laval G, Schneider S (2005) Arlequin (version 3.0): An integrated software package
for population genetics data analysis. Evolutionary Bioinformatics, 1, 47-50.

Ghafari J, Heeley JD (1982) Condylar adaptation to muscle alteration in the rat. Angle
Orthodontist, 52, 26-37.

Goheen JR, Swihart RK, Robins JH (2003) The anatomy of a range expansion: changes in
cranial morphology and rates of energy extraction for North American red squirrels from
different latitudes. Oikos, 102, 33-44.

Hall BK, Miyake T (2000) All for one and one for all: condensations and the initiation of skeletal
development. Bioessays, 22, 138-147.

Hanken J, Hall BK (1993) The Skull. University of Chicago Press, Chicago, Illinois.

133

Herring SW (1993) Epigenetic and Functional Influences on Skull Growth. In: The Skull (eds.
Hanken J, Hall, BK), pp. 153-206. University of Chicago Press, Chicago, Illinois.

Holbrook SJ (1982) Ecological inferences from mandibular morphology of Peromyscus
maniculatus. Journal of Mammalogy, 63, 399-408.

Kantomaa T, Rönning O (1997) Growth Mechanisms of the Mandible. In: Craniofacial Growth
(eds. Dixon AD, Hoyte DAN, Rönning O), pp 189-204. CRC Press, Boca Raton.Florida.

Katsnelson A (2010) Epigenome effort makes its mark. Nature, 467, 646.

Kendall DG (1977) Diffusion of shape. Advances in Applied Probability, 9, 428-430.

Kiliaridis S, Thilander B, Kjellberg H, Topouzelis N, Zafiriadis A (1999) Effect of low
masticatory function on condylar growth: A morphometric study in the rat. American
Journal of Orthodontics and Dentofacial Orthopedics, 116, 121-125.

Klingenberg CP, Badyaev AV, Sowry SM, Beckwith NJ (2001a). Inferring developmental
modularity from morphological integration: Analysis of individual variation and
asymmetry in bumblebee wings. American Naturalist, 157, 11-23.

Klingenberg CP, Leamy LJ, Routman EJ, Cheverud JM (2001b) Genetic architecture of
mandible shape in mice: Effects of quantitative trait loci analyzed by geometric
morphometrics. Genetics, 157, 785-802.

Klingenberg CP, Mebus K, Auffray JC (2003) Developmental integration in a complex
morphological structure: how distinct are the modules in the mouse mandible? Evolution
& Development, 5, 522-531.

Krzanowski WJ (2000) Principles of multivariate analysis: a user's perspective. Oxford
University Press, Oxford.

Lande R (2009) Adaptation to an extraordinary environment by evolution of phenotypic
plasticity and genetic assimilation. Journal of Evolutionary Biology, 22, 1435-1446.

134

LeBoulenge E, Legendre P, delaCourt C, LeBoulengeNguyen P, Languy M (1996)
Microgeographic morphological differentiation in muskrats. Journal of Mammalogy, 77,
684-701.

Liu ZJ, Ikeda K, Harada S, Kasahara Y, Ito G (1998) Functional properties of jaw and tongue
muscles in rats fed a liquid diet after being weaned. Journal of Dental Research, 77, 366376.

Maki K, Nishioka T, Shioiri E, Takahashi T, Kimura M (2002) Effects of dietary consistency on
the mandible of rats at the growth stage: Computed X-ray densitometric and
cephalometric analysis. Angle Orthodontist, 72, 468-475.

Manly BF (2007) Randomization, bootstrap and Monte Carlo methods in biology. CRC Press.,
Boca Raton, Florida.

Marquez, E. J. 2006. SemiThinner. University of Michigan, Ann Arbor, Michigan.

Michaux J, Chevret P, Renaud S (2007) Morphological diversity of Old World rats and mice
(Rodentia, Muridae) mandible in relation with phylogeny and adaptation. Journal of
Zoological Systematics and Evolutionary Research, 45, 263-279.

Millien V (2006) Morphological evolution is accelerated among island mammals. Plos Biology,
4, 1863-1868.

Monteiro, LR, Bonato V, dos Reis SF (2005) Evolutionary integration and morphological
diversification in complex morphological structures: mandible shape divergence in spiny
rats (Rodentia, Echimyidae). Evolution & Development, 7, 429-439.

Monteiro LR, Nogueira MR (2010) Adaptive radiations, ecological specialization, and the
evolutionary integration of complex morphological structures. Evolution, 64, 724-743.

Myers P, Lundrigan BL, Gillespie BW, Zelditch ML (1996) Phenotypic plasticity in skull and
dental morphology in the prairie deer mouse (Peromyscus maniculatus bairdii). Journal
of Morphology, 229, 229-237.

135

Myers P, Lundrigan BL, Hoffman SMG, Haraminac AP, Seto SH (2009) Climate-induced
changes in the small mammal communities of the Northern Great Lakes Region. Global
Change Biology, 15, 1434-1454.

Nei M (1973) Analysis of gene diversity in subdivided populations. Proceedings of the National
Academy of Sciences of the United States of America, 70, 3321-3323.

Nei M (1987) Molecular evolutionary genetics. Columbia University Press, New York.

Nogueira M R, Monteiro LR, Peracchi AL, de Araujo AFB (2005) Ecomorphological analysis of
the masticatory apparatus in the seed-eating bats, genus Chiroderma (Chiroptera :
Phyllostomidae). Journal of Zoology 266, 355-364.

Nogueira MR, Peracchi AL, Monteiro LR (2009) Morphological correlates of bite force and diet
in the skull and mandible of phyllostomid bats. Functional Ecology, 23, 715-723.

Perez SI, Diniz JAF, Rohlf FJ, Dos Reis SD (2009) Ecological and evolutionary factors in the
morphological diversification of South American spiny rats. Biological Journal of the
Linnean Society, 98, 646-660.

Pergams ORW, Ashley MV (1999) Rapid morphological change in channel island deer mice.
Evolution, 53, 1573-1581.

Pergams ORW, Lacy RC (2008) Rapid morphological and genetic change in Chicago-area
Peromyscus. Molecular Ecology, 17, 450-463.

Pergams ORW, Lawler JJ (2009) Recent and Widespread Rapid Morphological Change in
Rodents. Plos One, 4, e6452

Pfennig DW, Wund MA, Snell-Rood EC, Cruickshank T, Schlichting CD, Moczek AP (2010)
Phenotypic plasticity's impacts on diversification and speciation. Trends in Ecology &
Evolution, 25, 459-467.

Renaud S (2005) First upper molar and mandible shape of wood mice (Apodemus sylvaticus)
from northern Germany: ageing, habitat and insularity. Mammalian Biology, 70, 157-170.

136

Renaud S, Auffray JC (2010) Adaptation and plasticity in insular evolution of the house mouse
mandible. Journal of Zoological Systematics and Evolutionary Research, 48, 138-150.

Renaud S, Auffray JC, de la Porte S (2010) Epigenetic effects on the mouse mandible: common
features and discrepancies in remodeling due to muscular dystrophy and response to food
consistency. BMC Evolutionary Biology, 10, 28.

Renaud S, Michaux J (2004) Parallel evolution in molar outline of murine rodents: the case of
the extinct Malpaisomys insularis (Eastern Canary Islands). Zoological Journal of the
Linnean Society, 142, 555-572.

Renaud S, Michaux JR (2007) Mandibles and molars of the wood mouse, Apodemus sylvaticus
(L.): integrated latitudinal pattern and mosaic insular evolution. Journal of Biogeography,
34, 339-355.

Rencher A (2002) Methods of multivariate analysis. Wiley, New York.

Reynolds J, Weir BS, Cockerham CC (1983) Estimation of the co-ancestry coefficient - basis for
a short-term genetic-distance. Genetics, 105, 767-779.

Rinker GC (1954) The comparative myology of the mammalian genera Sigmodon, Oryzomys,
Neotoma, and Peromyscus (Cricetinae) with remarks on their intergeneric relationships.
University of Michigan Press, University of Michigan. Ann Arbor, Michigan.

Rohlf FJ (2006) tpsdig2.1. State University of New York, Stony Brook, New York.

Rohlf FJ, Slice D (1990) Extensions of the procrustes method for the optimal superimposition of
landmarks. Systematic Zoology, 39, 40-59.

Samuels JX (2009) Cranial morphology and dietary habits of rodents. Zoological Journal of the
Linnean Society, 156, 864-888.

Schlichting C, Pigliucci M (1998) Phenotypic evolution: a reaction norm perspective. Sinauer,
Sunderland, Massachussets.

137

Schlichting CD (2004) The role of phenotypic plasticity in diversification. In: (eds. DeWitt TJ,
Scheiner SM). Phenotypic Plasticity. Functional and Conceptual Approaches, pp.191200. Oxford University Press, New York.

Sheets HD (2009) SemiLand6. Canisius College Buffalo, New York.

Sokal RR, Rohlf FJ (1995) Biometry: the principles and practice of statistics in biological
research. Freeman, New York.

Steegman.AT, Platner WS (1968) Experimental cold modification of cranio-facial morphology.
American Journal of Physical Anthropology, 28, 17-30.

Taylor AB (2002) Masticatory form and function in the African apes. American Journal of
Physical Anthropology, 117, 133-156.

Ulgen M, Baran S, Kaya H, Karadede I. 1997. The influence of the masticatory hypofunction on
the craniofacial growth and development in rats. American Journal of Orthodontics and
Dentofacial Orthopedics, 111, 189-198.

Varrela J (2006) Masticatory Function and Malocclusion: A Clinical Perspective. Seminars in
Orthodontics, 12, 102-109.

Watt DG, Williams CHM (1951) The effects of the physical consistency of food on the growth
and development of the mandible and the maxilla of the rat. American Journal of
Orthodontics and Dentofacial Orthopedics, 37, 895-928.

West-Eberhard MJ (2003) Developmental plasticity and evolution. Oxford University Press,
Oxford; New York.

Wolff JO, Dueser RD, Berry KS (1985) Food-habits of sympatric Peromyscus leucopus and
Peromyscus maniculatus. Journal of Mammalogy, 66, 795-798.

Young RL, Badyaev AV (2010) Developmental Plasticity Links Local Adaptation and
Evolutionary Diversification in Foraging Morphology. Journal of Experimental Zoology
Part B-Molecular and Developmental Evolution, 314B, 434-444.
138

Zelditch ML, Swiderski DL, Sheets HD, Fink WL (2004) Geometric morphometrics for
biologists: a primer. Elsevier Academic Press, San Diego, California; Amsterdam.

Zelditch ML, Wood AR, Bonett RM, Swiderski DL (2008) Modularity of the rodent mandible:
Integrating bones, muscles, and teeth. Evolution & Development, 10, 756-768.

Zelditch ML, Wood AR, Swiderski DL (2009) Building Developmental Integration into
Functional Systems: Function-Induced Integration of Mandibular Shape. Evolutionary
Biology, 36, 71-87.

139

CHAPTER IV

MORPHOMETRIC DIFFERENTIATION OF EXPANDING POPULATIONS OF PEROMYSCUS LEUCOPUS IN
THE NORTHERN GREAT LAKES: ENVIRONMENTAL AND GEOGRAPHICAL FACTORS

Abstract
Geographical range expansion associated with recent climate change is usually regarded as a
simple process of habitat tracking. Yet, the invasion of new environments is also expected to
produce phenotypic changes as populations adjust to the new local conditions. The identification
of the factors causing phenotypic variations is very challenging because of the complexity of
variables and interactions that mediate phenotypic change. Understanding the association
between morphological variation and environmental factors is crucial to improve our ability to
anticipate the evolutionary consequences of predicted climatic change scenarios. In this study,
patterns of genetic and morphometric variation were examined in light of geographical and
environmental variables, with the objective of identifying factors associated with both variation
in mandible shape among populations, and differentiation between old and new populations.
Distance-based methods are used to test for broad associations among morphological, genetic
and environmental variables, and then PLS regression was used to find actual directions of
variation responsible for these associations. The results of this study suggest that cold
temperatures and snow fall are the main factors underneath the differentiation between old and
new populations, while land cover contributes with variation that is not related to the expansion
process.

Introduction

140

It is often assumed that the main effect of current trends in climate change on the geographic
range of populations experiencing range expansion is to increase their available habitat, so that
organisms simply track the habitat change (Hughes 2000; Thomas et al. 2001; Walther et al.
2002; Salamin et al. 2010). However, this might not be possible for populations in temperate
zones, which are subject to seasonal changes and presumably adapted to multiple environmental
factors that are not necessarily directly affected by climate (e.g., photoperiod). The decoupling
between those environmental cues and the actual environmental conditions is, perhaps, one of
the most devastating effects of climate change, because it leads to the wrong interpretation of
seasonal cues, which may have detrimental effects on an organism’s fitness (Bradshaw &
Holzapfel 2006; Bronson 2009; Schlaepfer et al. 2002). Although climate is a key factor in
determining species geographical range (Pigot et al. 2010), changes in biotic interactions may
also play a crucial role in species range shifts (Van der Putten et al. 2010). In general, range
expansion associated with current climate change should not be viewed as a passive process of
tracking habitat, but as a complex process involving evolutionary mechanisms that would allow a
population to establish successfully in a new locality.
It would expected that the invasion of new environments is accompanied by changes in
phenotype, as populations readjust to new optima under the novel ecological circumstances,
either by local adaptation or plasticity (Chevin et al. 2010). Those phenotypic changes are
mediated by a number of factors, including genetic, developmental, and ecological, and the
complex interactions among them, as well as deterministic and stochastic evolutionary forces
(Monteiro et al. 2003; Poroshin et al. 2010; Straney & Patton 1980). In the case of small
mammals, these phenotypic changes in response to novel environmental conditions can take
place very rapidly (Goheen et al. 2003; Pergams & Ashley 1999; Pergams & Lacy 2008;

141

Pergams & Lawler 2009; Renaud & Auffray 2010), probably because of their short generation
length and lifespan.
Identifying the key factors underlying observed phenotypic variation is very challenging.
Most studies focus on the geographical variation in variables associated with climate (e.g.,
temperature, precipitation), topography (e.g., elevation), geographical distribution (e.g., latitude,
longitude), and/or vegetation type (e.g., forest vs. grassland), aiming to find the basis of the
morphological variation (Cardini et al. 2007; Monteiro et al. 2003; Poroshin et al. 2010; Rychlik
et al. 2006). Because phenotype is the final product of the interaction of evolutionary forces,
genes, and environment, each of those factors is just a piece of a very complex puzzle. A better
understanding of the association between morphological variation and ecogeographical factors
would be of value in climate change research, as it would improve our ability to anticipate some
of the evolutionary consequences of predicted climatic change scenarios.
The white-footed mouse, Peromyscus leucopus, has experienced a remarkable northward
range expansion in the Great Lakes region in recent decades. Expanding populations are not only
facing photoperiodic changes, but also changes in their environment, the most notable reflected
in the forest habitat, which transitions from deciduous in the south, to boreal in the north (Myers
et al. 2009). Morphometric analyses suggest that there are significant differences in mandible
shape between recently expanded populations (new populations) and long established ones (old
populations). These differences involve not only the magnitude of shape variation, but also its
dimensionality, indicating that rapid expansion has been accompanied by marked phenotypic
change (Chapter III). The pattern of shape variation is incongruent with the phylogeographic
pattern observed in these populations, suggesting a potential role of selection on the observed
differences in morphology. Results also indicate that a large proportion of the variation in

142

mandible shape is related to differences among populations. Regardless of whether these
phenotypic differences reflect plasticity or adaptation, or both, ecogeographical factors likely
underlie both the variation in mandible shape among populations, and differentiation between
old and new populations. In this study, I examine the association among morphometric, genetic,
and ecogeographical factors in expanding populations of Peromyscus leucopus and their putative
progenitors in the northern Great Lakes. My results clearly show that statements suggesting that
expanding populations are merely ―tracking‖ their habitat are not only oversimplifying a rather
complex process but also ignoring the fact that, in a strict sense, no two habitats are the same.

Materials and Methods
Samples
Twenty-four populations comprising 12 localities from Wisconsin (WI), 5 localities from
Michigan’s Lower Peninsula (LP), and 7 localities from Michigan’s Upper Peninsula (UP) were
sampled (Fig II.1, Chapter III). These populations were sorted by lineage, i.e., LP-EUP and
WUP-WI, according to results previously obtained from SAMOVA and haplotype genealogy
analyses (Chapter II). Within each lineage, populations were classified as ―old‖ or ―new‖
depending on whether a population has been long or recently established, respectively (Table
A.IV.1). Morphometric data were obtained from the mandibles of approximately 10 specimens
per locality, for a total of 227 specimens (Table A.III.1). Genetic data were obtained from 9-20
individuals per locality (Table II.1, Chapter II). DNA sequences and landmark configurations
were obtained as described elsewhere (Chapters II and III). Trapping sites within a 25 km radius
and located within similar environments were considered to represent a single population. When
more than one location satisfied these criteria, the geographical midpoint of the combined sites

143

was used as their common geographical position; this point was computed using the web-based
tool at http://www.geomidpoint.com. In the text, populations are also referred to as localities.
All animals trapped for this study were handled according to protocols approved by the
Michigan State University Institutional Animal Care and Use Committee, application # 08/04108-00, in compliance with the American Society of Mammalogists guidelines for the use of
wild mammals in research (Gannon et al. 2007).

Ecological Data
Land Cover/Use
Using geographical midpoints as reference locations, custom-made land cover/ use maps for
each locality were obtained from aerial photographs (Fig. IV.1). These maps were elaborated by
a GIS/Remote Sensing Analyst (from Michigan State University Remote Sensing & Geographic
Information Science, Research and Outreach Services) by overlaying a 454 m-radius
(approximately 0.65 km2) circle on aerial images from the National Agriculture Imagery
Program (NAIP). NAIP is administrated by the Farm Service Agency of the Department of
Agriculture and charged with acquiring aerial imagery during agricultural growing seasons in the
continental United States. The maps used for this study were based on photographs taken in
2007, and the radius of the focal circle was chosen to match the home range of the white-footed
mouse (Burt 1940).
Land cover/use polygons were identified and drawn on aerial photographs for each
locality, and then classified according to the Michigan Resource Information System (MIRIS)
2007 land cover/use classification system, which is based on the classification system by
Anderson et al. (1976) for use with remote sensing data. In this way, the relative area covered by

144

each land cover category found in each locality was obtained. This system has four levels of
classification depending on the resolution (i.e., altitude and scale) of the data: Level I contains
satellite image data (i.e., LANDSAT) and offers the most comprehensive coverage; Level II
consists of high-altitude data (12,400 m) with a scale below 1:80,000; Level III provides
medium-altitude data (3,100 – 12,400 m), spanning a scale between 1:20,000 and 1:80,000; and
Level IV contains low-altitude data (i.e., aerial photographs taken below 3,100 m) with a scale
above 1:20,000, resulting in highly detailed information (Anderson et al. 1976). For this study,
the land cover categories (herein also referred to as environments) Urban, Agricultural, Grass
and Shrub, Water, and Barren were classified up to Level II, and the categories Forest and
Wetland were classified up to Level III.

Climate Data
The closest climate station to each of the sample localities was located and historical climate data
for the period of 1971 to 2000 were obtained from the Midwest Regional Climate Center
(http://mrcc.isws.illinois.edu/index.jsp). The following climate variables were analyzed in this
study: mean annual, winter, spring, summer and fall temperatures; annual snow precipitation,
and annual rain precipitation (Table A.IV.1). Seasonal temperatures were calculated as the
average of the mean monthly temperatures during the months of December, January, and
February for Winter; March, April, and May for Spring; June, July, and August for Summer; and
September, October, and November for Fall.

Index of Environmental Similarity
To compare localities in terms of their environmental compositions, Morisita’s Index was

145

computed for each pair of localities (Brower et al. 1990). This index measures similarity among
communities as the probability that two distinct environments sampled in two localities will be
the same, equaling zero when two communities are totally dissimilar and one when they are
completely similar. It is computed using

⌈

∑

⌉

where xi is the area of environment i in locality 1; yi is the area of environment i in locality 2; N
is the total area of each locality; and l1 and l2 are Simpson’s indices of dominance for localities 1
and 2, respectively. Simpson’s index incorporates information about the relative abundance of
environment types, and is computed using the following equation:

∑

where xi is the area of environment i and N is the total area of each locality (Table A.IV.1). This
index can be interpreted as the probability that two randomly selected environments at one
locality will be the same. It varies from zero, indicating maximum diversity, to one, indicating
absence of diversity, i.e., only one type of environment at that locality (Brower et al. 1990).
To convert these measurements into ecological distances, a matrix of pairwise distances
was constructed with elements equal to (1 - IM), with values ranging from zero, indicating that
the two localities possess identical environments, to 1, indicating that the two localities have no
environments in common (Table A.IV.2).

Statistical Methods
The following methods were used to assess statistical associations between geographical and

146

environmental variables, and observed patterns of genetic and morphometric variation.

Association between shape, geographical, ecological, and genetic distances
Partial Mantel tests (Manly 1986) were used to assess associations among shape, geographical,
ecological, and genetic distances among sampled populations. This test extends the basic Mantel
test for correlation between two matrices by using the rationale of multiple regression (Thorpe et
al. 1995). In a partial Mantel test, an independent, response variable matrix is simultaneously
compared to multiple dependent predictor matrices, after each has been controlled statistically
for the effects of the others. This approach was used to evaluate (1) the regression of
morphometric distances (both mandible shape and size) on geographical, ecological, and genetic
distances, and (2) the regression of genetic distances on geographical and ecological distances.
For these tests, shape distances were measured as pairwise Euclidean distances between
population means computed from average Procrustes residuals projected onto a space of
adequate dimensionality (see Chapter II for details), and size distances were measured as
pairwise Euclidean distances between population mean centroid sizes. Geographical distances
between the geographical midpoints of populations were calculated as geodesics, i.e., shortest
distances over the Earth’s surface (Table A.IV.3). Ecological distances were computed as 1- IM,
where IM is the Morisita’s index of community similarity described above (Table A.IV.2). Two
different genetic distance metrics were considered in these analyses, pairwise Nei’s distances D
(Nei 1987) and Fst estimates (Reynolds et al. 1983), which were computed using Arlequin 3.11
(Excoffier et al. 2005). Significance levels were adjusted by Bonferroni correction for multiple
comparisons. Additionally, dendrograms were constructed using an average distance-based
linkage function (UPGMA method) to graphically visualize distance patterns in these matrices.

147

Effect of climate, geography, and land cover on mandible shape
The pattern of covariance between mandible shape and ecogeographical variables (i.e., climate,
geography, and land cover) was investigated using two-block Partial Least Squares (2B-PLS).
Two-block PLS is a regression technique in which two multivariate sets of variables are
decomposed into pairs of linear combinations that maximally account for the covariance between
the original variables (Rohlf & Corti 2000). PLS resembles ordinary regression in that both
methods produce coefficients that can be used to describe associations between two sets of
variables. However, unlike regression, PLS does not specify a predictor-response structure
between the sets of variables, but instead treats them symmetrically (Zelditch et al. 2004).
Computationally, PLS is also similar to Principal Component Analysis (PCA) in that both
methods perform an eigendecomposition of the variance-covariance matrix, and thus produce a
series of orthonormal axes that each account for a progressively smaller proportion of the
original information. The two methods differ, however, in that in PCA, the variance-covariance
matrix among all variables is decomposed, thus producing a single set of eigenvectors, whereas
PLS produces more than one set of eigenvectors, as it is based on the matrix of cross-covariances
between variables in multiple data blocks (Zelditch et al. 2004). In two-block PLS, pairs of axes
are produced, each with an associated singular value (λi) that measures the proportion of the
squared covariance between the two data blocks that is captured by those axes. To assess
consistency between blocks, correlations between corresponding PLS axes can be measured
using a standard correlation coefficient. In each analysis, interpretation was restricted to the set
of PLS axes capturing a cumulative total of 90% of the square covariance between the blocks.
To investigate the covariation between mandible shape and climate, two PLS analyses

148

were conducted: first, shape variables were regressed on a composite matrix comprising annual,
winter, spring, summer, and fall average temperatures; and second, shape variables were
regressed on a composite precipitation matrix containing total annual snow- and rainfall records.
These two sets of climatic variables were analyzed separately because temperature and
precipitation are expressed in different units, and standardization would not preserve their
variation structure.
To investigate the pattern of covariance between mandible shape and geography and
environment, two additional PLS analyses were conducted. Here shape variables were regressed
on the geographical coordinates (i.e., latitude and longitude) of the midpoint of each population,
and on the population environmental variable (defined by the relative abundance of the 20 land
cover/use types identified in Fig. IV.1).
Analyses were performed in Matlab R2006a (The MathWorks, Inc.), using a combination
of tools for multivariate and graphical analysis included in Matlab’s Statistics Toolbox.

Results
Association between morphometric, geographical, ecological, and genetic distances
Partial Mantel tests reveal different patterns of correlations between mandible shape and
geographical, ecological, and genetic factors for the two lineages analyzed in this study (Table
IV.1). Because of a very high positive correlation between the two metrics of genetic distance
computed herein, i.e., Nei’s D and Fst (r = 0.98, P < 0.001), only the latter is reported. In the LPEUP lineage, shape distances between populations are not significantly correlated with
geographical, ecological, or genetic distances when the others factors are held constant (P >
0.016; Table IV.1). In contrast, shape distances between populations in the WUP-WI lineage are

149

significantly correlated with all three of these factors when the others are held constant (P <
0.003, Table IV.1). Pooling the two lineages and repeating these tests results in high positive
correlations between shape distances and both geographical and ecological distances (P < 0.001),
and a non-significant positive correlation between shape and genetic distances (P = 0.03). The
differences between lineages found here might be explained in part by the much larger WUP-WI
sample (18 populations versus 6 for LP-EUP), but might also reflect historical and ecological
differences between lineages.
The relationships between mandible size distances and geographical, ecological, and
genetic distances are broadly consistent across lineages (Table IV.2). Size distances are
significantly correlated with ecological distances (P < 0.001), but not with geological or genetic
distances. This result suggests that mandible size, and possibly overall body size, is responsive to
ecological differences, after geographical and genetic differences have been taken into account.
Despite significant positive correlations, there are no obvious patterns of correspondence
between dendograms representing the three sets of distances for either lineage (Fig. IV.2), i.e.,
dendograms based on morphometric (shape), genetic, and ecological distances do not resemble
one another.

Association between genetic and geographical and ecological distances
Geographical distances are positively and significantly correlated with genetic distances (based
on pairwise Fst). This result is concordant with the pattern of isolation by distance suggested by
genetic analyses (Chapter II). In contrast, pairwise genetic distances are not significantly
correlated with ecological distances, after controlling for geographical distances, suggesting that
environmental similarity does not affect gene flow in a predictable manner (Table IV.3).

150

Effect of climate on shape
Table IV.4 shows results from PLS analyses of multivariate covariation between shape and
temperature variables within each lineage. For the LP-EUP lineage, PLS produces an essentially
one-dimensional solution, where a single axis (labeled PLS1) explains 96% of the squared
covariance between the two sets of variables, although only a small proportion of the square
covariation is captured by those axes (λ1 = 0.0320). This dimension receives similar
contributions from most temperature variables, and therefore can be broadly interpreted as an
axis of variation in overall temperature. The largest contribution is from winter temperature,
which shows a slightly higher coefficient along this vector than the other temperature variables.
A plot of the temperature dimension (y-axis) against the shape dimension (x-axis, Fig. IV.3a),
contrasts mandible shape in localities with high temperature values (e.g., Charlevoix, CX) to that
in localities with relatively low temperature values (e.g., Chippewa, CW). Increasing
temperatures are associated with expansion of the angular and condyloid processes, contraction
of the caudal end of the masseteric ridge, and slight contraction of the incisor alveolus.
For the WUP-WI data, PLS supports a two-dimensional solution representing covariation
between shape and temperature variables, with both axes being required to account for 90% of
the squared covariance (Table IV.4). The square covariation captured by those axes is also small
(λ1, LP-EUP = 0.0105 and λ1, WUP-WI = 0.007). The first of these axes is very similar to PLS1 for
the LP-EUP lineage, whereas the second contrasts temperature in winter to that in the warmer
seasons. A plot of the temperature dimension against the shape dimension for PLS1 (Fig. IV.3b)
contrasts mandible shape in localities with high temperature values (e.g., Buffalo, BF) to that in
localities with relatively low temperature values (e.g., Gogebic, GO). Increasing temperatures are
151

associated with a slight expansion of the condyloid; and contraction of the coronoid, the caudal
margin connecting condyloid and angular processes, the masseteric ridge, and the molar
alveolus. The temperature axis for PLS2 (Fig. IV.3c) contrasts localities where winters are
milder and springs and summers are generally cooler (e.g., Marquette, MA), with localities
where winters are severe and spring and summer temperatures are warmer (e.g., Buffalo, BF).
Shape changes in transitioning from the former to the latter include the molar and incisor alveoli,
and the caudal margin connecting condyloid and angular processes, the latter resembling
variation in the coronoid process captured by PLS1 (Fig. IV.3c).
Analyses of covariation between shape and precipitation variables suggest that snowfall
has consistently had a much stronger effect on mandible shape than rainfall in both the LP-EUP
and WUP-WI lineages (Table IV.5). Increasing snowfall in LP-EUP populations is associated
with contraction of the anterior angular and condyloid processes, expansion of the coronoid
process, and expansion of the caudal end of masseteric ridge relative to the base of the mandible
(Fig. IV.4). The opposite pattern is found in the WUP-WI, where increasing snowfall is
associated with expansion of the condyloid process and molar alveolus, and contraction of the
caudal end of masseteric ridge. A larger proportion of the square covariation between shape and
precipitation is captured by these axes, in comparison to the covariation between shape and
temperature captured by PLS (λ1, LP-EUP = 0.2413 and λ2 = 0.1247).

Effect of geographical location on shape
PLS analyses of covariation between shape and geographical variables produced contrasting
results for the two lineages; in the LP-EUP, shape covaries almost exclusively with latitude,
whereas in the WUP-WI, shape covaries strongly with longitude, and only weakly with latitude

152

(Table IV.6). The square covariation captured by those axes is also small (λ1, LP-EUP = 0.0029;

λ1, WUP-WI = 0.0048 and λ2, WUP-WI = 0.0030). As expected, the effect of latitude on mandible
variation in LP-EUP populations is very similar to the effect of temperature (Figs. IV.3a and
IV.5a). In WUP-WI, movement from east to west (i.e., PLS1, Fig. IV.5b) is associated with
narrowing and straightening of the coronoid process and sharpening of the curvature of the
caudal margin connecting the angular and condyloid processes. Unlike what is observed in LPEUP populations, the effect of latitude in the WUP-WI lineage (i.e., PLS2, Fig. IV.5c) does not
resemble the effect of temperature, but is instead characterized by a reduction in curvature of the
caudal margin connecting the angular and condyloid processes, dorsal expansion of the
condyloid process, anterior displacement of the caudal masseteric ridge, and expansion of the
molar alveolus in transitioning from southern (e.g., Adams, AD) to northern (e.g., Marquette,
MA) localities (Fig. IV.5c).

Effect of land cover on shape
Results from a PLS analysis of covariation between shape and land cover variables are shown in
Table IV.7. The 90% threshold of the squared covariance between the two factors for the LPEUP lineage was obtained in two dimensions. The environments mostly responsible for the
association between shape and land cover on PLS1 are northern hardwoods and pine, which
appear to have opposite effects on mandible shape. On PLS2, northern hardwoods and pine both
have high positive loadings that contrast with high negative loadings for spruce-fir and lakes.
The 90% threshold of the squared covariance between the two factors for the WUP-WI
lineage was obtained in four dimensions. The dominant types in the first two dimensions (i.e.,
PLS1 and 2) are northern hardwoods and conifers, which show a similar pattern of covariation to
153

that found in the LP-EUP data. The third and fourth dimensions contrast central hardwoods to all
other environment types (PLS3), and pine to lowland conifers and scrubs/ shrubs/wetlands
(PLS4).
Plots of the land use dimension (y-axis) against the shape dimension (x-axis) suggest that
vegetation type has a substantial effect on mandible shape (Fig. IV.6). In the LP-EUP lineage,
there is a nearly linear gradient of environment-shape covariation for PLS1 (Fig. IV.6a), from
localities characterized primarily by northern hardwoods (e.g., Charlevoix, CX) to those where
pine is prevalent (e.g., Crawford, CR). The former habitats are associated with mandibles that are
more robust, with more prominent processes, and a proximal contraction of the mandibular
alveoli. PLS2 (Fig. IV.6b) contrast regions where hardwoods and pine are abundant (e.g.
Chippewa, CW) to regions characterized by lakes and spruce-fir forests (e.g., Cheboygan, CH).
Shape changes in transitioning from the former to the latter include expansion of the condyloid
and angular processes and contraction of the caudal end of the masseteric ridge and coronoid
process.
In the WUP-WI lineage, PLS1 (Fig. IV.6c) contrasts localities in which conifers are the
prevalent vegetation and hardwoods are scarce (e.g., Marquette, MA) to those with the opposite
configuration (e.g. Ontonagon, ON). Shape change along this axis is characterized by contraction
of the angular and condyloid processes, posterior displacement of the caudal end of the
masseteric ridge, contraction of the molar and ventral displacement of the frontal end of the
masseteric ridge. PLS 2 (Fig. IV.6d) contrasts environments where pine is abundant (e.g.,
Ontonagon, ON, and Marquette, MA) to localities where northern hardwoods and spruce-fir
prevail (e.g., Schoolcraft, SC). The change in mandible morphology along this axis includes
expansion of the dorsal margin connecting the coronoid and the condyloid processes and the

154

caudal end of the masseteric ridge, and contraction of the condyloid, and molar and incisor
alveoli.
The next dimension (PLS3, Fig. IV.6e) separates Marinette (MR) from all other
populations, as the only one where central hardwoods are abundant. The major effects on
mandible shape associated to this type of forest are increased curvature of the caudal margin
connecting condyloid and angular processes, and contraction of the caudal end of the masseteric
ridge, and of the molar and incisor alveoli, yielding a relatively shorter mandible. The last
dimension (PLS4, Fig. IV.6f) separates populations where pine forests are more abundant (e.g.,
Schoolcraft, SC) from those where lowland conifers, scrubs/shrubs/wetland predominate (e.g.,
Outagamie, OU, and Buffalo, BF). Morphological changes along this axis include expansion of
the coronoid process and of the curvature of the caudal margin connecting the angular and
condyloid process, and contraction in the angular process.
In all cases, the square covariation captured by those axes is also small (λ1, LP-EUP =
0.0025 and (λ2, LP-EUP = 0.0011; λ1, WUP-WI = 0.0012 and λ2, WUP-WI = 0.0010, λ3, WUP-WI =
0.0006, λ2, WUP-WI = 0.0005).

Discussion
As shown in previous analyses, as much as 85% of the variation observed in the shape of
mandible within each lineage was explained by differences among populations, including
genealogy, expansion history (i.e., old established vs. newly founded), and characteristics of the
populations, such as local ecology (Chapter III). In this Chapter, some aspects of the local
ecology were analyzed to assess their role in the variation in shape observed across populations.
As expected, results from this study show a consistent association between shape differentiation
155

and latitudinal distribution and temperature and precipitation variables, whereas land cover is
unrelated to variation accumulated during the expansion process. The following discussion
focuses on the relationship between those geographical and environmental variables and
observed changes in the shape of the mandible, with particular emphasis on the potential role of
these factors in explaining the differentiation between old and new populations.

Distance-based comparisons between shape and geographical, ecological, and genetic variables
Partial correlations between morphometric distances and geographical, ecological, and genetic
distances among populations in the LP-EUP lineage were not significant in any case, which
suggests that similarity in mandible shape among populations is not strongly influenced by
geographical proximity, environmental similarity, or gene flow. Alternatively, lack of statistical
power due to small sample sizes could hinder the discovery of patterns of association among the
examined variables, which remains plausible given the statistical significance of the same
associations in the better-sampled WUP-WI lineage.
Results based on WUP-WI data suggest greater similarity in mandible shape among
specimens from localities separated by shorter distances, after controlling for the effects of
genetic and ecological distances. This pattern of positive correlation between morphometric and
geographical distances suggests that gene flow has played a role in the process of phenotypic
divergence, which is consistent with expectations of a model of isolation by distance, and further
supported by observed positive correlations between morphometric and genetic distances (see
below). Mice from ecologically similar localities also have mandibles that are more similar than
mice from ecologically dissimilar localities, both in shape and size, after controlling for the
effects of genetic and geographical distances. A number of factors might explain this association.

156

For example, mice inhabiting similar habitats are likely to exploit similar resources, which could
lead to morphological similarity due to either shared epigenetic input and interactions during the
development of the mandible (Atchley 1993; Atchley & Hall 1991). Similarity could be due to
developmental plasticity (Schlichting and Pigliucci 1998), shared selective pressures (Rychlik et
al. 2006), or both.
Morphometric and genetic distances in the WUP-WI lineage appear to be significantly
correlated after partialing out the effects of ecological and geographical distances. This contrasts
with previous results from pairwise correlation analyses, which showed these distances to be
uncorrelated (Fig. IV.2; see also Chapter III). Similarly, pairwise correlations between
morphometric and ecological distances are non-significant before partialing out the effects of
genetic and geographic variation (r = 0.27, P = 0.33; data in Appendices IV and VI). In both
cases, these results suggest a complex pattern of association between causal factors and
phenotypic variation, which is likely a consequence of the high dimensionality of the latter. As
demonstrated by PLS analyses (see below), different dimensions of mandible shape variation
may be associated with disparate factors, and patterns of distances among populations may
reflect the added effect of multiple causal agents. Thus, whereas shape differences over all
possible dimensions may be uncorrelated with genetic or ecological factors, a multiple-predictor
analysis such as the partial Mantel tests carried out herein suggests that there are dimensions of
shape variation that, in fact, can be explained by either of these factors. Therefore, these results
complement the findings reported in Chapter III by suggesting that at least some aspects of
variation in mandible shape are consistent with expectations from phenotypic drift, whereas
others are consistent with ecological differentiation.

157

Genetic and geographical distances show a significant positive partial correlation, after
controlling for the effect of ecological distance, in agreement with previous results (Chapter II).
Because genetic distances are based on mtDNA polymorphism, which is maternally inherited,
this pattern of isolation by distance might be the consequence of the limited dispersal capabilities
and philopatry of females in Peromyscus leucopus, which supports results from field studies that
have found that dispersal in this species is male biased (Jacquot & Vessey 1995), as it is usually
seen in mammals (Johnson 1986). The pattern of isolation by distance associated with
morphometric distances is caused by a different process. For all morphometric analyses, males
and females were pooled together, since there were no differences between sexes (Chapter III).
Thus, that pattern is not caused by female philopatry, but might be explained by the fact that, in
general, Pln does not typically disperse over great distances (Krohne et al. 1984; but see Maier
2002 for a report of an exceptional dispersion distance in this mouse).
The absence of association between genetic and ecological distances, after controlling for
the effect of geographical distance, suggests that gene flow among these populations might not
be limited by the characteristics of environments of this regions. This result is supported by
evidence that Peromyscus leucopus is able to disperse over great distances (Maier 2002), and
presumably exchange genes, regardless of whether the habitat is continuous or fragmented
(Mossman& Waser 2001). However, this is not true for urbanized areas, like New York City,
where populations have very limited dispersal ability and are genetically differentiated (MunshiSouth & Kharchenko 2010).

Effect of climate, geography, and land cover on shape
In this study, Partial Least Squares regression was used to find the dimensions of shape variation
158

within each of the sampled lineages that most strongly covary with a number of geographical and
environmental variables. These dimensions can be interpreted as components of shape variation
that have been chosen to maximize the explanatory power of specific extrinsic predictors, in
contrast to the distance matrix approaches discussed above in which associations are studied
without prior decomposition.
Thus, even though comparisons based on distance matrices in the LP-EUP lineage
returned no significant associations, PLS decomposition found directions in shape space that are
highly correlated with climatic and environmental factors. Three of the considered predictors,
temperature, precipitation, and geographical location, capture similar signals, with CW
population bearing the most extreme scores for snowfall, average annual temperature, and
latitude. Along these dimensions, the temperature and snowfall axes describe similar patterns of
contraction of the angular and condyloid processes, accompanied by an anterior expansion of the
alveoli and a dorsal expansion of the masseteric ridge at colder and high-precipitation localities.
The latitude axis reflects some of the same patterns, i.e., contraction of the angular process
coupled with dorsal displacement of the masseteric ridge at northernmost localities, whereas
others, e.g., contraction of the condyloid process, were not seen.
The contraction/expansion of one or more mandibular processes is a recurrent pattern
captured by these predictors. The mandibular processes were the regions of the mandible that
captured most of the differences between old (LP) and new (CW) populations (Chapter III).
Together, these results suggest that climate rather than vegetation may be the most important
driver behind such divergence. In contrast, land cover (Table 7) captures the contrast between a
few major types of environment and is not strongly associated with geographical location. The
extreme value obtained for the CW population along the second axis, however, deserves further

159

attention upon acquisition of additional samples from the UP region, which might reveal a more
discrete separation between old versus new populations than is presently possible.
As previously discussed, distance analyses based on WUP-WI mandible data suggest that
only some dimensions of shape variation are correlated with geographical distribution (i.e.,
geographical distances), whereas other dimensions are correlated with environmental features
(i.e., ecological distances). PLS complements these results by providing estimates for some of
these dimensions. The novel application of fine-grained land cover/use maps to explore the
association between patterns of morphometric and habitat similarity has proven a valuable tool
for gaining insight regarding details of habitat structure that could be responsible for gross
patterns of similarity suggested by correlation analyses in this study.
Unlike PLS results for the LP-EUP lineage, PLS axes based on the WUP-WI shape data
do not show a predominant contrast among populations or a recurrent motif of shape variation.
Instead, predictors are associated with different, albeit partly overlapping, aspect of shape
variation. For instance, annual snowfall is associated with variation in the posterior masseteric
ridge, where average temperature has little effect, instead showing a marked influence on the
coronoid and angular processes; on the other hand, both of these factors exert an effect on the
molar alveolus.
The comparison between the mean shapes of old and new populations in the WUP-WI
lineage (Chapter III, Fig. III.6) indicates that the most prominent changes in the new populations
correspond to an expansion of the condyloid process, and a contraction of both the angular
process and the molar alveolus. These changes closely correspond to dimensions of variation
captured by oscillations in annual temperatures (Fig. IV.3a) and snow precipitation (Fig. IV.4a),
and are only partly related to patterns of shape variation associated with latitude (Fig. IV.5c). A

160

simple ANOVA testing the similarity between old and new populations for the projections of
morphometric data on the PLS axis of temperature, precipitation, and geographic predictors
reveals that PLS1 significantly accounts for differences between these groups (F = 6.93, P =
0.018), whereas PLS2 does not (F = 1.88, P = 0.190). Snowfall (captured by axis 2 of
precipitation-shape PLS) fails to reject the null hypothesis of equality of means (F = 4.02, P =
0.062), possibly due to the presence of localities with both extremely high and extremely low
snowfall rates among the new populations, whereas most of the populations that define this axis
(e.g., MA, ON, GO, AS, IR), are in fact new. Finally, latitude (F = 10.78, P = 0.005), but not
longitude (F = 0.05, P = 0.827), shows an association with the expansion process. Shape
variation associated with land cover/use, although displaying a contrast among environment
types (i.e., hardwood vs. conifer forests), is largely unrelated to the expansion process of the LPEUP or WUP-WI lineages, a result confirmed by ANOVAs on PLS scores (P ≥ 0.05 in all axes).
Cold temperature and snow fall have also been causally linked to changes in the shape of
the skull and mandible in the common shrew Sorex araneus over a period of only 27 years,
presumably these reflect changes in foraging strategies (Poroshin et al. 2010). The effect of cold
temperature on morphology has been investigated in several previous studies, especially the
association between temperature gradients (or latitude) and postcranial morphology (Ashton et
al. 2000; Heath 1984; Rae et al. 2006; Steegman & Platner 1968; Weaver & Ingram 1969).
Craniofacial features change in response to cold temperature (Nicholson & Harvati 2006; Perez
& Monteiro 2009; Poroshin et al. 2010; Rae et al. 2006; Steegman & Platner 1968). Those
morphological changes can be induced very rapidly, as has been shown experimentally by
Steegmann and Platner (1968). In their study, animals reared under cold stress (5º C) showed a
reduction in the surface of insertion of masticatory muscle in comparison to animals grown under

161

ambient temperature (22º C). Although those results might have limited explanatory power due
to the simplicity of the analysis, they suggest possible changes in the mandible associated with
colder temperatures. The results obtained in this study also show a reduction in mandibular
processes, which are the places of insertion of masticatory muscles, in the new populations of
Pln (from both lineages) that are found in colder localities.
Peromyscus leucopus is an omnivore that feeds mostly on seeds and arthropods (Wolff et
al. 1985), and its dietary breadth is determined not only by the abundance (Ivan & Swihart
2000), but also the nutritional value of seeds (Lewis et al. 2001). Animals in colder environments
have higher metabolic requirements that demand eating (and chewing) more often, and because
of the scarcity of resources during the colder seasons, these animals might increase their diet
breadth and/or feed on less-preferred resources during periods when their preferred foods are
unavailable (i.e., fallback foods ; Lambert 2009; Marshall & Wrangham 2007). Hence, a possible
explanation for the reduction in mandibular processes observed in specimens from new
populations in colder climates could be that these animals are feeding on softer items (e.g.,
insects, softer-shell seeds) that might offer a higher net nutritious value. There is experimental
evidence showing that animals fed on a soft diet show a significant reduction in mandibular
processes reflecting the lower masticatory strength required to chew softer foods (Corruccini &
Beecher 1982; Maki et al. 2002; Renaud et al. 2010), which is consistent with patterns observed
in this study.
Even though land cover is a poor predictor of differences between old and new
populations, more detailed community characterizations could help to expose causal local factors
behind observed morphological divergence. An obvious candidate is diet, which is considered a
determinant factor in the diversification of primates (Anapol and Lee 1994), bats (Nogueira et al.

162

2009; Monteiro and Nogueira 2010), and rodents ((Michaux et al. 2007; Samuels 2009).
Incorporating detailed dietary information with geographical expansion and seasonal and nonseasonal climatic fluctuations, might significantly improve our understanding of the causal
factors behind phenotypic divergence via local adaptation/acclimation to novel environmental
regimes.

Conclusions
Analyses carried out in this study attempted to dissect patterns of morphometric and genetic
variation in the light of geographical and environmental variables. Distance-based methods were
used to test for broad associations among these variables, and then PLS regression was used to
find actual directions of variation responsible for these associations. Together, these analyses
suggest that among the considered factors, climatic variables, specifically temperature and snow
precipitation, are the best predictors for the divergence between old and new populations.
However, it remains unclear whether these factors have a direct effect on morphology, or if
alternatively, their association with shape reflects the latitudinal component of the geographic
expansion. Shape deformations suggest that the latter could be partly the case, but the majority of
mandible changes associated with temperature and precipitation are unrelated to latitude. Besides
these climatic factors, vegetation and land use variables failed to show an association with the
expansion, suggesting that they instead contribute to non-geographical variation. Previous results
showed that differences among populations explain over 80% of the overall variation in shape;
the results reported here indicate that this component of variation might be partly accounted for
by the differences in land features mapped in this study. These effects, however, account for
disparate directions of variation. Given that the expansion process is a relatively recent

163

phenomenon, it is plausible that variation in shape is in fact highly noisy and difficult to predict
based on regional-scale factors. Careful modeling and increasing the level of detail, for instance
by adding species-level community information, could be a promising approach to explain
morphological differences among populations.
I found substantial differences between the two lineages in their responses to
geographical and environmental factors, further supporting previous findings that the expansion
process has differed between them (Chapters 1 and 2). In the LP-EUP lineage, the observed
association with climate variables is largely due to the extreme differences observed between the
new (CW) and old populations (LP), not only in terms of genetic and morphological
differentiation (see Chapters II and III), but also in relation to geographical and climatic factors.
In contrast, in the WUP-WI lineage, old and new populations are not sharply differentiated,
perhaps because they are more continuously distributed that in the LP-EUP lineage. In this case,
the association between morphological and genetic variation and geographical and climatic
factors is more gradual.
Have populations of Pln in the northern Great Lakes been tracking their habitat during
their recent geographical expansion? These results suggest that morphological variation in the
mandible has responded to climatic shifts encountered in the newly occupied localities, which
counters the notion of habitat tracking. Similarly, other environmental factors, such as land
cover/use, contribute to population differences, although such differences are unrelated to
geographical expansion. Combined, these findings are consistent with the fact that both intrinsic
(phenotypic, genetic) and extrinsic (geographical, environmental) factors comprise multifactorial
phenomena, whereas the notion of habitat tracking is better suited for one-dimensional responses
to one-dimensional environmental shifts. For instance, whereas climate encompasses such

164

meteorological elements as temperature, precipitation, atmospheric pressure, wind, and humidity,
phenotypes are also multidimensional and therefore can simultaneously respond to disparate
selective pressures along multiple directions of variation, such as those found by PLS analyses.
Therefore, the suggestion of some authors that range expansion in response to current climate
change has been possible due to populations tracking the expansion of their habitat (Hughes
2000; Thomas et al. 2001; Walther et al. 2002; Salamin et al. 2010) is an oversimplification of a
complex and multidimensional process for which passive tracking may be a reasonable
assumption only for a limited set of traits.

165

Tables

Table IV.1. Relationship between mandible shape distances and geographical (Geo), ecological
(Eco), and genetic (Fst) distances based on partial Mantel tests. See text for details about how
distances were calculated. The correlation (r) is between shape distance and the corresponding
distance variable after controlling for the other two factors. Tests of significance were based on
1,000 permutations. The Bonferroni-adjusted alpha level is 0.006. Significant correlations are
indicated in boldface.

Factor
Geo
r

Eco
P

r

Fst
P

r

P

LP-EUP

0.128

0.319

0.615

0.016

0.324

0.120

WUP-WI

0.506

< 0.001

0.483

< 0.001

0.317

0.003

Pooled

0.349

< 0.001

0.582

< 0.001

0.185

0.030

166

Table IV.2. Relationship between mandible size distances and geographical (Geo), ecological
(Eco), and genetic (Fst) distances base on partial Mantel tests. See text for details about how
distances were calculated. The correlation (r) is between size distance and the corresponding
distance variable after controlling for the other two factors. Significance tests were based on
1,000 permutations. The Bonferroni-adjusted alpha level is 0.006. Significant correlations are
indicated in boldface.

Factor
Geo

Eco

r

P

r

LP-EUP

0.154

0.284

0.870

WUP-WI

0.140

0.075

Pooled

0.128

0.097

Fst
r

P

0.002

0.015

0.450

0.386

< 0.001

0.178

0.073

0.409

< 0.001

0.092

0.191

167

P

Table IV.3. Relationship between genetic distances based on pairwise Fst and geographical
(Geo) and ecological (Eco) distances base on partial Mantel tests. . See text for details about how
distances were calculated. The correlation (r) is between genetic distance and the corresponding
distance variable after controlling for the other factor. Significance tests were based on 1,000
permutations. The Bonferroni-adjusted alpha level is 0.0083. Significant correlations are
indicated in boldface.

Geo
r

Eco
P

r

P

LP-EUP

0.912

0.004

-0.559

0.029

WUP-WI

0.524

0.001

0.104

0.231

Pooled

0.767

0.001

-0.121

0.151

168

Table IV.4. Results of a partial least-square analysis of the covariation between mean mandible
shape and temperature for populations of Peromyscus leucopus in the LP-EUP and WUP-WI.

LP-EUP

WUP- WI

Temperature variables
PLS1

PLS1

PLS2

Annual Temp

0.442

0.442

0.063

Winter Temp

0.573

0.563

-0.676

Spring Temp

0.398

0.378

0.556

Summer Temp

0.350

0.392

0.470

Fall Temp

0.443

0.437

-0.096

0.0320

0.0104

0.007

95.95

68.90

29.18

0.94

0.66

0.75

Singular value
%Explained covariance
Correlation

169

Table IV.5. Results of a partial least-square analysis of the covariation between mean mandible
shape and precipitation for populations of Peromyscus leucopus in the LP-EUP and WUP-WI.

LP- EUP

WUP-WI

PLS1

PLS1

Precipitation variables
Snowfall

0.998

1.000

Rainfall

0.069

-0.002

0.2413

0.1247

99.91

99.87

0.89

0.64

Singular value
%Explained covariance
Correlation

170

Table IV.6. Results of a partial least-square analysis of the covariation between mean mandible
shape and geographical coordinates for populations of Peromyscus leucopus in the LP-EUP and
WUP-WI.

LP-EUP

WUP- WI

Geographical variables
PLS1
Latitude

PLS1

PLS2

1.000

0.081

-0.997

Longitude

-0.023

0.997

0.081

Singular value

0.0029

0.0048

0.0030

90.75

72.05

27.95

0.82

0.80

0.79

%Explained covariance
Correlation

171

Table IV.7. Results of a partial least-square analysis of the covariation between mean mandible shape and land cover/use for
populations of Peromyscus leucopus in the LP-EUP and WUP-WI.

Land use variables

LP- EUP
PLS1

WUP-WI

PLS2

PLS1

PLS2

PLS3

PLS4

Residential

0.499

3.730

0.675

-1.984

-14.543

-1.656

Commercial, Services and Institutional

0.000

0.000

-3.303

-1.682

0.163

2.032

-1.074

-0.323

0.000

0.000

0.000

0.000

Extractive

0.000

0.000

0.657

0.372

-1.174

0.429

Open Land and Other

0.000

0.000

-8.492

-4.922

-0.322

4.701

Cropland

0.000

0.000

-1.694

-10.246

-0.563

-0.429

Permanent Pasture

1.651

3.178

-1.266

-7.513

1.635

1.900

Other Agricultural Land

0.000

0.000

-0.265

-1.229

0.374

0.260

-11.255

0.957

0.161

-1.381

-5.031

0.751

Shrubs

-1.977

-0.595

1.227

-2.590

-4.418

-7.295

Lakes

1.623

-42.156

0.243

0.746

-2.912

-0.105

77.360

44.481

70.939

53.334

22.530

-14.725

0.000

0.000

2.797

10.953

-86.732

2.004

-12.178

-27.590

-4.875

-13.975

22.076

-15.762

0.000

0.000

3.734

-1.984

-3.762

9.870

-60.907

58.357

8.499

-53.552

10.158

-61.212

Other Upland Conifers (Spruce-Fir)

3.667

-45.025

-66.245

57.744

17.249

-28.484

Lowland Conifers

2.590

4.987

-16.611

-20.725

20.233

48.194

Transportation, Communication and Utilities

Grasses and Forbs

Northern Hardwoods
Central Hardwoods
Aspen-Birch
Lowland Hardwoods
Pine

172

Table IV.7. Continued.
Land use variables

LP- EUP
PLS1

WUP-WI

PLS2

PLS1

PLS2

PLS3

PLS4

Scrub Shrub Wetland

0.000

0.000

7.060

-1.100

18.661

48.389

Emergent Wetland

0.000

0.000

6.763

-0.268

6.379

11.140

0.0025

0.0011

0.0012

0.0010

0.0006

0.0005

75.87

16.86

42.39

31.85

11.82

8.58

0.96

0.90

0.63

0.73

0.68

0.70

Singular value
%Explained covariance
Correlation

Note: PLS vector coefficients multiplied x100 to facilitate interpretation. The largest loadings within vectors are indicated by
boldface.

173

Figures

Wisconsin

Michigan’S Upper Penisula (UP)

Michigan’s Lower peninsula (LP)

Figure IV.1. Land Cover/Use of populations sampled in this study. Each circle covers a radius of 454 m from the geographical
midpoint of each locality. See text for additional information.
174

Morphometric Distances LP-EUP

Figure IV.2. UPGMA dendrograms based on morphometric, genetic, and ecological distances in LP-EUP and WUP-WI lineages.

175

Figure IV.2. Continued

Morphometric Distances WUP-WI

176

Figure IV.2. Continued

Genetic Distances LP-EUP

177

Figure IV.2. Continued

Genetic Distances WUP-WI

178

Figure IV.2. Continued

Ecological Distances LP-EUP

179

Figure IV.2. Continued

Ecological Distances WUP-WI

.
180

Figure IV.3a. Correlation between temperature and shape on PLS1 for LP-EUP, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.
181

Figure IV.3b. Correlation between temperature and shape on PLS1 for WUP-WI, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.
182

Figure IV.3c. Correlation between temperature and shape on PLS2 for WUP-WI, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.
183

Figure IV.4a. Correlation between precipitation and shape on PLS1 for LP-EUP, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.
184

Figure IV.4b. Correlation between precipitation and shape on PLS1 for WUP-WI, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.

185

Figure IV.5a. Correlation between latitude and shape on PLS1 for LP-EUP, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.
186

Figure IV.5b. Correlation between longitude and shape on PLS1 for WUP-WI, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.
187

r = 0.79

Figure IV.5c. Correlation between latitude and shape on PLS2 for WUP-WI, and the corresponding deformation grid. Two-letter
acronyms are as in Table A.IV.1.
188

r = 0.96

Figure IV.6a. Correlation between land cover/use and shape on PLS1 for LP-EUP, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.
189

Figure IV.6b. Correlation between land cover/use and shape on PLS2 for LP-EUP, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.

190

r = 0.63

Figure IV.6c. Correlation between land cover/use and shape on PLS1 for WUP-WI, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.
191

r = 0.73

Figure IV.6d. Correlation between land cover/use and shape on PLS2 for WUP-WI, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.
192

r = 0.68

Figure IV.6e. Correlation between land cover/use and shape on PLS3 for WUP-WI, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.
193

r = 0.70

Figure IV.6f. Correlation between land cover/use and shape on PLS4 for WUP-WI, and the corresponding deformation grid. Twoletter acronyms are as in Table A.IV.1.

194

APPENDIX

195

Table A.IV.1. Temperature and precipitation variables and coordinates of the geographical midpoint for all populations in this study.
Temperatures are in Fahrenheit degrees; snow and rain precipitations in mm2. In Timing, O and N stand for old and new, respectively.
l is Simpson’s Index of Dominance.
Population
Adams
Ashland
Bayfield
Buffalo
Charlevoix
Cheboygan
Chippewa
Crawford
Delta
Emmet
Gogebic
Iron
Marathon
Marinette
Marquette
Menominee
Oneida
Ontonagon
Otsego
Outagamie
Portage
Schoolcraft
Taylor
Wood

ID
AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

Timing
O
N
N
O
O
O
N
O
N
O
N
N
O
O
N
O
N
N
O
O
O
N
O
O

Lat
43.9012
46.1568
46.3705
44.5892
45.2615
45.5634
46.4320
44.7888
46.0788
45.5431
46.3437
46.3579
44.8175
45.5900
46.8736
45.1580
45.8709
46.7364
45.0549
44.3961
44.5412
46.0886
45.3127
44.3785

Long
-89.8676
-90.4755
-91.3376
-91.8022
-84.8037
-84.6637
-84.7266
-84.7340
-86.8345
-85.0478
-89.3678
-90.3774
-89.7509
-87.9464
-87.8958
-87.7000
-89.6527
-89.7394
-84.5688
-88.6663
-89.5654
-86.2278
-90.6268
-90.1090

Annual T
44.50
40.90
41.60
45.70
44.80
42.50
39.00
41.90
42.00
43.70
39.00
39.60
43.60
41.30
43.00
44.60
39.20
42.40
43.30
44.20
43.10
43.20
41.10
44.30

196

Winter T
18.30
13.77
14.53
18.23
22.93
19.87
15.83
18.50
18.83
22.10
12.83
13.43
16.90
15.63
20.87
19.83
12.97
18.40
20.00
18.10
16.43
19.23
14.13
18.20

Spring T
44.63
41.20
41.30
46.43
42.67
39.03
36.40
39.87
38.67
40.13
38.03
38.80
43.67
40.73
40.50
43.00
38.43
41.17
41.60
43.90
43.07
41.40
41.83
44.77

Summer T
68.10
65.10
66.00
69.93
65.70
64.67
61.00
64.30
64.30
64.73
63.00
63.13
67.83
64.90
64.13
68.00
63.13
64.67
65.40
68.10
67.43
65.67
64.73
67.53

Fall T
46.87
43.33
44.57
48.10
48.03
46.50
42.90
44.80
46.17
47.67
42.30
43.03
46.00
43.77
46.40
47.67
42.23
45.53
46.07
46.80
45.40
46.47
43.80
46.87

Table A.IV.1. Continued.
Population
Adams
Ashland
Bayfield
Buffalo
Charlevoix
Cheboygan
Chippewa
Crawford
Delta
Emmet
Gogebic
Iron
Marathon
Marinette
Marquette
Menominee
Oneida
Ontonagon
Otsego
Outagamie
Portage
Schoolcraft
Taylor
Wood

ID
AD
AS
BY
BF
CX
CH
CW
CR
DE
EM
GO
IR
MH
MR
MA
MN
OD
ON
OT
OU
PO
SC
TA
WO

Timing
O
N
N
O
O
O
N
O
N
O
N
N
O
O
N
O
N
N
O
O
O
N
O
O

Snow
50.10
41.20
68.10
47.80
106.50
90.00
187.80
104.70
46.70
120.80
98.30
139.40
59.10
53.10
120.90
53.70
58.00
187.80
149.60
42.50
44.50
97.90
56.90
41.30

Rain
33.28
32.09
34.25
33.56
32.14
29.56
35.80
33.42
28.53
31.15
30.34
34.40
33.36
29.83
30.03
32.40
32.09
33.56
36.59
30.82
32.13
30.76
32.81
31.92

l
0.3068
0.7482
0.5013
0.3938
0.9761
0.3664
0.4908
0.3550
0.7972
0.5026
0.2533
0.3088
0.2216
0.6876
0.8937
0.3139
0.7490
0.9672
0.3178
0.6242
0.2480
0.5195
0.3712
0.1808

197

Table A.IV.2. Matrix of pairwise ecological distances, (1-IM), were IM is Morista’s Index of community similarity.
Adams
Ashland
Bayfield
Buffalo
Charlevoix
Cheboygan
Chippewa
Crawford
Delta
Emmet
Gogebic
Iron
Marathon
Marinette
Marquette
Menominee
Oneida
Ontonagon
Otsego
Outagamie
Portage
Schoolcraft
Taylor
Wood

0
0.45
0.46
0.82
0.48
0.52
0.13
0.28
0.46
0.45
0.64
0.53
0.81
0.98
1.00
0.34
0.98
0.48
0.46
0.94
0.49
0.36
0.45
0.29

0
0.07
0.97
0.01
0.22
0.17
0.87
0.00
0.07
0.49
0.29
0.76
1.00
0.91
0.15
0.87
0.01
0.19
0.95
0.28
0.99
0.21
0.71

0
1.00
0.13
0.12
0.21
0.87
0.09
0.07
0.28
0.22
0.76
1.00
0.63
0.17
0.59
0.13
0.21
1.00
0.26
0.99
0.25
0.70

0
1.00
1.00
1.00
1.00
0.99
1.00
0.88
0.86
0.94
1.00
0.99
0.90
0.97
0.99
0.93
0.49
0.86
1.00
0.70
0.99

0
0.29
0.20
0.88
0.01
0.11
0.58
0.36
0.78
1.00
0.99
0.21
0.95
0.00
0.24
1.00
0.34
0.99
0.27
0.74

0
0.32
0.89
0.24
0.20
0.30
0.29
0.78
1.00
0.63
0.24
0.57
0.29
0.30
1.00
0.33
0.99
0.33
0.67

198

0
0.41
0.18
0.20
0.57
0.37
0.76
0.98
0.97
0.17
0.93
0.20
0.27
1.00
0.34
0.50
0.32
0.38

0
0.87
0.78
0.85
0.72
0.92
0.97
1.00
0.74
0.99
0.88
0.79
1.00
0.86
0.06
0.89
0.30

0
0.08
0.51
0.31
0.77
1.00
0.94
0.15
0.90
0.01
0.19
0.94
0.29
0.99
0.21
0.72

0
0.36
0.10
0.75
1.00
0.84
0.16
0.80
0.11
0.20
1.00
0.26
0.87
0.24
0.70

0
0.23
0.86
1.00
0.37
0.44
0.33
0.59
0.46
0.73
0.46
0.90
0.38
0.82

0
0.79
1.00
0.73
0.30
0.69
0.36
0.33
0.93
0.35
0.77
0.34
0.74

Table A.IV.2. Continued
Adams
Ashland
Bayfield
Buffalo
Charlevoix
Cheboygan
Chippewa
Crawford
Delta
Emmet
Gogebic
Iron
Marathon
Marinette
Marquette
Menominee
Oneida
Ontonagon
Otsego
Outagamie
Portage
Schoolcraft
Taylor
Wood

0
0.99
1.00
0.66
0.99
0.78
0.62
1.00
0.73
0.97
0.78
0.61

0
1.00
0.99
1.00
1.00
0.98
1.00
1.00
0.97
1.00
0.94

0
0.99
0.01
1.00
0.99
0.96
0.94
1.00
0.94
1.00

0
0.96
0.20
0.07
0.75
0.16
0.88
0.09
0.55

0
0.96
0.96
0.96
0.91
1.00
0.91
0.98

199

0
0.24
1.00
0.34
0.99
0.27
0.74

0
0.73
0.19
0.99
0.12
0.61

0
0.80
1.00
0.44
1.00

0
0.99
0.21
0.61

0
0.99
0.41

0
0.72

0

Table A.IV.3. Matrix of pairwise geographic distances, in km.
Adams
Ashland
Bayfield
Buffalo
Charlevoix
Cheboygan
Chippewa
Crawford
Delta
Emmet
Gogebic
Iron
Marathon
Marinette
Marquette
Menominee
Oneida
Ontonagon
Otsego
Outagamie
Portage
Schoolcraft
Taylor
Wood

0.0
255.3
297.8
172.0
428.5
450.6
491.4
419.9
339.8
422.2
274.4
276.1
102.3
241.4
364.6
221.5
219.7
315.4
439.4
110.5
75.1
375.5
168.1
56.5

0.0
70.4
202.8
451.4
454.8
442.7
472.7
280.8
425.8
87.7
23.6
159.3
205.7
212.9
242.6
71.0
85.6
475.4
241.6
193.2
327.4
94.6
199.8

0.0
201.4
521.0
523.4
506.8
543.0
347.9
494.7
151.2
73.7
212.3
276.0
268.7
312.7
141.2
128.8
545.3
302.8
245.9
394.2
129.9
241.4

0.0
555.9
570.7
587.9
559.0
422.1
540.8
272.2
225.9
164.1
322.5
395.5
329.3
220.5
287.6
572.7
249.7
177.3
466.4
122.6
136.3

0.0
35.3
130.3
52.8
182.1
36.7
373.6
448.8
391.8
247.9
298.4
227.2
383.4
414.9
29.4
319.4
383.4
143.9
455.5
429.7

0.0
96.7
86.3
177.7
30.0
373.8
450.3
407.1
255.5
288.2
241.5
388.8
412.1
57.0
340.5
401.4
134.5
466.0
448.1

200

0.0
182.7
166.8
101.9
356.1
433.3
429.9
265.7
246.8
270.5
384.5
384.5
153.6
381.8
431.7
121.6
473.3
478.1

0.0
217.8
87.4
399.9
472.5
395.8
267.0
337.2
236.9
402.8
444.5
32.3
314.4
383.0
185.7
466.5
428.0

0.0
150.7
197.1
274.3
267.2
101.9
120.1
122.5
219.0
234.4
209.9
235.7
273.5
46.8
306.6
318.5

0.0
345.6
421.8
377.3
225.7
264.4
211.6
359.4
385.0
66.0
311.9
372.0
109.7
436.0
418.7

0.0
77.5
172.3
138.2
126.9
184.7
57.0
52.1
399.2
223.4
201.0
243.2
150.5
226.1

0.0
178.1
206.3
198.0
246.9
77.8
64.4
473.6
255.8
211.7
320.6
117.8
221.1

Table A.IV.3. Continued.
Adams
Ashland
Bayfield
Buffalo
Charlevoix
Cheboygan
Chippewa
Crawford
Delta
Emmet
Gogebic
Iron
Marathon
Marinette
Marquette
Menominee
Oneida
Ontonagon
Otsego
Outagamie
Portage
Schoolcraft
Taylor
Wood

0.0
165.4
270.0
165.7
117.4
213.4
408.7
97.8
34.0
309.0
88.1
56.5

0.0
142.8
51.7
136.1
187.9
270.7
144.3
172.5
144.2
211.3
217.0

0.0
191.4
174.9
141.1
327.1
281.9
289.9
154.7
272.9
326.4

0.0
171.6
235.9
246.0
114.0
162.3
154.3
229.8
209.0

0.0
96.5
406.7
181.3
148.0
265.7
98.0
169.8

0.0
441.6
273.3
244.5
278.7
172.5
263.8

201

0.0
331.9
398.3
172.9
475.5
444.1

0.0
73.2
268.1
185.1
114.7

0.0
312.6
119.7
46.8

0.0
352.3
358.4

0.0
111.6

0.0

BIBLIOGRAPHY

202

BIBLIOGRAPHY

Anapol F, Lee S (1994) Morphological adaptation to diet in Platyrrhine primates. American
Journal of Physical Anthropology, 94, 239-261.

Anderson JR, Hardy EE, Roach JT, Witmer RE (1976) A Land Use and Land Cover
Classification System for Use with Remote Sensor Data. (US Geological Survey
professional paper ; 964) (ed. Interior Do), p. iii + 29. United States Government Printing
Office, Washington.

Ashton KG, Tracy MC, de Queiroz A (2000) Is Bergmann's rule valid for mammals? American
Naturalist, 156, 390-415.

Atchley WR (1993) Genetic and Developmental Aspects of Variability in the Mammalian
Mandible. In: The Skull (eds. Hanken J, Hall BK), pp. 207-247. University of Chicago
Press, Chicago. Illinois.

Atchley WR, Hall BK (1991) A model for development and evolution of complex morphological
structures. Biological Reviews of the Cambridge Philosophical Society, 66, 101-157.

Bradshaw WE, Holzapfel CM (2006) Climate change - Evolutionary response to rapid climate
change. Science, 312, 1477-1478.

Bronson FH (2009) Climate change and seasonal reproduction in mammals. Philosophical
Transactions of the Royal Society B-Biological Sciences, 364, 3331-3340.

Brower JE, Zar JH, Ende CNv (1990) Field and laboratory methods for general ecology, 3rd
edn. W.C. Brown Publishers, Dubuque, Iowa.

Burt WH (1940) Territorial behavior and populations of some small mammals in southern
Michigan University of Michigan press, Ann Arbor, Michigan.

Cardini A, Jansson AU, Elton S (2007) A geometric morphometric approach to the study of
ecogeographical and clinal variation in vervet monkeys. Journal of Biogeography, 34,
1663-1678.
203

Chevin LM, Lande R, Mace GM (2010) Adaptation, Plasticity, and Extinction in a Changing
Environment: Towards a Predictive Theory. Plos Biology, 8, e1000357.

Corruccini RS, Beecher RM (1982) Occlusal variation related to soft diet in a non-human
primate. Science, 218, 74-76.

Excoffier L, Laval G, Schneider S (2005) Arlequin (version 3.0): An integrated software package
for population genetics data analysis. Evolutionary Bioinformatics, 1, 47-50.

Gannon WL, Sikes RS, Comm ACU (2007) Guidelines of the American Society of
Mammalogists for the use of wild mammals in research. Journal of Mammalogy, 88, 809823.

Gibson G, Wagner G (2000) Canalization in evolutionary genetics: a stabilizing theory?
Bioessays, 22, 372-380.

Goheen JR, Swihart RK, Robins JH (2003) The anatomy of a range expansion: changes in
cranial morphology and rates of energy extraction for North American red squirrels from
different latitudes. Oikos, 102, 33-44.

Heath ME (1984) The effects of rearing-temperature on body conformation and organ size in
young-pigs. Comparative Biochemistry and Physiology B-Biochemistry & Molecular
Biology, 77, 63-72.

Hedrick PW, Ginevan ME, Ewing EP (1976) Genetic-polymorphism in heterogeneous
environments. Annual Review of Ecology and Systematics, 7, 1-32.

Hughes L (2000) Biological consequences of global warming: is the signal already apparent?
Trends in Ecology & Evolution, 15, 56-61.

Ivan JS, Swihart RK (2000) Selection of mast by granivorous rodents of the central hardwood
forest region. Journal of Mammalogy, 81, 549-562.

Jacquot JJ, Vessey SH (1995) Influence of the natal environment on dispersal of white-footed
mice. Behavioral Ecology and Sociobiology, 37, 407-412.

204

Johnson CN (1986) Sex-biased philopatry and dispersal in mammals. Oecologia, 69, 626-627.

Lambert JE (2009) Summary to the Symposium Issue: Primate Fallback Strategies as Adaptive
Phenotypic Plasticity-Scale, Pattern, and Process. American Journal of Physical
Anthropology, 140, 759-766.

Lande R (2009) Adaptation to an extraordinary environment by evolution of phenotypic
plasticity and genetic assimilation. Journal of Evolutionary Biology, 22, 1435-1446.

Lewis CE, Clark TW, Derting TL (2001) Food selection by the white-footed mouse (Peromyscus
leucopus) on the basis of energy and protein contents. Canadian Journal of Zoology, 79,
562-568.

Maier TJ (2002) Long-distance movements by female White-footed Mice, Peromyscus leucopus,
in extensive mixed-wood forest. Canadian Field Naturalist, 116, 108-111.

Maki K, Nishioka T, Shioiri E, Takahashi T, Kimura M (2002) Effects of dietary consistency on
the mandible of rats at the growth stage: Computed X-ray densitometric and
cephalometric analysis. Angle Orthodontist, 72, 468-475.

Marshall AJ, Wrangham RW (2007) Evolutionary consequences of fallback foods. International
Journal of Primatology, 28, 1218-1235.

Michaux J, Chevret P, Renaud S (2007) Morphological diversity of Old World rats and mice
(Rodentia, Muridae) mandible in relation with phylogeny and adaptation. Journal of
Zoological Systematics and Evolutionary Research, 45, 263-279.

Monteiro LR, Nogueira MR (2010) Adaptive radiations, ecological specialization, and the
evolutionary integration of complex morphological structures. Evolution, 64, 724-743.

Monteiro LR, Duarte LC, dos Reis SF (2003) Environmental correlates of geographical variation
in skull and mandible shape of the punare rat Thrichomys apereoides (Rodentia :
Echimyidae). Journal of Zoology, 261, 47-57.

205

Mossman CA, Waser PM (2001) Effects of habitat fragmentation on population genetic structure
in the white-footed mouse (Peromyscus leucopus). Canadian Journal of Zoology, 79,
285-295.

Munshi-South J, Kharchenko K (2010) Rapid, pervasive genetic differentiation of urban whitefooted mouse (Peromyscus leucopus) populations in New York City. Molecular Ecology,
19, 4242-4254.

Myers P, Lundrigan BL, Hoffman SMG, Haraminac AP, Seto SH (2009) Climate-induced
changes in the small mammal communities of the Northern Great Lakes Region. Global
Change Biology, 15, 1434-1454.

Nei M (1987) Molecular evolutionary genetics Columbia University Press, New York.
Nicholson E, Harvati K (2006) Quantitative analysis of human mandibular shape using threedimensional geometric morphometrics. American Journal of Physical Anthropology, 131,
368-383.

Nogueira MR, Peracchi AL, Monteiro LR (2009) Morphological correlates of bite force and diet
in the skull and mandible of phyllostomid bats. Functional Ecology, 23, 715-723.

Pelabon C, Hansen TF, Carter AJR, Houle D (2010) Evolution of variation and variability under
fluctuating, stabilizing, and disruptive selection. Evolution, 64, 1912-1925.

Perez SI, Monteiro LR (2009) Nonrandom factors in modern human morphological
diversification: a study of craniofacial variation in southern south american populations.
Evolution, 63, 978-993.

Pergams ORW, Ashley MV (1999) Rapid morphological change in channel island deer mice.
Evolution, 53, 1573-1581.

Pergams ORW, Lacy RC (2008) Rapid morphological and genetic change in Chicago-area
Peromyscus. Molecular Ecology, 17, 450-463.

Pergams ORW, Lawler JJ (2009) Recent and Widespread Rapid Morphological Change in
Rodents. Plos One, 4, e6452.

206

Pigot AL, Owens IPF, Orme CDL (2010) The environmental limits to geographic range
expansion in birds. Ecology Letters, 13, 705-715.

Poroshin EA, Polly PD, Wojcik JM (2010) Climate and morphological change on decadal scales:
Multiannual variation in the common shrew Sorex araneus in northeast Russia. Acta
Theriologica, 55, 193-202.

Rae TC, Vioarsdottir US, Jeffery N, Steegmann AT (2006) Developmental response to cold
stress in cranial morphology of Rattus: implications for the interpretation of climatic
adaptation in fossil hominins. Proceedings of the Royal Society B-Biological Sciences,
273, 2605-2610.

Renaud S, Auffray JC (2010) Adaptation and plasticity in insular evolution of the house mouse
mandible. Journal of Zoological Systematics and Evolutionary Research, 48, 138-150.

Renaud S, Auffray JC, de la Porte S (2010) Epigenetic effects on the mouse mandible: common
features and discrepancies in remodeling due to muscular dystrophy and response to food
consistency. Bmc Evolutionary Biology, 10, 28.

Reynolds J, Weir BS, Cockerham CC (1983) Estimation of the co-ancestry coefficient - basis for
a short-term genetic-distance. Genetics, 105, 767-779.

Rohlf FJ, Corti M (2000) Use of two-block partial least-squares to study covariation in shape.
Systematic Biology, 49, 740-753.

Rychlik L, Ramalhinho G, Polly PD (2006) Response to environmental factors and competition:
skull, mandible and tooth shapes in Polish, water shrews (Neomys, Soricidae,
Mammalia). Journal of Zoological Systematics and Evolutionary Research, 44, 339-351.

Salamin, N., R. O. Wüest, S. Lavergne, W. Thuiller, and P. B. Pearman. 2010. Assessing rapid
evolution in a changing environment. Trends in Ecology & Evolution, 25, 692-698.

Samuels JX (2009) Cranial morphology and dietary habits of rodents. Zoological Journal of the
Linnean Society, 156, 864-888.

207

Schlichting C, Pigliucci M (1998) Phenotypic evolution : a reaction norm perspective. Sinauer,
Sunderland, Massachussets.

Schlaepfer MA, Runge MC, Sherman PW (2002) Ecological and evolutionary traps. Trends in
Ecology & Evolution, 17, 474-480.

Steegman AT, Platner WS (1968) Experimental cold modification of cranio-facial morphology.
American Journal of Physical Anthropology, 28, 17-30.

Straney DO, Patton JL (1980) Phylogenetic and environmental determinants of geographicvariation of the pocket mouse Perognathus goldmani Osgood. Evolution, 34, 888-903.

Thomas CD, Bodsworth EJ, Wilson RJ, et al. (2001) Ecological and evolutionary processes at
expanding range margins. Nature, 411, 577-581.

Van der Putten WH, Macel M, Visser ME (2010) Predicting species distribution and abundance
responses to climate change: why it is essential to include biotic interactions across
trophic levels. Philosophical Transactions of the Royal Society B-Biological Sciences
365, 2025-2034.

Walther GR, Post E, Convey P, et al. (2002) Ecological responses to recent climate change.
Nature, 416, 389-395.

Weaver ME, Ingram DL (1969) Morphological changes in swine associated with environmental
temperature. Ecology, 50, 710-713.

Wolff JO, Dueser RD, Berry KS (1985) Food-habits of sympatric Peromyscus leucopus and
Peromyscus maniculatus. Journal of Mammalogy, 66, 795-798.

Zelditch ML, Swiderski DL, Sheets HD, Fink WL (2004) Geometric morphometrics for
biologists : a primer. Elsevier Academic Press, San Diego California; Amsterdam.

208