i r .2 i 3.... v5.3.3}!- Eight: fig “WW-«mm. «.2 at" . 43 x . 3. . a?» r. 2.. “a... nu» , .afl in». 19:? am...» v.2... . r a»... .. ...R.»R ., uidfifit . um.» . r 2.5.... 3. 1.. 3 5.1.. . s... .. 2d a... I» H ; trs; ! 1. ii {ling 31:3: 25.0.1.3 1.2.44.2!!! .54.. 12:33.» l...3.:.:{. . 151.13.. . n . .1... .31 .24... t 5‘3 . D‘Itvf ¢f . 1.1.1.. 1 3.. . 9. n . 41.44.4233.” :2 . .. . 4):.2»....1.: V [12‘ ifs... sfiafinfi. Agriéwaégig ..,, 4......3EJ. ,WJEEV .u. .. ".3 / 4/; .. ' /’ ”WI-b 0cm , 41‘? LIg‘ATE UiIIVERSITY CHIGAN S EXIST LANSING, MICH 48824-1048 This is to certify that the dissertation entitled SPATIAL STRUCTURE AND INVASION SUCCESS OF POPULATIONS UNDERGOING RANGE EXPANSION: EFFECTS OF DISPERSAL STRATEGY AND LANDSCAPE HETEROGENEITY presented by Genevieve Marie Nesslage has been accepted towards fulfillment of the requirements for the Ph.D. degree in Department of Fisheries and Wildlife and Program in Ecology, Evolutionary Biology and Behavior MAW/vim» Major Professor’s Signature WW Xyflaé fl / Date MSU is an Affirmative Action/Equal Opportunity Institution _ - - -._. PLACE IN RETURN Box to remove this checkout from your record.\~ To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE JAN 1 81009 010509 ‘——_———— H ‘ 2/05 chlRC/DatoDueJndd-sz “ SPATIAL STRUCTURE AND INVASION SUCCESS OF POPULATIONS UNDERGOING RANGE EXPANSION: EFFECTS OF DISPERSAL STRATEGY AND LANDSCAPE HETEROGENEITY By Genevieve Marie Nesslage A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Fisheries and Wildlife and Program in Ecology, Evolutionary Biology and Behavior 2005 ABSTRACT SPATIAL STRUCTURE AND INVASION SUCCESS OF POPULATIONS UNDERGOING RANGE EXPANSION: EFFECTS OF DISPERSAL STRATEGY AND LANDSCAPE HETEROGENEITY By Genevieve Marie Nesslage This dissertation explores the effects of dispersal strategy and landscape heterogeneity on the spatial structure and invasion success of populations undergoing range expansion. First, the spatial structure of an invasive population was evaluated using semivariogram analysis (a characterization of the dependence of data on spatial location). Twelve years (1985—1996) of historic gypsy moth (Lymantria dispar) monitoring records from the Lower Peninsula of Michigan were analyzed. Seven models were statistically fit to annual plots of semivariance against distance at two different spatial extents. The model with the lowest Akaike’s Information Criterion value was chosen to characterize overall semivariogram behavior for each year and spatial extent. At large spatial extents, three distinct patterns of semivariogram behavior were observed as the invasion progressed. This three-stage semivariance progression pattern likely represents the invasion stages of establishment, expansion, and stabilization and may be indicators of important mechanistic changes in invasion dynamics. Second, the relationship between landscape heterogeneity and invasion success of the gypsy moth in Michigan was analyzed by regressing invasion success on seven measures of landscape fragmentation and by identifying potential threshold responses (sudden, nonlinear changes) in invasion success with habitat loss. Invasion success thresholds were then compared with that of seven landscape structure metrics. Both patch- and gap-based landscape metrics exhibited threshold responses. However, invasion success increased linearly with increasing proportion of habitat, indicating that the negative effects of habitat fragmentation do not compound that of habitat loss. The absence of detectable thresholds in invasion success suggests that aerial dispersers like the invasive gypsy moth may not be as strongly affected by habitat fragmentation as previously thought. Lastly, semivariogram analysis was used to evaluate patterns of spatial autocorrelation generated by simple and stratified dispersal models. Dispersal was simulated across a homogeneous landscape composed of 100% habitat and a real, heterogeneous landscape composed of 54% habitat. Population growth rate and stratified dispersal distributions were varied as well. Semivariograms were calculated from model simulations at every other time step and fit with four different models to characterize changes in spatial structure of the population over time. Both simple and stratified dispersal modeled across real, heterogeneous landscapes generated semivariogram progression patterns similar to that of gypsy moths. Models of dispersal across 100% habitat landscapes and models of stratified dispersal with a long-tailed dispersal distribution across real landscapes generated semivariance progression patterns dissimilar to that of gypsy moths. Thus, landscape heterogeneity (and not dispersal strategy) was a major driving force behind semivariance patterns observed during the gypsy moth invasion of Michigan. ACKNOWLEDGEMENTS I sincerely thank all those who have helped me in the pursuit of my doctorate. In particular, I would like to thank my advisor Dr. Brian Maurer for his encouragement of my research and career. I have been privileged to work with one of the most outstanding ecologists of our time. In addition, I owe the success of my dissertation to Dr. Stuart Gage who provided me with both the data and the knowledge of dispersal ecology necessary to conduct this research project. I am deeply indebted to Dr. Michael Jones and Dr. Richard Kobe for their thoughtful and constructive reviews of my research ideas and manuscripts. I also thank Dr. Michael Wilberg for providing me with endless statistical advice and personal support over the last five years. His exceptional intelligence and kindness have inspired me to be a better scientist and a better person. I would like to thank all those who offered me advice, support, and career- advancing Opportunities over the last five years: Michael Donovan of the Michigan Department of Natural Resources; Drs. William Taylor, Thomas Coon, Shawn Riley, Henry Campa, and Robert Batie in the Department of Fisheries and Wildlife; my lab mates Jennifer Skillen, Katherine Kahl, Anne Axel, Robert Lanfear, Hui-Yu Wang, and Smruti Damania; and my parents Lee and Arthur Nesslage. Funding for this project was provided in part by The Graduate School and the Ecology, Evolutionary Biology, and Behavior Program at Michigan State University. iv TABLE OF CONTENTS LIST OF TABLES ............................................................................................................ vii LIST OF FIGURES ......................................................................................................... viii INTRODUCTION .............................................................................................................. 1 CHAPTER 1 CHARACTERIZING THE INVASION PROCESS USING PATTERNS OF SEMIVARIANCE: A CASE STUDY OF THE GYPYSY MOTH IN MICHIGAN ........ 5 Introduction .................................................................................................................. 5 Methods ..................................................................................................................... 11 Results ........................................................................................................................ 15 Variogram behavior in Michigan’s Lower Peninsula .................................... 15 Variogram behavior at the ecoregion scale ................................................... 16 Discussion .................................................................................................................. 20 Variogram behavior in Michigan’s Lower Peninsula .................................... 20 Variogram behavior at ecoregion scale ......................................................... 21 Timing of shifts in variogram behavior ......................................................... 23 Use of variograms in invasion analysis and interpretation ............................ 25 CHAPTER 2 THE RELATIONSHIP BETWEEN INVASION SUCCESS AND LANDSCAPE HETEROGENEITY: A CASE STUDY OF THE GYPSY MOTH IN MICHIGAN ...... 35 Introduction ................................................................................................................ 35 Methods ..................................................................................................................... 39 Invasion data .................................................................................................. 39 Land cover data .............................................................................................. 39 Scale of analysis ............................................................................................ 4O Characterizing the relationship between invasion success and landscape structure ......................................................................................................... 41 Results ........................................................................................................................ 43 Relationship between invasion success and habitat loss ............................... 43 Relationship between landscape metrics and habitat loss ............................. 43 Relationship between invasion success and habitat fragmentation ............... 45 Discussion .................................................................................................................. 46 Effects of habitat loss and fragmentation on invasion success ...................... 46 Effect of scale of analysis and land cover map .............................................. 49 Management implications .............................................................................. 51 CHAPTER 3 PATTERNS OF SEMIVARIANCE GENERATED BY SIMPLE AND STRATIFIED DISPERSAL ..................................................................................................................... 68 Introduction ................................................................................................................ 68 Methods ..................................................................................................................... 70 Model structure .............................................................................................. 7O Characterizing changes in semivariance ........................................................ 72 Results ........................................................................................................................ 74 Discussion .................................................................................................................. 77 CONCLUSIONS ............................................................................................................... 96 Detection of emergent properties ............................................................................... 96 Semivariance analysis ................................................................................................ 98 Future directions ........................................................................................................ 98 LITERATURE CITED ................................................................................................... 100 vi LIST OF TABLES Table 1.1. Modeling results from variogram analysis of annual gypsy moth (Lymantria dispar) trap catch in the Lower Peninsula (LP) and 10 ecoregions of Michigan, 1985-1996. The best model (model with lowest Akaike Information Criterion value) of seven models fit to annual variograms of gypsy moth trap catch is listed as “Type”. Note that the exponential (Expon.), Gaussian (Gaus.), and spherical (Spher.) model names were abbreviated. Model shape (Shape) describes the overall shape of the variogram being fit (B = bell or wave shaped, L = linear, and U = untrended variogram with typical increase in semivariance to an asymptote). ........ 28 Table 2.1. Summary of critical thresholds in landscape metrics with habitat loss at three scales of analysis (i.e. analysis window sizes of 75, 45, and 15 km on a side) using two land cover maps. Dashes indicate no threshold response was detected. ........... 54 Table 3.1. Description of twelve model scenarios, including dispersal strategy modeled, landscape type, population growth rate (r), and dispersal distributions used in model simulations. Distribution 1 is a gamma distribution with rate parameter of 0.5 and shape parameter of 1.2; Distribution 2 is a long-tailed gamma distribution with rate parameter of 0.043 and Shape parameter of 0.5. ....................................................... 82 vii LIST OF FIGURES Figure 1.1. Map of Michigan’s Lower Peninsula showing its location in North America and ecoregion boundaries used for variogram analysis of gypsy moth (Lymantria dispar) catch between 1985 and 1996 ....................................................................... 30 Figure 1.2. Annual semivariograms of male gypsy moth (Lymantria dispar) trap catch in the Lower Peninsula of Michigan, 1985—1996. Circles represent the original data and Open squares represent predicted values generated by the best—fit models. Predicted values were generated by either the cubic (1985—1988), wave (1989), spherical (1990), power (1991—1995), or exponential (1996) models. Note the shift in variogram shape from bell—shaped (1985—1989) to linear (1990—1995) to an untrended pattern of increasing semivariance and leveling off at an asymptote. ..... 31 Figure 1.3. Total annual male gypsy moth (Lymantria dispar) catch in the Lower Peninsula of Michigan between 1985 and 1996. Arrows indicate years in which variogram patterns changed from bell—shaped to linear (right—facing arrow) and from linear to typical increase to asymptote (left—facing arrow). Note that inspection of the time series alone (without consideration of Spatial location of data points) would indicate that the invasion had stabilized by 1990. ......................................... 33 Figure 1.4. Total annual male gypsy moth (Lymantria dispar) trap catch in each often ecoregions of the Lower Peninsula of Michigan between 1985 and 1996. .............. 34 Figure 2.1. Range in years to colonization among traps located in the same analysis window declined with increasing proportion of habitat. Results are reported using Map 1 for each scale of analysis: a) 75 km, b) 45km, and c) 15km. ........................ 55 Figure 2.2. Thresholds in landscape metrics with increasing proportion of habitat calculated using land cover Map 1 and a 75 km analysis window. .......................... 56 Figure 2.3. Thresholds in landscape metrics with increasing proportion of habitat calculated using land cover Map 1 and a 45 km analysis window. .......................... 58 Figure 2.4. Thresholds in landscape metrics with increasing proportion of habitat calculated using land cover Map 1 and a 15 km analysis window. .......................... 60 Figure 2.5. Relationship between range in years to colonization among traps located in the same analysis window and landscape metrics calculated using land cover Map 1 and a 75 km analysis window. .................................................................................. 62 Figure 2.6. Relationship between range in years to colonization among traps located in the same analysis window and landscape metrics calculated using land cover Map 1 and a 45 km analysis window. .................................................................................. 64 viii Figure 2.7. Relationship between range in years to colonization among traps located in the same analysis window and landscape metrics calculated using land cover Map 1 and a 15 km analysis window. .................................................................................. 66 Figure 3.1. Two dispersal distributions used to determine the proportion of dispersers moving a given dispersal distance (number of cells in the landscape grid). The solid line represents a gamma distribution with shape = 1.2, and rate = 0.5. The dashed line represents a long-tailed gamma distribution with shape = 0.043, and rate = 0.5. ................................................................................................................................... 83 Figure 3.2. Change in variogram behavior for model scenario realsim1.7 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 84 Figure 3.3. Change in variogram behavior for model scenario realsim2.1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 85 Figure 3.4. Change in variogram behavior for model scenario realstrat1.7Dl characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. .................. 86 Figure 3.5. Change in variogram behavior for model scenario realstrat2.lD1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time Step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. .................. 87 ix Figure 3.6. Change in variogram behavior for model scenario realstratl .7D2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. .................. 88 Figure 3.7. Change in variogram behavior for model scenario realstrat2.1D2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (Open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios ................... 89 Figure 3.8. Change in variogram behavior for model scenario sim1.7 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 90 Figure 3.9. Change in variogram behavior for model scenario sim2.1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 91 Figure 3.10. Change in variogram behavior for model scenario stratl .7D2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 92 Figure 3.11. Change in variogram behavior for model scenario strat2.lD1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios ................................................. 93 Figure 3.12. Change in variogram behavior for model scenario stratl .7D1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 94 Figure 3.13. Change in variogram behavior for model scenario strat2.lD2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. ................................................ 95 xi INTRODUCTION Range expansion is a complex process that is affected by the interaction of many factors, including short and long distance dispersal (Clark et a1. 2001 , Trakhtenbrot et a1. 2005), landscape heterogeneity (Hengeveld and Van den Bosch 1997, With 2002), interspecific interactions (Richardson et al. 2000a, French and Travis 2001), and evolution (Carroll and Dingle 1996, Travis and Dytham 2002). Recent work in the areas of dispersal and invasion ecology indicates that environmental heterogeneity may be one of the most influential factors affecting dispersal success of invasive species (With 2002, Hastings et al. 2005). However, these studies do not synthesize all the spatial information potentially available to ecologists. In particular, changes in the spatial autocorrelation structure of populations expanding across heterogeneous landscapes have not been thoroughly examined. Brodie et a1. (1995) first noted changes in spatial autocorrelation patterns as a result of local dispersal in a study of the post—fire re— invasion of a Populus balsamifera clone. Tracking such changes in the Spatial structure may be a reliable indicator of important changes in population growth and dispersal patterns during the course of an invasion. However, changes in spatial population structure have not been identified in other environments or for species with difi‘erent dispersal behaviors. This dissertation contains three chapters that, respectively, address each of the following questions: 1. Can important spatial changes in population dynamics over the course of an invasion be characterized by geostatistics such as semivariance? 2. How is the movement of an invasion wavefront related to measures of landscape heterogeneity? 3. How do dispersal strategy and landscape heterogeneity affect the generation of semivariance progression patterns during range expansion? In my first chapter, I identified new applications of geostatistics, specifically semivariance, to the problem of detecting invasion stage transitions. The term “invasion stages” refers to changes in population growth and dispersal patterns during the course of an invasion that may be important indicators of mechanistic changes in the invasion process. The goal of this chapter was to make and test predictions about the types of semivariogram models that would best fit expected patterns of semivariogram behavior at each stage of the invasion process. My objectives were to compare variogram behavior at different spatial extents over the course of an invasion, and to identify potential benefits and drawbacks of using semivariogram analysis as a tool in analyzing the progress of biological invasions. To meet my Objectives, I analyzed monitoring records from the gypsy moth (Lymantria dispar) invasion of Michigan (Yang et al. 1998). The gypsy moth is an exotic insect native to Europe that has caused extensive forest defoliation across much of the eastern United States Since its introduction to Massachusetts in 1868 (Elkinton and Liebhold 1990). The first breeding population of gypsy moths in Michigan was the result of an independent introduction in 1954 (O‘Dell 1955, Hanna 1981). Gypsy moths quickly spread across the state, reaching most areas of the Lower Peninsula by the early 19905 and the western Upper Peninsula by the late 19903 (Lele et al. 1998, Yang et al. 1998). The gypsy moth monitoring program in Michigan involved the placement of pheromone-scented traps in sections eight and 26 of every township in Michigan (Yang et al. 1998). At each trap, the total annual male moth catch and the trap location in UTM coordinates were recorded. Traps were run between 1985 and 1996 in the Lower Peninsula (LP) and between 1986 and 1996 in the Upper Peninsula. The Michigan gypsy moth monitoring program records were ideal for analyzing changes in Spatial population structure during range expansion because they contained spatially explicit data collected at the same locations across a wide area for a long period of time. My first chapter analyses characterized spatial changes in population structure of the gypsy moth population over time, but those analyses did not consider the effects of habitat on the range expansion process. Therefore, in my second chapter I investigated the relationship between invasion success and measures Of landscape heterogeneity. To date, landscape ecologists have used simple models of biological invasions across computer-generated landscapes to make predictions about how the dispersal success is affected by landscape heterogeneity (With and Crist 1995, With and King 1999a, Collingham and Huntley 2000). Such models indicate that dispersal success declines with habitat loss, but a precipitous decline in dispersal success may be observed if habitat patches begin to break down and the negative effects of such habitat fragmentation are additive (With and Crist 1995, With and King 1999a, Collingham and Huntley 2000). To characterize the effects of landscape heterogeneity on invasion success, threshold responses (sudden, nonlinear changes) in dispersal success with habitat loss are compared to thresholds in landscape structure measurements. Landscape measurements that exhibit thresholds most Similar to thresholds in dispersal success are thought to be the best candidate measurements for predicting the effects of landscape structure on dispersal success (With 2002). The goal Of this study was to characterize the relationship between invasion success of a real species and the structure of a real landscape. My objectives were to compare thresholds in gypsy moth invasion success and the structure of Michigan’s landscape with thresholds produced by dispersal models run on computer- generated landscapes. Also, I assessed the effect of scale on my results by repeating invasion and landscape structure calculations at three scales of analysis using two different classified land cover maps of Michigan. The goal of my third chapter was to identify the underlying process or processes responsible for the generation of semivariance progression patterns Observed in my first chapter. I hypothesized that differences in either dispersal strategy (simple vs. stratified) or landscape heterogeneity would affect the spatial structure of a population as it underwent range expansion. To test this hypothesis, I assessed the ability of two dispersal models to generate distinct changes in semivariance patterns over time and compared model results with gypsy moths semivariance patterns (Chapter 1). My Objectives were to explore the effects of landscape heterogeneity, population growth rate, and the amount of long-distance dispersal on semivariance patterns. CHAPTER 1 CHARACTERIZIN G THE INVASION PROCESS USING PATTERNS OF SEMIVARIANCE: A CASE STUDY OF THE GYPYSY MOTH IN MICHIGAN Introduction Predicting the course of biological invasions has been an important question since the emergence of invasion ecology almost 50 years ago (Reichard and White 2003). However, a broad predictive framework for the invasion process has yet to be proposed (Williamson 1999). One problem with accurately predicting the course of an invasion is that data are not typically collected at a temporal and spatial extent large enough to document and analyze the ecological processes of interest, namely dispersal and range expansion. Even when the scale of data collection matches the scale of the invasion process, prediction is complicated by the fact that overall patterns of population growth and spread change as an invasion progresses (Townsend 1991). Changes in growth and dispersal patterns during the course of an invasion are often referred to as "invasion stages" (Hengeveld 1989, Brown 1993, Cousens and Mortimer 1995, Shigesada et al. 1995, Hastings 1996b). Typically, three distinct stages are recognized following “introduction”, or the initial arrival of an invasive species in a new area. The first of these three stages is called “establishment” (or “colonization”) and represents the formation of a self-sustaining population through successful survival and reproduction (Anderson and May 1986, Cronk and Fuller 1995, Shigesada and Kawasaki 1997, Richardson et al. 2000b). The second stage is called “expansion” (or “spread”), and represents the dispersal and potential formation of new self—sustaining populations (Hastings 1996a). A final stage called “stabilization” follows expansion and represents the stage at which an invader ceases to spread and becomes integrated into its new environment (Cronk and Fuller 1995, Vermeij 1996). Reliably identifying invasion stages is a critical first step in breaking down the invasion process and understanding how the effects of different ecological factors (e.g. landscape heterogeneity) change throughout the course of an invasion (With 2002). The most common method for identifying the stages of invasion in empirical data is analyzing temporal trends in population growth. The population trajectory of an invader is theorized to undergo an exponential increase until environmental factors become limiting (Ashton and Mitchell 1989). The stages of population growth are usually characterized by an initial period of low density during which numbers increase exponentially (Mack 1985, Cousens and Mortimer 1995). This stage is followed by a period of density dependent population growth that ends when the upper density limit of the population is reached. After these stages are complete, invaders may exhibit decline, damped oscillations, stable oscillations, a stable equilibrium, multiple equilibria, or chaos (May 1974, Brown 1993). Mack (1985) suggested that the three stages of population growth correspond with invasion stages (e. g. establishment, expansion, and stabilization). However, there are several problems with using density dependent phase transitions to define invasion stages. Population increase with and without spatial spread (establishment vs. expansion) cannot be distinguished without corresponding spatial information. Also, population growth during the expansion stage may not be density limited (Mack 1985). Kowarik (1995) proposed that a time lag of slow, constant population growth signifies a clear transition between the conceptual stages of establishment and expansion. However, low initial population densities can make detection of early invasion stages difficult and apparent time lags should be interpreted with care. Cousens and Mortimer (1995) caution that it is extremely difficult to distinguish between a time lag and the early stages of exponential growth with a constant, proportional rate of increase. At the beginning of an invasion, low initial population size or low intrinsic growth rate may give the impression of a time lag when there is none (Case 2000). Additional problems arise when identifying the transition between expansion and stabilization. Curve—fitting an incomplete time series may lead to false conclusions about the status of an invasion if sampling error or year—to—year variation are significant (Case 2000). Given the difficulties of identifying invasion stages in time series data, invasion dynamics may be too variable among organisms for a common set of invasion stages to be reliably identified. An alternative explanation may be that stages of invasion could be more clearly characterized when spatial dynamics are explicitly considered. If this is true, a more spatially explicit approach is needed to identify the stages Of an invasion. Ecological data are typically spatially autocorrelated (Rossi et al. 1992, Legendre 1993) so that data collected at locations close to one another are more similar than data collected at locations farther apart (Cliff and 0rd 1981, Griffith 1987). Analysis of spatially autocorrelated data collected from permanent locations over a large area relative to the area occupied by an invasive species may provide reliable indications of transitions between invasion stages. Spatial autocorrelation often violates crucial assumptions of many traditional statistical analyses because sample values are not independent (Legendre and Legendre 1998). However, patterns of spatial autocorrelation result from important ecological processes and should be considered a fundamental component of the ecology of a Species (Legendre 1993). Careful inspection of these patterns may be critical in identifying the stages of invasion. For example, detailed information on patterns of spatial autocorrelation has been used to identify the stages of post-fire re— invasion of a Populus balsamifera clone using correlograms, or plots of Moran’s I (an index of spatial autocorrelation) against distance (Brodie et al. 1995). Three distinct patterns of spatial autocorrelation were observed for each correlogram representing a different age class, suggesting three separate stages of clonal development: post-fire colonization, consolidation, and directional expansion (Brodie et al. 1995). These stages may be roughly analogous to the invasion stages of introduction, establishment, and expansion. A related method for identifying changes in spatial autocorrelation patterns over time is semivariance analysis. Semivariance is a spatial statistic used to characterize dependence of sample values across space. The graph of semivariance against distance class is called a semivariogram (or simply variogram) and depicts the average dissimilarity between values based on separation distance between sampling points. The standardized semivariance is equal to 1 minus the correlation when the population mean and variance remain constant across the study area (Rossi et a1. 1992). Semivariogram models are often used to interpolate between sample locations within a study area (e. g. kriging, Burgess and Webster 1980), so variograms are most often treated as a means to an end. Variograms by themselves, however, provide important information about the invasion process that other methods such as time series analysis cannot (Rossi et al. 1992). For instance, semivariance is typically more sensitive to proportional effects (correlation between local means and variances) than other measures of spatial correlation or covariance; proportional effects cause “first—order” or linear trends in a variogram and can result in overestimation of the range of distances at which data are correlated (Rossi et al. 1992). Generally, this property of variograms is highly undesirable for spatial modeling. For the purposes of identifying large—scale invasion stage transitions, however, variogram analysis may be superior to other spatial analyses because first—order trends in variograms may be indicative of important mechanistic changes in invasion dynamics. Given a general understanding of invasion dynamics, specific predictions can be made regarding the pattern of semivariance expected at each invasion stage (establishment, expansion, and stabilization). During the establishment stage, the invasion is likely restricted to a small area relative to its potential range. Therefore, only local, small-scale patterns of semivariance should be observed. If sampling occurs over an area equal to or larger than the area occupied, a bell- or wave-shaped variogram may result (Villard and Maurer 1996, Pebesma 1999). In this case, low semivariance (high spatial autocorrelation) is observed not only between pairs of sample points located close together, but also between pairs of sample points located far apart at opposite edges of species’ range. Movement of an invasion wavefront across the study area during the stage of expansion is a process that may result in first—order variogram trends. A first— order trend indicates that a pattern is emerging in the data that is larger than the lag classes modeled; therefore, local means and variances will differ in different locations and directions (Rossi et al. 1992). During the expansion stage Of invasion, a gradient of high to low population index values may occur that would appear in a variogram as a straight line (Rossi et al. 1992, Bailey and Gatrell 1995, Legendre and Legendre 1998). During stabilization, populations begin to settle and resume normal fluctuations with no new large additions to the range unless a species is migratory or overcomes a large physical barrier to expansion (Isard and Gage 2001). Therefore, at this stage one might expect to observe an untrended variogram pattern in which semivariance increases until it reaches the maximum distance at which data are spatially dependent and then levels off. Obtaining the data necessary to conduct variogram analysis and test these predictions requires careful, long—term monitoring over a large area. One such monitoring program tracked the invasion of the gypsy moth (Lymantria dispar) across Michigan (Gage et al. 1990, Yang et al. 1998). The gypsy moth is an exotic insect native to Europe that has caused extensive defoliation across much of the eastern United States since its introduction to Massachusetts in 1868 (Elkinton and Liebhold 1990). The first breeding population of gypsy moths in Michigan was the result of an independent introduction in 1954 and a second possible reintroduction occurred near Midland in the 19803 (O'Dell 1955, Hanna 1981). Gypsy moths quickly spread across the state, reaching most areas of the Lower Peninsula by the early 19903 and the western Upper Peninsula by the late 19903 (Lele et al. 1998, Yang et al. 1998). In this chapter, my goal was to make and test predictions about the types of models that would best fit expected patterns of variogram behavior at each stage of the invasion process. My objectives were to (1) model and characterize variogram behavior over the course of an invasion, (2) compare variogram behavior at different spatial 10 extents over the course of an invasion, and (3) present the potential uses of variogram analysis as a tool in monitoring the progress of additional biological invasions. Methods Pheromone—baited traps were placed in sections eight and 26 of every township in Michigan over 12 years from 1985 to 1996 (Gage et al. 1990, Yang et al. 1998). At each trap, the total annual male moth catch and the trap location in UTM coordinates were recorded. Traps were monitored between 1985 and 1996 in Michigan’s Lower Peninsula (LP) and between 1986 and 1996 in Michigan’s Upper Peninsula. A subset of 1,090 traps was selected for this analysis from over 3,000 placed in a regular grid across the state. Only traps Operated all 12 years at the same location were selected to avoid change of support in the analysis (Isaaks and Srivastava 1989). Traps located in the Upper Peninsula were excluded because the invasion was still in the early stages by 1996. Variogram analysis was conducted at two different scales for this study. First, semivariance was calculated for traps located across the entire LP. Semivariance was calculated as: NU?) 2N(h),= _—_le(xi_ xi+h) 7(h)— - where y(h) is the semivariance for the distance interval h, N(h) is the total number of sample pairs for the distance interval h, Xi is the measured sample value at location i, and Xi +h is the measured sample value at point i + h (Isaaks and Srivastava 1989). Second, data were grouped and analyzed by ecoregion using subsection delineations (Albert et al. 1986) that were modified to ensure continuous yet distinct groupings (Figure 1.1). 11 Ecoregions were chosen as units Of analysis because they represent areas of similar habitat and climate without regard for political boundaries. Isotropic variograms were calculated for each unit of analysis (entire LP or individual ecoregion) in each year of the study using the geoR package in R (Ribiero and Diggle 2001, R Development Core Team 2004). Semivariance was calculated for 13 distance intervals or lag classes; the number of pairs in each lag class ranged from 9,643 to 59,594 pairs for the LP analysis and 16 to 1,983 pairs for the ecoregion analysis. For the LP analysis, the maximum distance was set at 300,000 111, roughly representing the north—south length of the LP; for ecoregion analyses, the maximum distance was set at 100,000 m, roughly representing the length of most ecoregions. Large maximum distances were used to observe variogram behavior across the extent of the unit being analyzed. To characterize the shape of each variogram, seven models were fit to each variogram using maximum likelihood estimation (Bumham and Anderson 2002). Although maximum likelihood may produce biased variogram parameter estimates in cases of small to moderate sample size (Cressie 1993), this study included a sufficient number of data points (>1,000) to avoid such problems. Each model included 3—4 of the following parameters: range (a), the maximum distance at which data are spatially dependent; sill (C1), the semivariance value at the which the range is reached; nugget (Co), the component of the variance caused by local variability at scales smaller than the sampling interval; or smoothness (K), a shape parameter signifying the order of the modified Bessel functions of the third kind. 12 As noted above, a bell— or wave—shaped variogram would be expected during the initial stage of establishment. Such variogram patterns were predicted to be best fit with either the cubic model 2 3 5 7 )7(h) = C0 + C1[7(£) —8.75[£) +3.5[flj — 075(2) ] if h < a; a a a a otherwise y(h) = C O or wave (also known as cardinal sine or hole effect) model m) = C0 + C,[1- {5’— sin[11-D] h a ' Linear, first—order trends are expected during the expansion stage that would be best fit with the power model (exponent equal to l) 79(11)==C0 +C1h",00 or possibly the exponential model a }9(h) = C0 +C1I:1——exp [—3—]:H Finally, a variogram was predicted to exhibit untrended behavior once the stage of stabilization is reached. In this situation, semivariance would increase until an asymptote is reached. Such a pattern would be best fit with either the power (exponent < 1), exponential, spherical 3 )9(h) =CO +Cl[1.5fi—0.5[£] Jitha; 7(h) =CO ifh>a a a Gaussian 13 . h2 y(h) = C0 +C,|:l—exp [—3a—z—H or generalized Cauchy model 2 }9(h)=CO+C1 1-— 1+6?) whereK>0,a>0. These seven models were used to characterize variogram behavior because each characterizes a slightly different pattern of spatial autocorrelation. The first five models are bounded models for which the semivariance increases until the range is reached and the variogram levels off (Goovaerts 1997, Chiles and Delfiner 1999). The spherical and cubic models reach their sill at the range. The exponential, Gaussian, and generalized Cauchy models approach their sill asymptotically. The spherical and exponential models exhibit linear behavior near the origin, whereas the cubic, Gaussian, Cauchy, and wave models exhibit parabolic behavior near the origin. The exact behavior of the generalized Cauchy is further governed by the smoothness parameter (K). The power model is an unbounded model that is linear when a = 1, curves upward when a > 1, and curves downward when a < 1; the behavior of the power model at the origin also depends on the value of the exponent. Note that the sill and range in the power model equation are not equivalent to their counterparts in other models because this equation models an unlimited dispersion process and is not intrinsically stationary (Chiles and Delfiner 1999). The best—fit model for each year and unit of study was chosen by selecting the model with the lowest Akaike Information Criterion (AIC) value (Bumham and Anderson 2002). Comparing differences in magnitude of AIC values among models 14 allowed us to evaluate the amount of evidence in favor of each model. Finally, the timing of any observed changes in variogram shape were compared with temporal patterns of total trap catch in each unit of study. Results Variogram behavior in Michigan’s Lower Peninsula Annual patterns of spatial autocorrelation, quantified as variograms of gypsy moth trap catch, changed markedly between 1985 and 1996 in the LP (Figure 1.2, Table 1.1). In the first four years of the gypsy moth invasion (1985—1988), variograms displayed a bell—shaped curve and the best—fit model was cubic, as predicted. In 1989, the bell—shaped variogram was best fit with the wave model, but the cubic model exhibited a very similar AIC value. From 1990 to 1992, variograms were linear in shape and the best—fit model was either the power (1991—1992; a = 1), as predicted, or the spherical (1990) model. By 1993, first order trends were no longer prevalent in the data and semivariance increased steadily to an asymptote. The best—fit models for these untrended variograms were the power (1993—1995; exponent < 1) or the exponential (1996), as predicted. Changes in variogram behavior corresponded only partially with temporal changes in trap catch (Figure 1.3). Variogram patterns changed from bell— shaped to linear in 1990 when the temporal patterns in total catch began to level Off. However, the switch from a linear to an untrended variogram in 1993 did not correspond with a change in temporal trends. 15 Variogram behavior at the ecoregion scale At the ecoregion scale, variogram behavior was more variable and predictions for variogram patterns at each stage of the invasion were only supported weakly by the three largest ecoregions (Table 1.1). Variograms for the High Plains ecoregion were linear for the first five years of the study. The best fit for these variograms was provided by the Gaussian (1985—1986) or wave (1987—1989) model due to a small downward curve in the first few lag classes. By 1990, a small wave pattern emerged that was, not surprisingly, best fit by the wave model. Between 1991 and 1996, an untrended variogram pattern was Observed and the best-fit model was either the Gaussian (1991 and 1994), the cubic (1992), the wave (1993 and 1996), or the power (1995, exponent < 1). The observed switch in behavior of the variograms did not correspond with patterns observed in the time series. Trap catch in the High Plains ecoregion grew exponentially until 1990 when catch began to decline then level off by 1993 (Figure 1.4). Between 1985 and 1993, variograms in the Ionia ecoregion displayed a slight bell—shaped pattern that was best fit by the wave (1985, 1990—1991) cubic (1988—1989) or Gaussian (1986—1987) model. By 1992, an untrended variogram emerged and the Cauchy (1992), wave (1993), spherical (1994—1995) or the power (1996, exponent < l) was the best—fit model. Trap catch in the Ionia ecoregion did not correspond with changes in variogram behavior. Trap catch increased over the course of the study with a slight dip reported in 1992 (Figure 1.4). Variograms for the Washtenaw ecoregion were largely linear in nature until the last year of the study. Four of the first five years (1985, 1986, 1988, and 1989) were best fit with the wave model. In 1987, and 1990—1992, the Gaussian model provided the best 16 fit. Between 1993 and 1995, the power model (exponent ~1) best fit the variogram. By 1996, however, semivariance increased then leveled off; in 1996, the best—fit model was the generalized Cauchy. Inspection of the time series for the Washtenaw ecoregion indicated the population was growing in an exponential fashion over the 12 years of the study (Figure 1.4). Variogram behavior over the study period for the Saginaw Bay ecoregion was highly variable. In the first two years of the study, a bell—shaped curve was observed that was best fit with the cubic (1985) and wave (1986) models. In 1987 and 1988, a linear pattern emerged that was best fit with the wave model. By 1989, the bell—shaped curve was observed again and best fit with the cubic model. Between 1990 and 1992, linear variograms were observed and best fit with either the Gaussian (1990), spherical (1991), or wave (1992) model. In 1993 and 1994, the bell shape emerged again and was best fit using the wave (1993) and cubic (1994) models. In 1995, another linear variogram was best fit with the Gaussian model. Finally, in 1996, an variogram emerged that was linear with a negative slope; none of the seven models fit this pattern well with the possible exception of the Cauchy model. Variogram patterns did not correspond with temporal trends in trap catch. In general, trap catch in the Saginaw Bay ecoregion increased until 1989 when it dropped back to levels Observed in 1985 and then leveled off (Figure 1.4). Variogram behavior for the Huron ecoregion oscillated between a bell—shaped curve and a linear pattern. The bell—shaped curve was best—fit with the cubic (1985 and 1989) and the wave (1986—1987 and 1992) models. Linear variograms were best fit with either the Gaussian (1988 and 1996), power (1990 and 1993), or the wave (1991, 1994— 1995) model. Trap catch in the Huron ecoregion did not correspond with changes in 17 variogram behavior. Trap catch in Ecoregion 2 increased exponentially until 1989 when it began to level off (Figure 1.4). Variograms for the Arenac ecoregion were mostly linear with the exception of 3 years in which an untrended variogram pattern was observed. Linear patterns were best fit with either the wave (1985, 1987, 1989-1990, and 1994), power (1988 and 1991; exponent ~ 1), or Gaussian (1993 and 1996) model. In 1986, 1992, and 1995, untrended variogram patterns were observed and best fit with the wave (1985), Gaussian (1992), and cubic (1995) models. Trap catch in the Arenac ecoregion did not correspond with changes in variogram behavior. Trap catch increased exponentially until 1989 when it dropped and leveled off by 1992 (Figure 1.4). Variograms for the Presque Isle ecoregion exhibited a bell—shaped curve pattern until 1996 when an untrended variogram pattern was observed. Bell—shaped variograms were best fit with either the wave (1985, 1988—1994) or Gaussian (1986—1987 and 1995) model. The Gaussian was also the best—fit model in 1996 when semivariance increased and leveled off at the sixth lag class. Changes in trap catch did not correspond with changes in variogram behavior. Trap catch in the Presque Isle ecoregion exhibited a possible time lag between 1985 and 1988 then increased exponentially until 1992 when catch dropped and leveled off (Figure 1.4). Variogram patterns for the Kalamazoo ecoregion were dominated with a bell-— shaped curve until 1996 when an untrended variogram pattern was observed. The bell— shaped pattern was best fit with either the wave (1985 and 1994), cubic (1986, 1990— 1991, and 1993—1994), generalized Cauchy (1987), or Gaussian (1988—1989 and 1992). In 1996, the untrended variogram shape was best fit with the spherical model. The 18 switch between bell—shaped and untrended variogram patterns was not reflected in the graph of trap catch that increased exponentially throughout the course of the study (Figure 1.4). Variogram behavior for the North Allegan/Manistee ecoregion was mostly characterized by a bell—shaped curve from 1985—1992; one exception was the largely flat variogram in 1986. Variograms were best fit with either the cubic (1985, 1988, and 1991—1992), Gaussian (1987), or wave (1989—1990) model. Between 1993 and 1995, an untrended variogram pattern emerged that was best fit with the power (1993,exponent <1), Gaussian (1994), and cubic (1995) model. In 1996, a linear pattern was observed and best fit with the wave model. Again, changes in total trap catch did not correspond temporally with changes in variogram behavior. Trap catch in the North Allegan/Manistee ecoregion exhibited a possible time lag from 1985 to 1988 then increased exponentially until 1990 when it began to level Off (Figure 1.4). The South Allegan ecoregion exhibited largely erratic variogram behavior. The 1985 variogram exhibited a bell-shaped curve that was best fit with the cubic model. Over the next 11 years, however, variograms oscillated between untrended and linear behavior. Untrended variograms were best fit with the exponential (1986), the wave (1987, 1990, and 1995-1996), and spherical (1989 and 1991). More linear variograms were best fit with the exponential (1992 and 1994), wave (1988), and power (1993, exponent ~1). Time series data did not reflect the instability observed in variogram analysis. The South Allegan ecoregion exhibited a possible time lag of nine years followed by exponential growth and a leveling Off after the peak catch in 1993 (Figure 1.4). 19 Discussion Variogram behavior in Michigan ’3 Lower Peninsula The three, distinct variogram patterns displayed by gypsy moth trap catch data across the LP from 1985—1996 likely represented expected patterns of spatial autocorrelation for the invasion stages of establishment, expansion, and stabilization because sampling occured across a large enough spatial scale. During the early stages of invasion, bell—shaped variograms best fit by the cubic and wave models were observed that may signify the establishment stage of invasion during which the population spreads very little as it increases locally. Next, linear variograms were observed that were best fit typically by the power model (exponent = 1); this linear behavior was likely caused by expansion of the gypsy moth as it disperses across the state, resulting in a gradient of high to low trap catch in the north—south direction (Yang et al. 1998). From a purely spatial perspective, gypsy moths could be characterized as having reached the stage of stabilization by 1993 in the LP of Michigan. From 1993—1996, local patterns of spatial autocorrelation were no longer masked by first order trends and an untrended variogram emerged. Gypsy moths were well established across most of the LP by 1993 and the gradient in trap catch was not as strong (Yang et al. 1998). A similar progression of variogram patterns was reported by Liebhold et al. (1991) in their study of historical egg mass distributions over a large portion of northeastern Massachusetts, southeastern New Hampshire, and southern Maine. Liebhold et al. (1991) reported a slightly bell— or wave—shaped variogram in 1910, largely linear variograms between 1911 and 1914, and an untrended variogram in 1915. Given that gypsy moths were introduced to Massachusetts in the late 18803, it is not 20 unreasonable to assume that the gypsy moth invasion of New England had begun to transition between stages over the next 20 years. The observed variogram shifts over time potentially represent changes in invasion stages, but additional studies are needed on new invasive Species using appropriate spatial sampling protocols to unravel this complex process. Variogram behavior at ecoregion scale While variogram analysis for the entire LP showed three distinct behavioral switches, analyses conducted at the ecoregion extent were not as straightforward to interpret. Overall, predictions of variogram shape and model fit at each stage of the invasion were only supported in part by the three largest ecoregions. A variogram progression similar to that of the LP was displayed by the High Plains ecoregion data. Also, the Ionia ecoregion exhibited linear variograms followed by untrended variograms in the early 19903 in a fashion similar to the LP. These observed variogram patterns make sense given the invasion history of the area; gypsy moths partially invaded this ecoregion before the monitoring period began, then halted or slowly spread across it, establishing a gradient in trap catch until the early 19903 when the invasion began to spread more extensively across south—central Michigan (Yang et al. 1998). The Washtenaw ecoregion was the last ecoregion in the LP to be completely invaded by gypsy moths; this ecoregion contains the city of Detroit and the surrounding suburbs and, in general, is a highly developed urban area with little suitable gypsy moth habitat. Until the last year of the study, variograms were largely linear in nature, indicating that the invasion wavefront was late in arriving in this region of the state. The remaining 7 21 ecoregions exhibited variogram behavior changes over time that were dissimilar to that of the LP. The Saginaw Bay, Huron, and Arenac ecoregions contained the original establishment area for gypsy moths in Michigan (Hanna 1981), so trap catch in these ecoregions reflect numbers Observed after the wavefront has passed. Not surprisingly, variogram patterns in these ecoregions were different fiom patterns observed in ecoregions colonized later in the study period. The alternation between bell—shaped and linear (Saginaw Bay and Huron) or between linear and untrended (Arenac) variograms may reflect changes in spatial autocorrelation that accompany natural cycling patterns of gypsy moths after population stabilization. Other patterns of variogram change over time were more difficult to explain. Both the Presque Isle and Kalamazoo ecoregions exhibited bell—shaped variograms until an untrended variogram pattern emerged in 1996. However, the invasion history of these two ecoregions is strikingly different. The Presque Isle ecoregion was invaded within a two—year period (1990—1991), whereas the invasion wavefront did not hit Kalamazoo until 1993. It is unclear why both these ecoregions with dissimilar histories and shapes would exhibit similar variogram patterns over time. Another mixed series of all three variograms types were displayed in the North Allegan/Manistee ecoregion, perhaps because two somewhat dissimilar, narrow ecoregion subsections were inappropriately combined to create the area of analysis. In addition, the South Allegan ecoregion alternated erratically between linear and untrended variograms over the study period. This ecoregion is also long and narrow and may not be an appropriate shape for isotropic variogram analysis. 22 In summary, only the three largest ecoregions (High Plains, Ionia, and Washtenaw) exhibited variogram patterns over time that were consistent with their invasion history, my initial predictions, and the patterns observed at the LP analysis extent. This Observation is also supported in part by an egg mass study conducted by Liebhold et al. (1995). Of the 4 large spatial datasets analyzed, only the dataset largest in area (the tri—state dataset) exhibited the three—stage variogram patterns observed in the Michigan LP analysis; an additional dataset collected over three years across the state of Massachusetts was too short in duration for fair comparison. Therefore, the use of variogram analysis to monitor invasion progress may need to be conducted across relatively large areas (e.g. area of the High Plains ecoregion is 21,000 kmz). These results reaffirm the fact that variograms are designed to statistically model structural spatial dependence, or large—scale trends using many data locations simultaneously (Rossi et al. 1992). However, the exact size of the monitoring region will likely depend on the dispersal ability and other life history characteristics of the species in question. Timing of shifts in variogram behavior One important clue that variogram analysis may be a vital component in the analysis of invasion progress is that the timing of variogram behavior transitions did not correspond with changes in total trap catch, an index of population growth (Figure 1.3). Without consideration Of spatial patterns, inspection of the time series alone would have indicated that the invasion had stabilized by 1990. A similar discrepancy between variogram patterns and trends in total catch over time was observed at the ecoregion extent. For example, the time series for the Washtenaw ecoregion did not show any 23 stabilization in trap catch, whereas the 1996 variogram hinted at the possibility of a more stable pattern of semivariance late in the study period. A similar disparity between the timing of population growth and range expansion was observed by Maurer er al. (2001) in their study of the European starling (Sturnus vulgaris) invasion of North America. In this case, population size lagged behind geographic range size (kmz) such that marked population grth did not occur until the invasion wavefront began to slow upon arrival at the west coast of North America. Because shifts in variogram behavior over time do not appear to coincide with major changes in observed time series, one or more unidentified processes may be driving the spatial dynamics of an invading species. For example, the final stage of the gypsy moth invasion (stabilization) appeared to coincide with two major changes: 1) a temporal leveling off of total trap catch evident in the time series, and 2) the switch from linear to untrended variograms. While the concept that an invading population will eventually reach its maximum is widely accepted, the concept that local spatial autocorrelation is necessary for invasion stabilization is not as intuitive. One possible explanation may be that stabilization of an invasion requires the establishment of metapopulations (Moller 1996). A metapopulation is a set of local populations within some larger area that experiences inter—population migration (Hanski and Gilpin 1997). Moller (1996) suggested that metapopulation structure may be an important mechanism for buffering the effects of local demographic and environmental fluctuations on the overall invasion. However, barriers to interpopulation dispersal may prevent stabilization of the invasion process (Anderson and May 1986, Moller 1996). Stabilization of invasion dynamics is thought to occur in part by immigration between populations 24 subject to different environmental conditions. Brown and Kodric—Brown (1977) described this phenomenon as the rescue effect because influx of individuals fi'om neighboring populations may “rescue” a population from extinction. Stabilizing metapopulation dynamics may be especially important for small invader populations or populations located in less than optimal conditions (Roughgarden 1986, Stacey and Taper 1992). However, the extent to which metapopulation dynamics or other potential mechanisms contribute to the invasion process is still unknown. Use of variograms in invasion analysis and interpretation If spatial analysis of other invasion data reveals similar spatial patterns that correspond to invasion stages, the use of variogram analysis should prove insightful in the improved monitoring and prediction of invasive spread. However, the use of variograms in this manner will require a change in the way ecologists typically collect, display, and analyze patterns in spatial data. To obtain enough data for proper variogram analysis, extensive long—term, spatially explicit monitoring programs for multiple invaders must be developed that include a static sampling design (i.e. sampling method or location of traps does not change significantly over the study period). Any method of identifying invasion stages will falsely identify stage transitions if there is a change in monitoring effort or sampling methodology (Cousens and Mortimer 1995). Given the ever—increasing number of new invasions and billions of dollars in damage that result (Pimentel et al. 2000), applying more resources and effort into large—scale monitoring programs will likely prove to be a good investment. 25 The second change in approach to invasion analysis involves plotting variograms of data collected across or even beyond the scope of the spatial process (i.e. beyond the invasion wavefront, Rossi et al. 1992, Legendre and Legendre 1998). Even when all zero counts are excluded, a similar bell—shaped curve results that is best fit using either the wave or cubic models. Although these variograms may be informative for exploring invasion stage transitions, bell—shaped variograms should not be used to for interpolative modeling unless localized kriging is employed with an appropriately small search radius (Pebesma 1999). The third change in approach to spatial data will require ecologists to tolerate the presence of first order trends in a variogram. First order trends typically linearize a variogram and are best fit using the power or exponential models. First order trends imply that the attribute value measured at a given point depends on the location of that point in the study area, a clear violation of the assumption of stationarity (no proportional effects) required for many methods of interpolative modeling (Isaaks and Srivastava 1989). Despite the drawbacks of trended data, such patterns are indicative of a species that is undergoing range expansion and should not be ignored. Examining variograms for first order trends may be useful as an important invasion analysis tool. If interpolative modeling or estimating spatial process parameters (e.g. range, nugget effect) are of interest, first order trends should be modeled, if possible, and the residuals used to calculate a variogram and krige (Cressie 1993, Chiles and Delfiner 1999). Variogram analysis has great potential for detecting invasion stages. If an invasion is monitored at an appropriately large spatial extent, changes in variogram behavior over time may provide a quick identification of invasion stage. Invasions 26 contained in the establishment stage are more likely to result in successful eradication or reasonable control; therefore, managers may decide to concentrate efforts on populations identified by semivariance analysis to be in the early stages of invasion. Given that timing of changes in semivariogram behavior did not coincide with major temporal changes in total trap catch, semivariance may be the most reliable method available for identifying invasion stages. Alternatively, a lag between arrival of the invasion wavefront and changes in the number of adult male moths caught may have resulted in the timing discrepancy between spatial and temporal patterns. Further examination of when and why invasions transition between stages may improve prediction of invasion dynamics if important mechanistic changes in population dynamics prove to be involved. 27 Table 1.1. Modeling results from variogram analysis Of annual gypsy moth (Lymantria dispar) trap catch in the Lower Peninsula (LP) and 10 ecoregions Of Michigan, 1985- 1996. The best model (model with lowest Akaike Information Criterion value) of seven models fit to annual variograms of gypsy moth trap catch is listed as “Type”. Note that the exponential (Expon.), Gaussian (Gaus.), and spherical (Spher.) model names were abbreviated. Model shape (Shape) describes the overall shape of the variogram being fit (B = bell or wave shaped, L = linear, and U = untrended variogram with typical increase in semivariance to an asymptote). LP High Plains Ionia Washtenaw Saginaw Bay Huron Year Type Sh. Type Sh. Type Sh. Type Sh. Type Sh. Type Sh. 1 985 Cubic B Wave L Wave B Wave L Cubic B Cubic B 1 986 Cubic B Gaus. L Gaus. B Wave L Wave B Wave B 1 987 Cubic B Wave L Gaus. B Gaus. L Wave L Wave B l 988 Cubic B Wave L Cubic B Wave L Wave L Gaus. L 1 989 Wave B Wave L Cubic B Wave L Cubic B Cubic B 1990 Spher. L Wave B Wave B Gaus. L Gaus. L Powerl L 1991 Power1 L Gaus. U Wave B Gaus. L Spher. L Wave L 1992 Power‘ L Cubic U Cauchy U Gaus. L Wave L Wave B 1993 Power2 U Wave U Wave U Powerl L Wave B Powerl L 1994 Power2 U Gaus. U Spher U Power‘ L Cubic B Wave L 1995 Power2 U Power2 U Spher. U Power‘ L Gaus. L Wave L 1996 Expon. U Wave U Power2 U Cauchy U Cauchy L Gaus. L 1 exponent = 1 2 exponent < l 28 Table 1.1 (cont’d). Arenac Presque Isle Kalamazoo North Allegan South Allegan Year Type Sh. Type Sh. Type Sh. Type Sh. Type Sh. 1985 Wave L Wave B Wave B Cubic B Cubic B 1986 Wave U Gaus. B Cubic B Expon. B Expon. U 1987 Wave L Gaus. B Cauchy B Gaus. B Wave U 1988 Power L Wave B Gaus. B Cubic B Wave L 1989 Wave L Wave B Gaus. B Wave B Spher. U 1990 Wave L Wave B Cubic B Wave B Wave U 1991 Power' L Wave B Cubic B Cubic B Spher. U 1992 Gaus. U Wave B Gaus. B Cubic B Expon. L 1993 Gaus. L Wave B Cubic B Power2 U Powerl L 1994 Wave L Wave B Cubic B Gaus. U Expon. L 1995 Cubic U Gaus. B Wave B Cubic U Wave U 1996 Gaus. L Gaus. U Spher. U Wave L Wave U ' exponent = 1 2 exponent < 1 29 North Allegan! Manistee Allegan Washtenaw . Kalamazoo Figure 1.1. Map of Michigan’s Lower Peninsula showing its location in North America and ecoregion boundaries used for variogram analysis of gypsy moth (Lymantria dispar) catch between 1985 and 1996. 30 200 56000 48000 , 40000 - 32000 , 24000 - 16000 8000 0.7,. ,7 0 100 200 1987 300 1989 300 1986 35000 30000 25000 20000 15000 10000 5000 1988 96000 1990 84000 ‘ 72000 . 60000 48000 36000 24000 . 12000 O . Lag distance (km) Figure 1.2. Annual semivariograms of male gypsy moth (Lymantria dispar) trap catch in the Lower Peninsula of Michigan, 1985—1996. Circles represent the original data and Open squares represent predicted values generated by the best—fit models. Predicted values were generated by either the cubic (1985—1988), wave (1989), spherical (1990), power (1991—1995), or exponential (1996) models. Note the shift in variogram shape from bell—shaped (1985—1989) to linear (1990—1995) to an untrended pattern of increasing semivariance and leveling off at an asymptote. 31 56000 : 49000 - 42000 ' 35000~" 28000 - 21000 . 14000 - 7000 a 0--,,- 0 100 1991 35000* 30000* 25000 ‘ 20000 ~ 15000 10000 - 5000 “ 1993 Senfivafiance 0 100 21000 18000 15000 4 12000 « 9000 A 6000 J 3000 - 1995 0 100 Figure 1.2 (cont’d). 120000 ~ 105000 . 90000 . 75000 ~ 60000 . 45000 - 30000 ~ 15000 - 0--- -- A ~- 200 300 0 100 200 300 1992 40000 35000 - 30000 - 25000 - 20000 ~ 15000 - 10000 . 5000 - O - ,,_, . .. if z. 200 300 0 100 200 300 1994 35000 ' 1996 30000 ~ 25000 — 20000 ~ 15000 . 10000 ~ 5000 ~ 200 300 0 100 200 300 Lag distance (km) 32 300000 250000 200000 ~ 150000- 100000 50000 ~ 0 no, 1985 1987 1989 1991 1993 1995 Ybar Total trap catch Figure 1.3. Total annual male gypsy moth (Lymantria dispar) catch in the Lower Peninsula of Michigan between 1985 and 1996. Arrows indicate years in which variogram patterns changed from bell—shaped to linear (right—facing arrow) and from linear to typical increase to asymptote (left—facing arrow). Note that inspection of the time series alone (without consideration of spatial location of data points) would indicate that the invasion had stabilized by 1990. 33 120000 +HighPlains 100000 « 80000 4 60000 - 40000 , 1985 1987 1989 1991 1993 1995 _= . .3 35000 -- ° sag” Bay § 30000 - ”A Ammc 3 25000 a +Hm°n a; 20000 . Q A 35 15000 ‘ .§ 10000 ~ ‘j ‘ 2A ‘ 4, ‘ 5 5mm— // a - , +5 0 A A _, _ 7, z 7 “,7 a- - ”-2 ,, - [-1 1985 1987 1989 1991 1993 1995 ' —o—North Allegan 40000 I .g—t—SouthAlleJ: 30000 ‘ —x—Kalarmzoo 20000 - 10000 a I O -‘ ___. __" 7 A A A 7., 1985 1987 1989 1991 1993 1995 Year Figure 1.4. Total annual male gypsy moth (Lymantria dispar) trap catch in each often ecoregions of the Lower Peninsula of Michigan between 1985 and 1996. 34 CHAPTER 2 THE RELATIONSHIP BETWEEN INVASION SUCCESS AND LANDSCAPE HETEROGENEITY: A CASE STUDY OF THE GYPSY MOTH IN MICHIGAN Introduction Dispersal success is a measure of the ability of an organism to move away from its site or origin and arrive in a patch of suitable habitat. The success of dispersal events is thought to be affected in part by the structure and composition of the landscape across which organisms move (Doak et al. 1992, King and With 2002). Individual-based simulation models Of dispersal have shown that as much as 89% Of the variation in dispersal success can be accounted for by differences in size and isolation of habitat patches (Gustafson and Gardner 1996). The relationship between dispersal success and landscape structure is of great conservation importance because habitat loss and fragmentation may inhibit the dispersal or migration of organisms across the landscape (Schwartz 1992, Malanson and Cairns 1997, Collingham and Huntley 2000, Mennechez et a1. 2003). Two different methods have been used to quantify the relationship between dispersal success and landscape heterogeneity. The first method tests the linear relationship between dispersal success and landscape structure metrics in order to identify negative effects of habitat fragmentation on dispersal success. Either correlation coefficients are calculated directly, or dispersal success is regressed on several landscape structure metrics individually (Schumaker 1996, Li and Wu 2004). In general, patch- based metrics of landscape structure such as connectivity and patch cohesion are 35 positively correlated with dispersal success (Taylor et al. 1993, Schumaker 1996), whereas metrics that characterize the gap structure of the landscape such as lacunarity (Plotnick et al. 1993, Plotnick et al. 1996, Dong 2000) are negatively correlated with dispersal success. The second method involves identifying potential thresholds (sudden, nonlinear changes) in dispersal success and landscape structure with habitat loss; the levels of habitat loss at which dispersal and landscape structure thresholds occur are then compared (With and Crist 1995, With and King 1999a, Collingham and Huntley 2000). Landscape metrics that exhibit thresholds most similar to thresholds in dispersal success are thought to be the best measurements for predicting the effects of habitat fiagmentation on dispersal success (With 2002). Dispersal success thresholds have been observed in model simulations and empirical data (Wiens and Milne 1989, With and Crist 1995, Schumaker 1996). Landscape structure measurements have been calculated using either real landscapes (classified land cover maps, Gardner et al. 1987) or computer- generated neutral landscapes (With and King 1999a). Neutral landscapes are typically binary (suitable or unsuitable) habitat maps generated using simple random or fractal algorithms; these landscapes are called “neutral” because they do not contain the effects of topography, disturbance history, climate, or ecological processes (Gardner et al. 1987, With and King 1999a, With 2002). Although dispersal success and many landscape structure metrics have been shown to display threshold responses to habitat loss, the critical proportion of habitat at which such thresholds occur (pcrit) is unknown (Andren 1994, With and Crist 1995, With and King 1999b). Percolation models, based on the physics of liquid flow (Stauffer 36 1985), have been used to simulate random dispersal of organisms across a random landscape. Percolation models predict that when proportion of habitat (p) is reduced below pcrit = 0.59, a small loss of habitat will result in a disproportionately large decline in dispersal success (With and Crist 1995). In other words, when < 59% of the landscape is composed of suitable habitat, dispersal success drops dramatically because large habitat patches that span most of the landscape and aid dispersal begin to break down into many small habitat patches (O'Neill et al. 1988, Andren 1994, With and Crist 1995). Thus, habitat fragmentation has a relatively small effect on dispersal success when p > 0.59; once pcrit is reached, habitat fragmentation compounds the negative effects of habitat loss on dispersal success (Wiens et al. 1997, With and King 1999b). However, a much lower value ofpcrit = ~0.05-0.1 has been predicted by models of dispersal across more complex neutral landscapes (Gardner et al. 1987, With and King 1999a, With 2002). Many landscape metrics also exhibit threshold responses to habitat loss that are comparable to thresholds in dispersal success (With and King 1999a, With 2002). With and King (1999a) used dispersal models across neutral landscapes to compare thresholds in simulated dispersal success with thresholds in six landscape metrics: landscape connectivity, average distance between patches, size of largest patch, total number of patches, total length of edges, and lacunarity. The only landscape metric that displayed a critical threshold value similar to that of dispersal success was lacunarity (landscape “gappiness”) which increased suddenly below pcrit = ~0.05-0. 1. The remaining patch- based metrics exhibited either a peaked response (total number of patches and total length of edges) or critical threshold values at p > 0.3 (landscape connectivity, average distance 37 between patches, and size of larges patch). Thus, the authors concluded that lacunarity (gap structure) may be more important than patch structure in determining dispersal success (With and King 1999a). To date, the relationship between real landscape structure and the actual dispersal success of an invasive species has yet to be investigated. Although computer simulations and experimental studies of individual dispersal movements are important to the development of dispersal theory, such data are not easy to obtain across a large spatial extent or over long time periods and, thus, are not practical for use in invasion management. In order to study the effects of landscape structure on biological invasions, measurements of invasion progress must be collected systematically across real landscapes throughout the course of an invasion. One such monitoring program tracked the invasion of gypsy moths (Lymantria dispar) across Michigan (Gage et al. 1990, Yang et al. 1998). The gypsy moth is an exotic insect native to Europe that has caused extensive defoliation across much of the eastern United States since its introduction to Massachusetts in 1868 (Elkinton and Liebhold 1990). The first breeding population of gypsy moths in Michigan was the result of an independent introduction in 1954 and a second possible reintroduction occurred near Midland in the 19803 (O'Dell 1955, Hanna 1981). Gypsy moths quickly spread across the state, reaching most areas of the Lower Peninsula by the early 19903 and the western Upper Peninsula by the late 19903 (Lele et al. 1998, Yang et al. 1998). The goal of this chapter was to quantify the relationship between progression of the gypsy moth invasion wavefront and the structure of Michigan’s landscape. My objectives were to (1) test the linear relationship between dispersal success and seven 38 landscape structure metrics, (2) identify critical thresholds in invasion success and landscape structure with habitat loss, (3) compare empirical critical threshold values to the predictions of percolation and neutral landscape theory (p = 0.59 and p = 0.05-0.1, respectively), (4) explore the effect of scale on my results by repeating all analyses at three spatial scales, and (5) compare analyses that used landscape metrics calculated from two different classified land cover maps of Michigan. Methods Invasion data Pheromone—baited traps were placed in sections eight and 26 of every township in Michigan over 12 years from 1985 to 1996 (Gage et al. 1990, Yang et al. 1998). At each trap, the total annual male moth catch and the trap location in UTM coordinates were recorded. Traps were monitored between 1985 and 1996 in Michigan’s Lower Peninsula and between 1986 and 1996 in Michigan’s Upper Peninsula. A subset of 1,090 traps was selected for this analysis from over 3,000 placed in a regular grid across the state. Only traps Operated all 12 years at the same location were selected to avoid change of support in the analysis (Isaaks and Srivastava 1989). Traps located in the Upper Peninsula were excluded because the invasion was still in the early stages by 1996. Land cover data Analyses of landscape structure were repeated using two 30 m resolution raster images representing the land cover of Michigan. The first map (Map 1) was the Michigan Resource Information System statewide land cover classification (Michigan 39 Department of Natural Resources. 1999). Map l was derived from 1978 color-infrared aerial photographs and depicts 52 categories of urban, agricultural, wooded, wetland, and other land cover types. The second map (Map 2), the 2001 Michigan Gap Analysis Project land cover image created for the Michigan Department of Natural Resources (Donovan et al. 2004), was derived from the classification of Landsat Thematic Mapper 5 and 7 imagery collected during spring, summer, and fall from 1999 to 2001. This map depicts 32 categories of urban, agricultural, wooded, wetland, and other land cover types. Gypsy moths are polyphagous herbivores that prefer oaks (Quercus spp.) and aspens (Populus spp.) but will eat a wide variety of other deciduous tree species as well (Elkinton and Liebhold 1990). Therefore, each land cover map was reclassified so that all deciduous forest cover classes were combined to represent gypsy moth habitat. All remaining types of land cover were considered unsuitable for gypsy moths. This simple reclassification allowed me to explore broad patterns of structure in gypsy moth habitat across the landscape without complicating the analysis with more detailed (and Often less accurately classified, Donovan et al. 2004) distinctions among deciduous tree communities or species (Li and Wu 2004). Scale of analysis I repeated analyses of landscape structure and invasion success at three different spatial extents to assess the effect of scale (Gamder et al. 1989, Doak et al. 1992, Li and Wu 2004). Because female gypsy moths are incapable of flight and male moths generally do not disperse beyond 1 km/year (White et al. 2003), it was assumed that the invasion wavefront was driven largely by wind-dispersed larvae that are passively 40 transported up to ~40 km away from their hatching site (Elkinton and Liebhold 1990). Although a small percentage of larvae may actually survive such a long trip, occasional long-distance dispersal events (in combination with human-aided dispersal of egg masses) likely explain the high observed expansion rates of 20-40 km per year in some areas of the country (Taylor and Reling 1986, Liebhold et al. 1992). Therefore, analysis windows of 75, 45 and 15 km on each side were chosen, representing spatial extents that were slightly larger, slightly smaller, and roughly equal to the spatial scale at which the invasion process was likely occurring. At each scale of analysis, the Lower Peninsula of Michigan was clipped into multiple subsections using square-shaped analysis windows (the shape required for lacunarity calculations). At each scale of analysis, the maximum number of windows that could fit inside the Lower Peninsula of Michigan were created and aligned so as to maximize the number of traps included in the analysis. Altogether, 12 boxes 75 km on each side, 37 boxes 45 km on each side, and 207 boxes 15 km on each side were created and used to clip each land cover map. Characterizing the relationship between invasion success and landscape structure Time series in total annual catch at each trap revealed a dramatic increase (from tens to hundreds or thousands) in the number of moths caught at most traps once about 25 individuals had been caught in a given year. Therefore, year of colonization was defined as the year in which total trap catch reached 25 or greater individuals. To obtain a rough estimate of how close each trap was located to the initial site of gypsy moth introduction in Michigan, the distance from each trap to a trap located near the Midland and Bay 41 County border was calculated. At each scale of analysis, traps were grouped by analysis window and summary statistics were calculated, including average years to colonization, invasion rate (distance from original introduction site/average years to colonization), and variance and range in years to colonization. These summary statistics were considered potential indices of invasion success. Analysis windows were used to clip out subsections of both land cover maps for landscape metric calculations. Each clipped land cover raster file was converted to an ascii text file and imported into the software package APACK Version 2.22 (Mladenoff and DeZonia 2002). Within APACK, proportion of gypsy moth habitat (p) and multiple landscape structure metrics were calculated for each subsection of the landscape including lacunarity, total number of habitat patches, size of largest habitat patch, total length of edges, average area per patch, fractal box dimension, and centroid connectivity. Lacunarity was calculated using moving window sizes of l, 5, 10, 50, 100, 150, 200, and 250 cells. Once all indices and metrics were calculated, invasion success was regressed on landscape structure metrics individually to test for linear relationships between invasion wavefront movement and measurements of habitat fragmentation. Next, all indices and landscape structure metrics were plotted individually against p to identify potential threshold responses to decreasing proportion of habitat in the landscape. To identify the presence of a threshold response, I looked for large declines or increases in the slope of the plot. Subsets of 5, 11, or 31 successive points (at the 75, 45, and 15 km scales of analysis, respectively) were created for each plot. I then Obtained a rough estimate of the derivative at the midpoint of each subset by estimating the slope of a simple linear 42 regression line through the subset of points. For all decreasing monotonic curves, thresholds were defined as the value of p corresponding with the subset midpoint at which the regression slope declined below 10% of the maximum. For all increasing monotonic curves, thresholds were defined as the value Of p corresponding with the subset midpoint at which the regression slope rose above 10% of the maximum. For all peaked curves, thresholds were defined as the value of p corresponding with the subset midpoint at which the sign of the regression slope switched consistently from positive to negative. Inverse predictions were made to calculate standard errors and 95% confidence intervals for each threshold (N eter et al. 1985). Results Relationship between invasion success and habitat loss Proportion of gypsy moth habitat (p) ranged from: 0.13-0.54 and 0.15-0.50 for Maps 1 and 2 in 75 km landscapes; 0.09-0.55 and 0.11-0.58 for Maps 1 and 2 in 45 km landscapes; 0004-071 and 0.02-0.70 for Maps 1 and 2 in 15 km landscapes. No indices of invasion success displayed a threshold response to decreasing proportion of habitat. Among the four indices of invasion success calculated, only range in years to colonization declined significantly (a = 0.05) with increasing p in each landscape (Figure 2.1). Because the pattern of increase in range in years to colonization with decreasing p was similar between land cover maps, only Map 1 results are presented. Relationship between landscape metrics and habitat loss Most landscape metrics exhibited a threshold response to changing p with both land cover maps and at all three scales of analysis (Table 2.1 and Figures 2.2-2.4). 43 Because the overall behavior of landscape metrics in response to increasing p was similar between land cover maps, only Map 1 results are presented. However, several large differences in computed thresholds were observed between the two maps for size of largest patch (at 75 and 45 km), average area per patch (at 75, 45, and 15 km), and centroid connectivity (at 75 and 15 km). Two-fold increases in patch metric thresholds were typically observed for the 75 and 45 km analyses, indicating that Map 2 displayed a sudden increase in habitat fragmentation at a higher value of p than Map 1. In contrast, thresholds occurred at a lower p for Map 2 than Map 1 at the 15 km scale of analysis. Several landscape metrics exhibited threshold responses that were less sudden but similar in location to thresholds predicted by neutral landscape models (pcrit = ~ 0.05- 0.1). Most analyses of size of largest patch and average area per patch displayed critical thresholds at or near p = 0.1 (Table 2.1). Also, at the 15 km scale of analysis, lacunarity exhibited critical thresholds of 0.1 and 0.153 for Maps 1 and 2, respectively; although only thresholds for box size 1 are graphed, similar thresholds were observed for box sizes 5, 10, and 50. Only centroid connectivity calculated using Map 2 at the 75 km scale of analysis exhibited a threshold (pcrit = 0.51) similar to that predicted by percolation theory (pcrit = 0.59); however, a low pcrit = 0.09 was also observed at the 15 km scale of analysis for Map2. All other landscape metrics either displayed thresholds at p > 0.1 or no detectable threshold behavior. Most lacunarity analyses either did not exhibit detectable threshold responses to habitat loss or exhibited a threshold (pcrit = 0.44 at the 45 km scale of analysis for Map 1) larger than that predicted by neutral landscape models. Total number of habitat patches on the landscape peaked as expected for both land cover maps between 44 p = 0.22 and 0.30 for all maps and scales of analysis with the exception of Map 2 at 15 km which exhibited no observable threshold. Total length of edges also peaked at about p = 0.26-0.37. Fractal box dimension began to level off at about p = 0.40-0.43 at the 15 km scale of analysis for both maps. However, this metric did not exhibit a critical threshold for either land cover map at the 75 and 45 km scales of analysis. Relationship between invasion success and habitat fragmentation Although no thresholds were detected in measures of invasion success, range in years to colonization is likely negatively affected by habitat fragmentation because it was linearly related with several measures of landscape structure (Figures 2.5-2.7). Similar results were observed for Map 2, so only Map 1 results are presented. As might be expected, range in years to colonization declined significantly (a = 0.05) with size of largest patch, average area per patch, fractal box dimension, and connectivity at the 75 km scale of analysis for Map 1; also, range in years to colonization increased significantly with increasing lacunarity as expected (Figure 2.5). Total number of patches and total length of edges did not exhibit a significant linear relationship with range in years to colonization as might be expected given their peaked response to declining habitat (Figure 2.2-2.4). Results for Map l at the 45 km scale of analysis (Figure 2.6) were similar to those obtained at 75 km. However, at this scale, range in years to colonization declined with increasing lacunarity; also, range in years to colonization was found to increase slightly with increasing total number of patches. Finally, results for Map 1 at the 15 km scale of analysis (Figure 2.7) were similar to those obtained at the 75 km scale of analysis with the exception that range in years to 45 colonization declined slightly with increasing total length of edges. Also, range in years to colonization was found to increase slightly with increasing total number of patches and decrease slightly with total length of edges. Although many of the 45 and 15 km regression results exhibited significant slopes, most linear relationships were weak (R2 < 0.2) and several of the 15 km analyses were likely influenced by extreme outliers. Discussion Efibcts of habitat loss and fragmentation on invasion success Contrary to the predictions of percolation and neutral landscape theory, the gypsy moth invasion wavefront exhibited a linear response to changes in landscape structure across the Michigan landscape. Percolation and neutral landscape models predict a sudden decline in dispersal success at p = 0.59 (O'Neill et al. 1988) and p = 0.05-0.1, respectively (With and King 1999a). However, range in years to gypsy moth colonization declined linearly with increasing proportion of habitat in the landscape (Figure 2.1), indicating that habitat fragmentation did not compound the negative effects of habitat loss. Although thresholds in invasion success were not detected, the amount of time necessary to complete the invasion of a given area was likely increased by habitat loss and fragmentation. Range in years to colonization was significantly correlated with several landscape metrics at all scales of analysis (Figures 2.1, 2.5-2.7). For example, landscapes with p > 0.4 in Figure 2.1a represent areas of the northern Lower Peninsula where the number of years to colonization were all small even though these traps were located relatively far away from the site of original gypsy moth introduction. Landscapes 46 with p < 0.4 represent areas from across the peninsula that were not invaded as uniformly; in other words, some traps in a given analysis window were invaded early in the monitoring period while others were not colonized for up to 11 years later. Therefore, time to complete colonization of a given area (i.e. uniformity of the invasion) was likely slowed by habitat loss. Also, invasion success decreased with increasing habitat fragmentation. Specifically, range in years to colonization increased as size of largest patch, average area per patch, connectivity, and fractal box dimension decreased; also, range in years to colonization typically increased with increasing lacunarity. The only evidence for a non-linear threshold response of invasion success to changes in habitat fragmentation may be the curvilinear decline of range in years to colonization with increasing landscape connectivity, average area per patch, and possibly size of largest patch at the 15 and 45 km scales of analysis. Overall, most relationships between invasion success and landscape metrics were linear. Thus, for passive dispersers like the gypsy moth, habitat fragmentation may slow the invasion wavefront but not cause a sudden, nonlinear decline below a critical level of habitat loss. One potential reason for the apparent discrepancy between these results and many published dispersal model predictions may be that this study measured large-scale movement of the invasion wavefront instead of relatively small-scale movements of individual dispersers. Most research conducted on the relationship between dispersal and landscape structure has involved either individual-based simulation models (Schwartz 1992, With and Crist 1995, Pitelka 1997, With and King 1999a, Matlack and Monde 2004) or experiments documenting short distance, terrestrial movements of beetles (Wiens and Milne 1989, Wiens et al. 1997). My study differs from all others conducted 47 to date in that it tracks the movement of the invasion wavefront and not the movement of individual dispersers; sudden increases in gypsy moth trap catch data indicate that the invasion wavefront (including larvae and non-vagile adult females) recently arrived into the area near that trap, not that a large number of adult male moths have recently flown into the region. Therefore, the index of invasion success used in this study may be displaying a macroecological phenomenon, an emergent property of the combined individual dispersal movements of all individuals in the population (Brown 1995, Maurer 1999). Emergent properties are common features of large, complex systems, but they are not observable in data collected at smaller scales. Therefore, data collected at the scale of individual dispersal movements may not be sufficient to predict the relationship between overall invasion success and landscape structure. Alternatively, the relationship between invasion success and landscape structure may be linear because it is strongly affected by non-random movement and environmental factors not accounted for in most simulation models. Percolation models assume organisms move randomly across a landscape composed of randomly distributed patches of suitable habitat (Wiens et al. 1997); similarly, neutral landscapes do not represent the true complexity of real landscapes and their associated dispersal models assume relatively simple movements (With and King 1999a). However, gypsy moths are carried passively by directed winds that do not allow for purely random dispersal movement (Elkinton and Liebhold 1990). Increasing habitat fragmentation will, of course, have the effect of increasing the probability that gypsy moth larvae will be deposited by atmospheric motion systems or rainfall events into unsuitable areas and fail to survive (Isard and Gage 2001). Anemochorous species like the gypsy moth (and, 48 perhaps, species that actively fly) may not be as negatively affected by the same loss of habitat and patch connectivity as active terrestrial dispersers. Also, percolation and neutral landscape models do not consider the effects of environmental conditions such as topography, disturbance history, climate, or ecological processes such as competition and predation (Gardner et al. 1987, With and King 1999a, With 2002). Therefore, real invasion data may not exhibit threshold responses because environmental factors may mediate the effects of habitat fragmentation. The exact response of a species to the landscape likely depends on a number of additional factors such as dispersal ability, degree of habitat specialization, and rate of habitat turnover in dynamic landscapes (With and Crist 1995, Isard and Gage 2001, King and With 2002, Matlack and Monde 2004). Effect of scale of analysis and [and cover map Increased variation in range of years to colonization (Figure 2.1) and weak correlations between invasion success and landscape structure (Figures 2.5-2.7) at the two smaller analysis scales (45 and 15 km) suggest that movement of the invasion wavefront is likely occurring at a large spatial extent compared to the majority of individual moth movements. By comparing results from different scales of analysis, a breakdown in ability to characterize invasion success was detected when analysis windows of < 75 km on a side were used. Even though most gypsy moth larvae do not disperse 40 km (Elkinton and Liebhold 1990), a small number may survive long distance dispersal events. Those individuals may then begin forming their own new colonies, or “nascent foci” (Moody and Mack 1988), ahead of the invasion wavefront that eventually merge with the wavefront. The formation of nascent foci typically results in much higher 49 rates of expansion (Yang et al. (1998) reported 6 km/year in Michigan) than would be expected given the average individual dispersal distance (Moody and Mack 1988, Shigesada et al. 1995). With (2004) has suggested that landscape pattern is less crucial for predicting colonization success if even a few long-distance dispersal events are successful. However, these results indicate gypsy moth invasion success responds to landscape structure at large spatial extents; thus, occasional long-distance dispersal events and their interplay with landscape structure may be vital to our understanding of the invasion process. If the effects of landscape structure on invasion progress are to be quantified, invasions may need to be monitored, and landscape metrics calculated, at the spatial extent of long-distance dispersal events. Observed thresholds in landscape metrics were remarkably similar between the two land cover maps of Michigan for several landscape metrics despite the fact that these maps were generated almost 20 years apart using different types of imagery and classification methodology; however, three patch-based metrics (size of largest patch, average area per patch, and centroid connectivity) differed greatly between maps. One reason these metrics are sensitive to choice of land cover maps is that they all measure the exact size of patches in the landscape and the two maps used in this study differed in their characterization of patch size. In general, Map 2 paints a much patchier picture of the Michigan landscape than Map 1. For example, the total number of patches/1 ,000 ranged from 6.34-9.45 for Map 1 and 33.8-102.2 for Map 2 on 75 km landscapes. This is likely due to the fact that Map 2 was generated from satellite imagery using a complex classification methodology and because habitat fragmentation likely increased in Michigan between 1978 and 2001. Threshold detection appeared to be reasonably 50 resilient to differences in the magnitude of other landscape metrics that are less dependent on patch size. However, my choice of land cover map would have had a more obvious effect had the landscape not been reclassified to reflect either habitat or non- habitat. Landscape metrics would have been quite different and more complicated to interpret had even a few of the original 52 land cover classes in Map 1 or 32 classes of Map 2 been retained in the analysis. Management implications Gypsy moths did not exhibit a threshold response to declining proportion of habitat, indicating that habitat fragmentation did not compound the negative effects of habitat loss for this species. Thus, gypsy moths may be more resilient to habitat fragmentation than expected. Although habitat fragmentation likely slowed the invasion, fragmentation did not prevent the invasion from reaching all areas of the Lower Peninsula in a relatively short period of time (<10 years in most cases). Therefore, using habitat fragmentation as a “fire-break” to limit invasive spread as suggested by Sharov and Leibhold (1998) and With (2004) will probably not be a successful long-term strategy for generalists like the gypsy moth that show great dispersal capacity. This study also indicates that lacunarity alone may not be sufficient for predicting the response of invasive species to landscape heterogeneity as has been suggested by With and King (1999a, With 2002). First, invasion success did not exhibit threshold behavior similar to that of lacunarity (p = 0.05-0.1), suggesting that invasion success is not affected by habitat loss in the same fashion as lacunarity. Second, my results show that lacunarity was highly sensitive to scale of analysis, indicating that the effects of 51 habitat gaps on invasion success may be different at different scales. Lacunarity only displayed thresholds at p = 0.05-0.1 at the 15 km scale of analysis; at larger spatial extents, thresholds in lacunarity were either not observed, or were detected at p = 0.44. The close correspondence between dispersal success and lacunarity thresholds observed by With and King (1999a) may only be generated by processes occurring at spatial extents similar to individual movements; this relationship may not “scale up” to the movement of an invasion wavefront (Li and Wu 2004). Another important result regarding lacunarity is that it was not the only metric that exhibited thresholds near the dispersal success threshold of p = 0.1 reported by With and King (With and King 1999a). Several patch-based metrics such as size of largest patch, average area per patch, and connectivity also exhibited thresholds at p = ~0.1. Regression analyses also indicated that indices of patch size and distribution may be just as good at describing the negative effects of habitat fragmentation on invasion success as lacunarity. Lacunarity should be used with caution given that the correlation between lacunarity and range in years to colonization vacillated from positive to negative as analysis window size decreased (Figures 2.5-2.7). In contrast, patch—based metrics such as size of largest patch, average area per patch, and connectivity were negatively correlated with range in years to colonization at all scales of analysis. Given the problems associated with the lacunarity metric (especially at larger scales of analysis), I suggest that lacunarity not be used as the primary or sole predictor of an invasive species’ response to landscape structure. This study also suggests that the scale of data collection and analysis must be carefully matched to the scale of the process of interest if invasion management is to be effective. Although studies of individual dispersal movements are important to the 52 development of dispersal theory, such data are not easy to obtain across a large spatial extent or over long time periods and, thus, are not practical for use in invasion management. Therefore, studies conducted on individual dispersal movements may not be applicable to processes occurring at much larger spatial extents and should be used with caution when planning management actions. 53 Table 2.1. Summary of critical thresholds in landscape metrics with habitat loss at three scales of analysis (i.e. analysis window sizes of 75, 45, and 15 km on a side) using two land cover maps. Dashes indicate no threshold response was detected. 75 km 45 km 15 km Map l Map2 Map 1 Map 2 Map 1 Map 2 . 0.44 0.1 0.15 Lacunamy ' ' i 0.02 i 0.001 1» 0.08 Total number 0.28 0.27 0.23 0.30 0.22 _ ofpatches a: 0.36 i 0.63 i 0.76 i 0.81 : 15.12 Size of largest 0.15 0.27 0.14 0.27 0.12 0.18 patch i 0.10 i 0.08 i 0.11 i" 0.13 i 0.002 :t 0.98 Total length of 0.28 0.27 0.32 0.37 0.315 0.26 edges 3: 1.4 i 2.2 i 1.84 x 1.68 i 0.36 i 5.12 Average area 0.15 0.22 0.14 0.28 0.22 0.09 per patch : 0.11 i 0.11 i 0.09 i 0.13 i 0.004 :t 0.08 Centroid 0.28 0.51 0.42 0.42 0.31 0.09 connectivity i 0.15 i 1.2 i 0.1 i- 0.12 :t 0.005 :r 0.11 Fractal box _ _ _ _ 0.43 0.40 dimension i 0.000] t 0.06 54 a) 75 km landscapes . y=-19.71x+ 12.18 R2 = O.7l,p = 0.001 0 0 - b) 45 km landscapes o o y= -6.42.x+ 5.9 . . . R2=0.19,p=0.004 Range in years to colonization (years) 0. 8- A v I 0 0.1 0.2 0.3 0.4 0. 0.6 c) 15 km landscapes 0 y=-3.13x+2.61 . . . . R2=0.14,p<0.0001 0 .0 0.2 0.4 0.6 0.8 Proportion of habitat Figure 2.1. Range in years to colonization among traps located in the same analysis window declined with increasing proportion of habitat. Results are reported using Map l for each scale of analysis: a) 75 km, b) 45km, and c) 15km. 55 Lacunarity H N b) «P U! ON \J 00 Total number of patches/1000 «BMQQOOOO o o 0.6 0.1 0.2 0.3 0.4 0.5 0.6 O o N .o . 0) O 4:. O u: O O 2 2500 1 $2 ‘55 2000 i 0 CL § . . L5 1500 .2 1;; 1000 1 ° . Q) 9 E 5001 . o g..— 0 O a 0.1 0.2 0.3 0.4 0.5 0.6 800 ~ , 1.4 1 750 l 1.21 ' ‘5 600 - ’ 5 '8‘ its 3 0.6 E 550 ~ g .. 500 4 U 0.4 450 ‘ 0.. 0'2 .~ . o o . 400 . . . . 0.0 e . . . . 0,1 0,2 0,3 0,4 0,5 0.1 0.2 0.3 0.4 0.5 0.6 Proportion of habitat Figure 2.2. Thresholds in landscape metrics with increasing proportion of habitat calculated using land cover Map 1 and a 75 km analysis window. 56 1600 ~ 1400 * o 1200 ‘ 1000 1 600 1 400 1 ‘ 200 1 .. . ° ' Average area per patch 00 o o . 0.1 0.2 0.3 0.4 0.5 0.6 1.95 1.90 l 0‘ 1.85 ‘ 1.80 < 1.75 * .' Fractal box dimension 0 1.70 . , . . a 0.1 0.2 0.3 0.4 0.5 0.6 Proportion of habitat Figure 2.2 (cont’d). 57 12 - g 3.5 . 10 ‘ . E 3.0 . o ' a .Q 8i \ 132.51 0" '. o a . LE.- 0 . o . g 6 1 ‘ O 2.0 ‘ '0 ‘ o .. h o . . (U 4 . B 15 4 o .. O .4 4 s. . E . . fl" 2 '-~.._§1.0‘ a... S o 0 1 . . . . ES 0.5 - . . r 4 . 0.1 0,2 0,3 0,4 0,5 01 0.2 0.3 0.4 0.5 0.6 O O S 1200 ~ g . .5 1000 ~ Q. o 2 800 < . . E 600* E20 4001 . ° 3 . - :5 200 l :' . . . 0&1) 0 -L~—‘le ', . '. . . . . E7: 0.1 0.2 0.3 0.4 0.5 0.6 350 ' 3.01 300 ‘ o 0 . 2.5 . . ° C . O . >.. .0 8 250 * 0 ' ° .0 o E 2'0 o 0 3.1.200 5:. ’ 8"5‘ . . '0 c m o. 8 1.0 i . ' 150 ~ . . . r ‘ ' 0.5 . g .- . 100* ' 1 , , . 00 -',"'"';° . . . . 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.6 Proportion of habitat Figure 2.3. Thresholds in landscape metrics with increasing proportion of habitat calculated using land cover Map 1 and a 45 km analysis window. 58 1800 . 1600 1400 1200 1000 800 < . .. ' . 600 . 400 200 0 Average area per patch 0 O _4_42 I O 1.95 « 1.90 4 4,...- l.85 - 0. 1.80 - '- 1.75 - . ' 1.70 « . " l.65« 1.60 a . . . . . 0.1 0.2 0.3 0.4 0.5 0.6 Fractal box dimension Proportion of habitat Figure 2.3 (cont’d). 59 250 g 0.6 200 1 230.5 J . 3 0.4 ...o 0 g 150 o. '. . '. .' ' a “03‘ *4? "° " o 100 802 .p g .l. .u :3 .8 . ‘3‘ ..:§ ’0'.” ' 50 , 3 50.1 1 ‘0'. o . 0 s": , . v . a . 0.0 02 04 06 00 0.2 04 06 08 8 200 ' 2123‘ = E 140 - .354. ‘3 120 ‘ 3' Q 00. o 1;; 100 ‘ o . . ‘ 15° 28‘ -'§"=-- c6 . o. 9 "" o. . 0.8' 95 40 0. .‘. Q) E CD 400 - Connectivity p—s O O 300 . 200 . 0.0 0.2 0.4 0.6 0.8 . MW”. 0.0 0.2 0.4 0.6 0.8 Proportion of habitat Figure 2.4. Thresholds in landscape metrics with increasing proportion of habitat calculated using land cover Map 1 and a 15 km analysis window. 60 50- 0000 4321 83835. 0.6 0.8 0.4 0.2 0.0 10000 1 8000 . 6000 * 4000 ~ 2000 ~ :83 Ba 83 owmco>< 0.4 0.6 0.8 0.2 0.0 l 0.8.6. 211 865:: 4 1 cxo 1.2 1 n F—d 0. 1 308m 8. O 0.2 0.4 0.6 0.8 Proportion of habitat 0.0 Figure 2.4 (cont’d). 61 12~ 10* Range in years to colonization (years) 12 « y = 0.74x + 0.20 010-50500 ° R2=0.06, =0.49 - 10 . p . 8 « o o 6 / 4 q o o y= l.l9x+ 1.12 2‘ . . R2=0.53,p=0.01 , . . . a a . 0 . . . 3 4 5 6 7 8 6 7 8 9 10 Lacunarity Total number of patches/1 000 12 1 )1); -0.004x+ 8.44 ' = .7, = .008 10 . . 05 p 0 8 . 6 . 4 2 . 0 . , . . . 0 500 1000 1500 2000 250( Size of largest patch / 1000 (m2) 12 y=-0.015x+ 15.34 1 R2 = O.27,p= 0.10 V 400 450 500 550 600 650 700 750 800 Total length of edges/1000 Figure 2.5. Relationship between range in years to colonization among traps located in the same analysis window and landscape metrics calculated using land cover Map 1 and a 75 km analysis window. 62 . y = -0.01x + 9.49 10 1 . R2=0.74,p=0.001 8 1 6 . 4 . 2 . 0 . . . . 0 400 800 1200 160( A Average area per patch (m2) g 12 ' )1); -6.78x + 8.09 .5 10 q . = 0.50, p = 0.02 E; o 'E 8 .9. 8 6~ 8 4 8 >\ 2 ~ 0 .E «in 0 , . . . . . —. a? 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Connectivity 12 ~ 10 1 8 . 6 . 4 4 y = -37.55x + 74.8 2 . R2 = 0.63, p = 0.004 0 fl 1.70 1.75 1.80 1.85 1.90 1.95 Fractal box dimension Figure 2.5 (cont’d). 63 10 - 10 * y=2.l9x—O.15 8. , ' .y=0.266x+2,74 84 R2=0.45,p<0.000'l .- , . R2=O.08,p=0.05 6. 4. 2. O p O p fl 0.0 0.1 0.2 0.3 0.4 0.5 0. r 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Lacunarity Total number of patches/ 1000 10 , , y = -0.003x + 4.87 8 R2 = 0.20, p = 0.004 o . . .r , . . . 0 200 400 600 800 1000120( Range in years to colonization (years) Size of largest patch/1000 (m2) 10 .y=-0.002x+4.38 .. 8 fl R2=0.003,p=0.25 6 . . . . 4 .W 2 . . .. . . . 0 e 100 150 200 250 300 350 Total length of edges/1000 Figure 2.6. Relationship between range in years to colonization among traps located in the same analysis window and landscape metrics calculated using land cover Map l and a 45 km analysis window. 64 y = -0.003x + 5.30 R2 = 0.27, p = 0.0006 0 . .. - . fl 0 400 800 1200 1600 200( Average area per patch (m2) y=-l.26x+5.00 8 . , 122:0.22, p= 0.002 Range in years to colonization (years) O . . = f . . . 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Connectivity 10 1 8 w . y=-9.26x+20.71 . . R2=0.12,p= 0.02 6 « o o o o 4 . 2 . 0 c 1 .6 l .7 1.8 1.9 2.0 Fractal box dimension Figure 2.6 (cont’d). 65 _ qy=6.85x+0.14 . . R2 = 0.18,p < 0.0001 6 1' 6 . o 4 4 + 2 y = 0.02): + 1.55 R2 = 0.03, p = 0.004 2 l 0 . ' ' . ' 0 ~ 0 50 100 150 200 250 .0 0.1 0.2 0.3 0.4 0.5 0.6 Lacunarity Total number of patches/1000 8 y=-0.le+2.l - R2=0.12,p< 0.0001 0 40 80 120 160 200 Size of largest patch/1000 (m2) Range in years to colonization (years) 8 - y=-0.03x+2.37 o o o 6 ‘ R2=0.02,p= 0.02 o o no 4 ~ 0 a a... o. o 0 10 20 30 40 50 Total length of edges/1000 Figure 2.7. Relationship between range in years to colonization among traps located in the same analysis window and landscape metrics calculated using land cover Map l and a 15 km analysis window. 66 y = -0.0006x + 1.97 R2 = 0.086, p < 0.0001 0 2000 4000 6000 8000 1000' Average area per patch (m2) y = -0.008x + 1.74 R2 = 0.06, p = 0.03 ‘- " fl 0 100 200 300 400 Range in years to colonization (years) Connectivity 8 - y = -3.27x + 7.41 6 R2 = 0.09, p < 0.0001 0.9 1 .2 1.5 1 .8 2.1 Fractal box dimension Figure 2.7 (cont’d). 67 CHAPTER 3 PATTERNS OF SEMIVARIANCE GENERATED BY SIMPLE AND STRATIFIED DISPERSAL Introduction Ecological variables may be correlated across space as a result of both population processes and environmental conditions (Legendre 1993, Koenig 1999). However, patterns of spatial autocorrelation are not static and may change over time as populations fluctuate and respond to a changing environment (V illard and Maurer 1996, Koenig 1999). Brodie et al. (1995) first noted changes in spatial autocorrelation patterns as a result of local dispersal in a study of the post—fire re—invasion of a Populus balsamifera clone. Correlograms, plots of the spatial autocorrelation statistic Moran’s I against distance (Legendre and Legendre 1998), revealed three distinct patterns of spatial autocorrelation representing three stages of clonal development: post—fire colonization, consolidation, and directional expansion (Brodie et al. 1995). Likewise, three general patterns of semivariance were observed during the gypsy moth (Lymantria dispar) invasion of Michigan (Chapter 1). Semivariance is a spatial statistic used to quantify the variance of the differences between all possible points located a given distance apart (Isaaks and Srivastava 1989); the standardized semivariance is equal to 1 minus the correlation when the population mean and variance remain constant across the study area (Rossi et al. 1992). In the gypsy moth invasion study, graphs of semivariance against distance called semivariograms (or simply variograms) were used to depict the average similarity between values based on 68 separation distance between sampling points (Chapter 1). Initially, bell— or wave—shaped variograms were obtained as the initial population became established. As the invasion wavefront began to move across the study area, a gradient of high to low gypsy moth trap catch occurred that produced linear variograms. Finally, as the invasion neared completion in the Lower Peninsula of Michigan, untrended (asymptotic) variogram patterns were observed in which semivariance increased until it reached the maximum distance at which data are spatially dependent and then leveled off. A similar progression of semivariance patterns was reported by Liebhold et a1. (1991) in their study of historical gypsy moth egg mass distributions over a large portion of northeastern Massachusetts, southeastern New Hampshire, and southern Maine. The fact that two species like P. balsamifera and the gypsy moth both exhibit substantial changes in spatial autocorrelation patterning is intriguing because each species has quite different dispersal strategies. The P. balsamifera clone studied by Brodie et al. (1995) expanded primarily via root suckering, an example of simple dispersal. Simple dispersal (also called neighborhood diffusion) involves the movement of dispersing individuals to areas adjacent to their point of origin (Hengeveld 1989). In contrast, gypsy moth larvae disperse primarily through the air via stratified dispersal, a combination of both short and rare, long-distance movements to areas far away from the point of origin (Hengeveld 1989). Stratified dispersal differs from simple dispersal in that occasional long-distance (or “jump”) dispersal events sometimes result in the establishment of new populations independent of the original population (Sharov and Leibhold 1998). To date, the way in which dispersal strategy (simple vs. stratified) affects the spatial structure of populations as they undergo range expansion has not been explored. 69 In this study, my goal was to determine if computer models of dispersal generate distinct patterns of semivariance that change over time as a population grows and spreads across the landscape. My objectives were to (l) quantify and compare changes in patterns of semivariance over time generated by models of both simple and stratified dispersal, (2) characterize the effects of different population growth rates and dispersal distributions on patterns of semivariance, and (3) compare model results to empirical patterns of semivariance observed during the gypsy moth invasion of Michigan (Chapter 1). Stratified dispersal models with long-tailed dispersal distributions typically generate the most realistic estimates of range expansion rate (Clark 1998, Higgins and Richardson 1999, Clark et al. 2001). Thus, I predicted such stratified dispersal models would also generate variogram progression patterns most similar to empirical patterns displayed by gypsy moths (bell-shaped, linear, and asymptotic). Methods Model structure Two models were created in R (R Development Core Team 2004) to simulate population dynamics with either simple or stratified dispersal movements. Population change in both models was simulated using a discrete, stochastic Ricker model with density dependent dynamics and non-overlapping generations Ni+l : Ni exp(r —(uNi ))exp(sz) where is N; is the population size at time step i, r is the maximum per capita rate of growth, u is the effect of intraspecific competition on population growth, .5 is a scale factor that represents the size of random fluctuations in r, and z is a standard normal 70 random variable with a mean of zero and a standard deviation of 1. For this study, u and s remained constant at respective values of 0.005 and 0.1. Each model was run twice with a different value for r. The first value (r = 1.7) represented a rate similar to that exhibited by gypsy moths (Sharov and Leibhold 1998); a second, higher value (r = 2.1) was also selected to generate population dynamics with a two-point limit cycle for comparison (Case 2000). Reproduction occurred at discrete time intervals, and each reproduction pulse was followed by a dispersal pulse. Models were initiated with N, at carrying capacity (K = r/u) in a suitable habitat cell at the center of the landscape. In both simple and stratified dispersal models, the number of individuals that would disperse a given distance away from each cell in each time step was determined, and then the direction of dispersal was selected. To simulate density dependent dispersal, dispersal began in each cell when population abundance reached 75% of K, allowing each population to build up substantially before emigration occurred. For the simple dispersal model, proportion of individuals dispersing out of a given cell was selected from a random uniform distribution between 0 and 0.1. Dispersing individuals were then moved one cell away in a random direction. For the stratified dispersal model, the proportion and number of individuals dispersing a given distance was determined by one of two gamma distributions (Figure 3.1). Distribution 1 used a rate parameter of 0.5 and a shape parameter of 1.2; this distribution was similar to half of a high-peaked normal distribution in that most dispersers moved short distances and few to none moved moderate to long distances (Clark 1998). Distribution 2 was a long-tailed gamma distribution with rate parameter (a) of 0.043 and a shape parameter (6') of 0.5; use of this dispersal distribution resulted in most dispersers moving short distances, but relatively more individuals made 71 moderate to long dispersal jumps across the landscape compared to Distribution 1 (Clark 1998). For both simple and stratified dispersal models, movement occurred in a random direction once dispersal distance was determined for each dispersing individual. Reflecting boundaries were employed such that no individuals were allowed to leave the landscape. Both models were run on two different landscapes consisting of 400 cells. The homogeneous landscape consisted of 100% suitable habitat and was chosen to examine the movement of organisms that were uninhibited by habitat loss or fragmentation. The heterogeneous landscape was a selected section (360 kmz) of the northern part of Michigan’s Lower Peninsula consisting of 54% suitable habitat; this landscape was taken from the 2001 GAP/IF MAP Michigan land cover image (Donovan et al. 2004) in which all deciduous hardwood types were reclassified with a value of “1 ” and all other land cover classes were reclassified as “0” or unsuitable habitat (Chapter 1). In summary, 12 scenarios were created by varying growth rates, landscapes, and dispersal distributions (for the stratified dispersal model) for each of the two dispersal models (Table 1). Each model scenario was run for 100 time steps and repeated 50 times. Characterizing changes in semivariance Semivariance was calculated as: 1 N(h) 2 707) - m #2109 _xi+h) where y(h) is the semivariance for the distance interval h, N(h) is the total number of sample pairs for the distance interval h, Xi is the measured sample value at location i, and 72 Xi +h is the measured sample value at point i+ h (Isaaks and Srivastava 1989). Variograms were created using all data from all cells in every other time step (total of 50/simulation) with 13 lag distance bins and a maximum distance equal to the length of one side of the landscape. To characterize the shape of each variogram, I fit four models to the data generated by each scenario (Table 1) using maximum likelihood estimation (Bumham and Anderson 2002). First, a quadratic model defined as flh) = ah2 + bh + c was fit to the data to characterize bell-shaped variograms; h represents the distance lag class, c the intercept, and a and b the slopes for the second- and first-order terms, respectively. To characterize linear variograms, I fit a linear model flh) = bh + c to the data. I also fit a horizontal line to the data in which the starting value for the constant c was set to the mean semivariance 701) = c in order to characterize variograms with no spatial autocorrelation. To characterize asymptotic (untrended) variograms, I fit the data with the exponential variogram model ;9(h) = C0 + Cl [1 —exp [—3 3]] a which included the following parameters: range (a), the maximum distance at which data are spatially dependent; sill (C 1), the semivariance value at the which the range is reached; and nugget (C 0), the component of the variance caused by local variability at 73 scales smaller than the sampling interval. Models were fit in AD Model Builder (Limited 2000) and the best fit model was chosen by selecting the model with the lowest Akaike Information Criterion (AIC) value (Bumham and Anderson 2002). Weight of evidence in favor of each model was assessed by computing model probabilities using Akaike weights (Bumham and Anderson 2002). Results As range expansion progressed, changes in variogram behavior were displayed by all 12 model scenarios. Four general types of changes were exhibited as indicated by changes in the model selected as best fit most ofien in each time step (Figures 3.2-3.13). The first group of model scenarios (realsiml .7, realsim2.l, realstrat1.7Dl, realstrat2.1Dl) displayed a clear switch from quadratic to exponential in the model most often chosen as best fit, implying that variogram shapes progressed from bell-shaped to asymptotic (Figures 3.2-3.5). One commonality among these four scenarios is that they all modeled dispersal across real landscapes composed of 54% suitable habitat (Table 1). Simple dispersal models (Figures 3.2 and 3.3) took longer to progress from bell-shaped to asymptotic than stratified models (Figures 3.4 and 3.5). The linear model was chosen as best-fit more often in model scenarios with r = 2.1 (Figures 3.3 and 3.5) than in models with r = 1.7 (Figures 3.2 and 3.4); otherwise, no major differences were observed between model scenarios with different population growth rates. Both stratified dispersal model scenarios in this group (realstratl .7D1 and realstrat2.]Dl) generated dispersal patterns using dispersal Distribution 1. In general, confidence intervals around probability estimates were narrow relative to the second and third groups of model 74 scenarios. Also, models chosen as best-fit were clearly the best model most of the time (i.e. had the highest probability of being the actual best model) except when variogram behavior progressed from bell-shaped to asymptotic. The second group of model scenarios (realstratl .7D2 and realstrat2.1D2 in Figures 3.6 and 3.7) displayed a clear switch in the best-fit model from quadratic to a mixture of exponential and (less ofien) linear models early in the time series, implying that variogram shapes progressed from bell-shaped to either asymptotic or linear. One commonality among these two scenarios is that they both modeled stratified dispersal with the long-tailed dispersal Distribution 2 (Figure 3.1). Also, both model scenarios were run on a real landscape composed of 54% suitable habitat. No major differences were observed between scenarios with different population growth rates. Although the exponential had a higher probability of being the actual best model than the linear or quadratic (at time steps > ‘~10), confidence intervals around probability estimates broadened once variogram behavior progressed from bell-shaped to asymptotic or linear. Such broad confidence intervals indicate that evidence in favor of the exponential model over the linear or quadratic models was not strong. The third group of model scenarios (siml .7 and sim2.1 in Figures 3.8 and 3.9) was best fit with the quadratic model early in the time series. However, later in the time series, the number of times the quadratic model was chosen as best-fit declined and the linear (and occasionally mean) model became more prevalent. Later in the time series, the mean (and exponential for sim2.1) model was occasionally chosen as best-fit as well but not as often as the quadratic or linear. These two scenarios both modeled simple dispersal on landscapes composed of 100% suitable habitat. Interestingly, the scenario 75 with a higher population growth rate (sim2.l) took longer for convergence between the top three models; also, confidence intervals around probability estimates for scenario sim2.1 were narrower than that of the model with a lower growth rate (siml .7). The fourth group of model scenarios (strat2.lD2, strat2.lDl , stratl .7D2, and stratl .7D1 in Figures 3.10 to 3.13) produced variogram patterns that were best-fit by the quadratic model a majority of the time. Occasionally, the mean or linear models were chosen as well for model scenarios strat1.7D1 and strat2.lD2 (Figures 3.12 and 3.13), respectively. However, the probabilities that either the mean or linear models were actually the best-fit model were quite low. All four of these model scenarios were stratified dispersal models run on landscapes composed of 100% suitable habitat. Contrary to my prediction, both simple and stratified dispersal models exhibited variogram progression patterns similar to that of gypsy moths in Michigan (Chapter 1) when dispersal was modeled across a real section of the Michigan landscape. Among all 12 model scenarios, scenarios in the first group exhibited variogram progression patterns most similar to those exhibited by gypsy moths (bell-shaped to linear to asymptotic). However, the linear model was not typically chosen as an intermediate best-fit model with the possible exception of scenarios realsim2.l (Figure 3.3) and realstrat2.1Dl (Figure 3.5) in which the linear model sometimes fit the data well. Another discrepancy between these results and the gypsy moth study is that changes in variogram behavior generated by dispersal models showed close correspondence with changes in population abundance over time. Progression in variogram behavior typically occurred when population abundance began to reach an asymptote. No apparent relationship was observed between the presence of limit cycles 76 (generated when r = 2.1) and higher variability in the probability of a given model actually being best-fit. Discussion The results of this study indicate that landscape heterogeneity is the major driving force behind semivariance patterns observed during the gypsy moth invasion of Michigan (Chapter 1). In other words, spatially autocorrelated habitat (not dispersal strategy) may be determining the way in which abundance is distributed across the landscape. Only models of dispersal across sections of a real landscape produced variogram progression patterns similar to that of the gypsy moths. A strong effect of landscape heterogeneity on the range expansion process has been observed in other studies of dispersal and biological invasions as well. Environmentally dependent rates of range expansion have been observed repeatedly in both empirical and simulation studies (Bergelson et al. 1993, Hengeveld and Van den Bosch 1997, Collingham and Huntley 2000, Hastings et al. 2005). Also, landscape heterogeneity may affect dispersal success (With and Crist 1995, Schumaker 1996, King and With 2002, Matlack and Monde 2004), although not necessarily to the extent suggested by many simulation studies (Chapter 2). It is unlikely that unique characteristics of the Michigan landscape caused such variogram progression patterns because similar patterns were observed during the gypsy moth invasion of the east coast of North America (Liebhold et al. 1991). Contrary to my predictions, variogram progression patterns similar to those of the gypsy moth do not appear to be generated by differences in dispersal strategy. Variogram progression patterns were remarkably similar between simple and stratified 77 dispersal models in the first group of model scenarios (realsim1.7, realsim2. 1 , realstrat1.7D1, realstrat2.1D1); the major observable difference was that simple dispersal models took longer to progress from quadratic to exponential. Thus, stratified dispersers like the gypsy moth are likely not the only organisms that produce a progression of variogram patterns from bell-shaped to asymptotic as they undergo range expansion. Although variograms were not calculated by Brodie et al. (1995), they did observe a distinct progression of correlogram change over time; this suggests that variogram analysis of their data might reveal that organisms with dispersal strategies quite different from gypsy moths also exhibit distinct variogram progression patterns. An alternative explanation for why I did not observe typical variogram progression patterns in many of the stratified dispersal model scenarios (group 4) may be that the size of my landscape was not large enough to allow for the process of long- distance dispersal to be modeled properly. The second group of long-distance stratified dispersal scenarios (realstrat2.1D2 and realstratl .7D2) exhibited a switch in variogram behavior from bell-shaped to either linear or exponential; if reflecting boundaries and a larger landscape had been used, this group of model scenarios may have exhibited variogram progression patterns most similar to that of gypsy moths. Schneider (2003) observed a linear relationship between mean movement distances of five butterfly species and the size of the study area across which mark-recapture studies were conducted. Therefore, the study area size may have a major impact on the movement distances observed for volant organisms that exhibit long-distance dispersal. I did not observe many occurrences when the linear model was best fit to variograms in the first group of model scenarios (realstrat1.7D1, realsim2. 1 , realsim2. 1 , 78 realstrat2.1Dl) possibly because of differences between the linear model used in this study and the power model used in Chapter 1. In Chapter 1, I did not restrict the exponent of the power model flh) =C0 +Clha,00 to equal 1; therefore, models that were primarily linear with a slight curve downward at the end of the time series were still best-fit by the power model with exponent < 1. In this chapter, I used a strict linear model, so only variograms with a straight line were best fit with the linear model. If the variograms in Chapter 1 were refit with the models used in this chapter (quadratic, linear, mean, and exponential), the gypsy moth variograms that were best fit by the power model with an exponent < 1 (Table 1.1) in years 1993 to 1995 (Figure 1.2) would likely be best fit by the exponential model and only two of the 12 variograms would have been characterized as linear. The close correspondence between timing of variogram progression patterns and changes in population dynamics in my simulation results indicates that a time lag may be affecting the timing of variogram progression in the gypsy moth monitoring data. Because female gypsy moths are incapable of flight and male moths generally do not disperse beyond 1 km/year (White et al. 2003), the gypsy moth invasion wavefront is driven largely by passively transported, wind-dispersed larvae (Elkinton and Liebhold 1990). However, the gypsy moth monitoring program in Michigan sampled the adult male moth population in a given area (Yang et al. 1998). Therefore, estimates of the timing of spatial changes in population structure may need to be adjusted to account for a time lag if larval movements drive those changes. Although inspection of spatial patterns in variogram behavior has potential implications for detecting stages in the invasion 79 process (Chapter 1), the biology of the invading organism should be taken into account when interpreting the timing of changes in both temporal and spatial population dynamics. Alternatively, the dispersal models used in these simulations may have tracked closely changes in spatial structure because they simulated density dependent dispersal. If gypsy moths do not exhibit density dependent dispersal, a time lag between leveling off of total population abundance and the switch from linear to asymptotic variograms would be expected in the gypsy moth variogram analysis (Chapter 1). Another aspect of my prediction that was not supported by these simulation studies was that gypsy moth-like variogram patterns would only be observed when dispersal was modeled using a long-tailed dispersal distribution. In fact, models that utilized a relatively high amount of long-distance dispersal generated less stable spatial structure (i.e. fluctuations between exponential and linear models) than simple dispersal models. Again, the second group of long-distance stratified dispersal scenarios (realstrat2. 1 D2 and realstratl .7D2) may have exhibited variogram progression patterns most similar to that of gypsy moths if reflecting boundaries and a larger landscape had been used. Although long-tailed dispersal distributions were not required to generate variogram progression patterns, they are often required to accurately predict the rate of range expansion for species with stratified dispersal strategies such as trees (Clark 1998, Higgins and Richardson 1999, Clark et al. 2001), insects (Kot et al. 1996, Lewis 1997, Bailey et al. 2003), and birds and mammals (Van den Bosch et al. 1992, Shigesada et al. 1995, Veit and Lewis 1996). Therefore, different aspects of the range expansion process may depend on different aspects of an organism’s dispersal strategy (Kinlan et al. 2005). Range expansion rate (Moody and Mack 1988, Clark 1998) and genetic connectivity 80 (Trakhtenbrot et al. 2005) may rely heavily on rare, long-distance events. However, maintenance of spatial structure within a species range may be driven by patterns of local dispersal (this paper) across a spatially autocorrelated environment (Brown 1984). Clark has shown that the use of mixed dispersal kernels (dispersal distributions split into local and long-distance components) in range expansion models produce the best estimates of range expansion rate (Clark 1998). Such an approach may be needed to tease apart the relative importance of each dispersal component to the overall range expansion process, and to identify the effect of landscape heterogeneity on each component. 81 Table 3.1. Description of twelve model scenarios, including dispersal strategy modeled, landscape type, population grth rate (r), and dispersal distributions used in model simulations. Distribution 1 is a gamma distribution with rate parameter of 0.5 and shape parameter of 1.2; Distribution 2 is a long-tailed gamma distribution with rate parameter of 0.043 and shape parameter of 0.5. Scenario Dispersal Landscape r Dispersal name strategy drstrrbutron siml .7 Simple 100% habitat 1.7 - sim2.1 Simple 100% habitat 2.1 - strat1.7Dl Stratified 100% habitat 1.7 Distribution 1 strat1.7D2 Stratified 100% habitat 1.7 Distribution 2 strat2.lDl Stratified 100% habitat 2.1 Distribution 1 strat2. 1 D2 Stratified 100% habitat 2.1 Distribution 2 realsim1.7 Simple 54% habitat (real) 1.7 - realsim2.] Simple 54% habitat (real) 2.1 - realstrat1.7Dl Stratified 54% habitat (real) 1.7 Distribution 1 realstratl .7D2 Stratified 54% habitat (real) 1.7 Distribution 2 realstrat2.1Dl Stratified 54% habitat (real) 2.1 Distribution 1 realstrat2.1D2 Stratified 54% habitat (real) 2.1 Distribution 2 82 0.4 - Probability density o .o N t» .O y—d 0.0 . . . . 5 10 15 20 Dispersal distance (number of cells) Figure 3.1. Two dispersal distributions used to determine the proportion of dispersers moving a given dispersal distance (number of cells in the landscape grid). The solid line represents a gamma distribution with shape = 1.2, and rate = 0.5. The dashed line represents a long-tailed gamma distribution with shape = 0.043, and rate = 0.5. 83 50 8e+4 g -o 0 2‘3. 40 . £- 6e+4 ‘fi § 30. 8 O 4e+4 E 8 20 - .§ gES 2e+4 5 10 - JD '1: g o '8 Z 0 ~ 0 g 20 40 60 80 100 8 2‘ 1.2 - — 8 +4 (b) e a E 1.0 i (D :2 - 6e+4 B 0.8 ~ .2 33 a 0'6 ‘ ~ 4e+4 A? :1 0.4 4 '8 '8 - 2e+4 a: 0.2 - 0.0 3 I I I I I 0 0 20 4O 60 80 100 Time step Figure 3.2. Change in variogram behavior for model scenario realsim1.7 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 84 - le+5 g "C d) a - 8e+4 E 8 f - 6e+4 O “O O a m ,“g’ 4e+4 «I: O 312' r 2e+4 E a z '2 0 a? 100 g' 1.2 ~ - le+5 g‘ (b) 5‘ .3 1.0 '1 _ 8 1: 8e+4 E ”,3 0.8 __ a - 6e+4 3 1 a 0.6 Lg 0.4 . .D 2 °‘ 0.2 1 - 2e+4 0.0 l 0 0 20 40 60 80 100 Time step Figure 3.3. Change in variogram behavior for model scenario realsim2.] characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 85 50 0009 00990 00900 90 0 000900 000000 90000 00 — 8e+4 S . . . . o 00 e 0 0 C3 1: .2 3 40 - l ,5 (a) 6e+4 3 .o T) 30 ‘ '8 E ~ 4e+4 3’3 20 - .22 C 2 ~ 2e+4 g 10 ~ g '0 z . '8 0 ” . A . . o g— 20 40 6O 80 100 § 1.2 e 8‘ E o. 10 v I‘ " 'yg ' ! "U'. 'U'." _ 6e+4 g . 1 A " . . .' ’. . . . .‘ '. . .‘.~',‘. 0 .'.'. .'. O o . O .'.'.'o.0.o.o . O o . 8 Probability model is best fit 0 20 40 60 80 l 00 Time step Figure 3.4. Change in variogram behavior for model scenario realstratl .7D1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 86 50 - 9009... ..o 09 o .9999 0.0 o o 090 90 — le+5 O . O. O . .0099 O O. s O O O . 8 . g 40 - 8e+4 E (a) ‘67: B 30 .1 ; ’ 6e+4 g l E g 20 ' ' - 4e+4 .._. l (I; O 33 10 0/ - 2e+4 E H Z 0}_ .‘1.AQ.AAA. OAAAIL. .. .l 0 g 20 4O 60 80 100 g 1.2 - 8 8- - 8e+4 E E it. .:.'.'9 .'.'9'.O. .:.|O.9..:.'......’.;59.009.'.OO.:. .:.O.9. g 8 w o . ('0 :2 0.8 ‘ ~ 6e+4 :g l (b) g 0.6 ~ ‘ .1? l " 4e+4 '2'; 0.4 - .o e .. D.- 0.2 7 A“: L ze+4 :12: o o ifi'“ 7777777 . 7‘4“" 375% 77 . 777‘; 1"? $3533? . . 7777777 7 t I l j 1 0 0 20 40 6O 80 100 Time step Figure 3.5. Change in variogram behavior for model scenario realstrat2.1Dl characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 87 50 1 - 8e+4 a g ( ) 13 Q) '5 40 ~ 7 E r 6e+4 E g 30 ' g - 4e+4 V) t... . S - 2e+4 fig 10 l l V z A/" . g 0 , r , . M— 0 E. 20 40 60 80 100 g S 1.2 1 (b) r 86+4 g. E E. 1.0 a E g ~ 6e+4 n .o ”‘2 E 8 E - 4e+4 .é‘ .5 8 g ‘ — 2e+4 0 100 Time step Figure 3.6. Change in variogram behavior for model scenario realstrat1.7D2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 88 50 ~ le+5 8 a g ( ) d.) 4: g 40 * 86+4 E E) 30 * 6e+4 o E 3 .§ 20 * 4e+4 c: o 0) . "g 10 r 28+4 Z 5 0 n‘ . A A . K IA . A 0 ST 20 40 60 80 100 g' 1.2 - g l se+4 : 10 - (b) g l... ,3; 0.8 - " 6e+4 E g 0.6 - __ ~ 4e+4 g .l- a: 555+“ :_ _ __" 1+: _ g 0.4 ~ " l "‘” ‘"“ " " ___ ‘ .0 ppm m ’ 8 RAM 5) 1 ‘7“ 0,2 - _ q_ ___.~,:-=‘~— —: ,hfrfiztvw E5377: ’ 2e+4 10000. F011); 1110001)"! 1 ~-‘ 00 ‘ ‘ minimum: 0 0 20 40 60 80 100 Time step Figure 3.7. Change in variogram behavior for model scenario realstrat2.1D2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 89 5° - l.6e+5 s (a) 8 - l.4e+5 D 40 - g * l.2e+5 ‘63 3 30 ~ l l.0e+5 a . . . ‘ o . . ”‘5 . 'of ~ 8.0e+4 a ”m 1W GE) 20 7 i . 7‘ ‘ r 6.0e+4 E 10 ” 4.0C+4 3 - g - 2.0e+4 é“ Z O oowooooooooooooo- w‘woov‘ooMoWooao‘. 00 9E; . 8' 20 4O 6O 80 100 3 cr 1'2 ' — l.6e+5 E (b) g 0 1.0 . 'rv~- ‘ r 1.4e+5 a t: . 1 7% l ’ - l.2e+5 B 0.8 ‘l '9‘ - 1.0e+5 g 0 6 o . - g r l‘ ‘t " ‘} "bn 8.0e+4 , . r .s‘ ‘ é 0'4 d a: ‘ AZA‘I‘ ‘} a.) ‘ nai‘ "' 6. 064-4 «s A . IE ‘ I I - 4.0e+4 a. 0.2 - I I I ‘1 k' I ~ 2.0e+4 OO . #9ch ‘o' o o o o o o o o o I a .6 o o o 9 onto 9 o o o'o‘oio..‘.l..,..'.','.'!.L.$.1”.” O O 0 20 4O 6O 80 100 Time step Figure 3.8. Change in variogram behavior for model scenario sim1.7 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 9O 5 0 ~ 2e+5 (6 g (a) t 2e+5 \ {‘9’ 4o - V ‘ * 26+5 7; r le+5 o f 30 _‘ " le+5 g L 1 +5 E e g 20 . ' 86+4 5;: . . )- 6e+4 g 10 - ’ 4C+4 . . “U 0 ~ 2e+4 o z ‘ 7 7 7 ‘ ‘ ‘3 ‘ 70.. £759.72 :70 7 o o IO 0 00300990999094“.oooo""‘. . .. O 93". l' O E 20 4O 6O 80 100 g g. 1.2 - (b) _ l.8e+5 g 10 * l.6e+5 g E ' - l.4e+5 i”; .23 0.8 ‘ ~ l.26+5 T) __ ,_ g 0.6 - . l.0e+5 é‘ 8.0e+4 = 0.4 i“ l E l H- e 6.0e+4 o 5: 0.2 ‘ 4.0e+4 ‘ 2.0e+4 O 20 40 60 80 100 Time step Figure 3.9. Change in variogram behavior for model scenario sim2.1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 91 50 1.6e+5 S (U 73 - 1.4e+5 5 40 . W E (a) ~ 1.2e+5 ‘53 g . 1.0e+5 O) '8 E L 8.0e+4 8 .g — 6.0e+4 “5 a; - 4.0e+4 E Z r 2.0e+4 BU .. '8 . 0.0 g 100 g' g. - l.4e+5 é 7:? 1'0 ‘ W v mm I 2’: I I t t I: g cm) a ‘.‘ u n u a t a a I,— l.2€+5 (D ,9 l I I U) -— 0.8 l _ T) i (b) l.Oe+5 ‘3 E 0.6 - ~ 8.0e+4 .Q‘ = II ~ 6.0e+4 g 0.4 - u ,5: 9' . 4.0e+4 N" 0.2 - "LL “ill . - 2.0e+4 0 0 M}; M . A in: .‘.‘!‘_‘J‘s'!a‘s"a‘u‘s‘ '4‘". , .v. .v. .v. in ‘1'". . 0.. QOOQOOOOOOOQOQO 906066000OO0,0.0,0000,0.00,000(60.0},0 00 O 20 40 60 80 100 Time step Figure 3.10. Change in variogram behavior for model scenario stratl .7D2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 92 Number of times model best fit the data aoucpunqe uouelndod Probability model is best fit Time step Figure 3.11. Change in variogram behavior for model scenario strat2.lD1 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 93 50 — l.6e+5 S .8 - l.4e+5 0) 5 40 - 3'3 (a) - 1.2e+5 ‘53 '8 30 - - l.0e+5 E o - . + E 8 0e 4 U) “8’ 20 7 ~ 6.0e+4 <3: 3 ~ 4.0e+4 3 10 « g ~ 2.0e+4 z 4 . a O .v‘A......................-;~Q96‘2‘9'.'OOOOO""‘999909.9....7 0.0 :8... 2° 40 60 80 100 g. :1 1.2 - g— - l.4e+5 g 1.0 ‘ ‘ ‘. 1 _. Z. r. I‘- E‘ E ’1 I ‘ * l.2e+5 8 273 73) 0'8 d l' (b) - l.Oe+5 T) “g": 0.6 + - 8.0e+4 g - 6.0e+4 35 0.4 - (U E r 4.0e+4 9‘ 0.2 - l ‘ i — 2.0e+4 .‘ "loll , 0.0 70.“..o’ooetoooo9099000060'...0’09ooo”"¢o€'¢o.oo0’°7 00 0 2° 40 60 80 100 Time step Figure 3.12. Change in variogram behavior for model scenario strat1.7Dl characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 94 50 ’ r 2e+5 :3 3 ~ ‘ H '. 2e+5 d.) A a . J; t g 40 4 WV ' W “w _ 264.5 3:: *" ~ le+5 § (a) T. 30 ‘ l le+5 B E \ - le+5 (I) g 20 ' I - 8e+4 “E ‘ an “E 10 ‘ . r 4e+4 Z J i - 2e+4 3" “U 0 ‘0'09009609699099000969000606099090000.99906006090 0 g 20 40 6O 80 100 g' 1.2 - - 1.8e+5 “g D. - 1.6e+5 g 1.0 " 8 as ’1 | l u ’ ~1.4e+5 ‘5 5 l t 3; 0'“ I, ll I VI II : : - 1.2e+5 6 ’ ,. ”8 0.64 t 1.0e+5 5. ~ 8.0e+4 E; 0.4 - i - 6.0e+4 '8 d: 0.2 - . - 4.0e+4 ‘ la“ 2.0e+4 0.0 - , . . . . 0.0 o 20 40 6O 80 100 Time step Figure 3.13. Change in variogram behavior for model scenario strat2.lD2 characterized by the total number of times each model was chosen as best-fit to the data in each time step (a) and the mean probability that each model was the actual best model in each time step (b). Mean probability estimates include 95% confidence intervals. Variograms were fit with a quadratic equation (closed circles), a linear model (open triangles), a model of the mean semivariance (closed squares), and an exponential model (open diamonds). The solid line and values on the right ordinate represent total population size averaged over all 50 scenarios. 95 CONCLUSIONS The results of this dissertation illustrate two ways in which landscape heterogeneity affects the range expansion process. First, the spatial distribution of abundance across the landscape appears to be a population response to the composition and configuration of habitat. For simulation models (Chapter 3) to produce semivariance progression patterns similar to those observed empirically (Chapter 1), dispersal had to be modeled across a real landscape that contained a mosaic of habitat and non-habitat patches. The second way in which landscape heterogeneity may affect the range expansion process is by reducing the uniformity and speed with which a species expands across a landscape. In my second chapter, gypsy moth invasion success was shown to decline linearly with habitat loss and fragmentation. Detection of emergent properties Spatial patterns observed in this analysis of gypsy moth range expansion across Michigan reflect the combined individual responses of numerous insects. As shown in Chapter 2, the population response to habitat loss and fragmentation was markedly different from simulation models of individual dispersal events. Therefore, the linear response of the gypsy moth invasion wavefront to landscape heterogeneity may be an emergent property of a complex system. In other words, individual dispersal success may exhibit a threshold (sudden, nonlinear) response to habitat loss and fragmentation that is not exhibited by the entire population when invasion success of the wavefront is measured. 96 In order to detect such emergent properties, large-scale monitoring programs are essential. The spatial sampling protocol used in the Michigan gypsy moth monitoring program was designed to identify large-scale processes such as movement of the invasion wavefront. Most management programs cannot invest resources in tracking individual dispersal events. Therefore, large-scale, long-term monitoring programs like the one used to track movement of the gypsy moth invasion wavefront in Michigan are vital if we are to improve the management of invasive species spread. The best method for spatially sampling a population so that large-scale spatial patterns (e. g. movement of an invasion wavefront or changes in spatial structure of a population over time) can be detected has yet to be determined. The gypsy moth monitoring program in Michigan was ideal for spatial analyses because most traps remained in the same location across the entire 12 years of data collection. Without permanent sampling locations, many spatial statistics cannot be compared among years because of change in support in the data (Isaaks and Srivastava. 1989). Another aspect of the gypsy moth monitoring program that was ideal for spatial analysis is that traps were placed across the entire state. Good spatial coverage across all areas in which an invasive species might spread is necessary if semivariance analysis is to be used to identify invasion stages. Without sampling points located along the edge or outside the area that a species occupies, the bell-shaped variograms that signal the establishment stage (Chapter 1 and 3) may not be detectable. In the future, I plan to explore the effect of different sampling schemes on the detection of changes in spatial structure of expanding populations. 97 Semivariance analysis By comparing methodologies from Chapter 1 with Chapter 3, I learned that non- traditional models of semivariance (e. g. quadratic and mean used from Chapter 3) should be used in place of traditional variogram models (e. g. exponential, spherical, Gaussian from Chapter 1) when semivariograms are being used as a tool in monitoring range expansion. Traditional semivariogram models are not adequate for describing trended data because they do not characterize well the wide range of semivariogram behavior exhibited by trended data (e. g. bell-shaped variograms). Although such trended data should never be used as the basis for interpolative modeling (i.e. kriging) because they do not provide reliable estimates of the range and sill of a variogram, trended semivariograms appear to be useful tools in identifying significant changes in spatial structure of a population undergoing range expansion Future directions This dissertation suggests that landscape heterogeneity is a driving force behind spatial patterns in abundance for populations undergoing range expansion. Therefore, future research in dispersal and invasion ecology should concentrate more closely on the relationship between landscape structure and dispersal. Dennis et al. (1998) and Lele et al. (1998) have found marked improvements in population modeling when spatial variation in environmental disturbance, population growth rates, and dispersal are included. In the future, I will build upon their results by explicitly modeling stochastic, density dependent population grth across an array of different heterogeneous landscapes. I will then compare semivariogram progression patterns generated by 98 different dispersal scenarios across landscapes with substantially different levels of habitat composition and fragmentation (a combination of approaches from my second and third chapters). This approach will allow me to identify what types of structured landscapes generate specific types of semivariogram patterns and also the level of habitat loss at which semivariance progression patterns emerge (Chapter 3). I also plan to test for improvements in range expansion model predictions when realistic population responses to landscape structure are explicitly modeled. I expect that such an approach will greatly improve predictions of the movement of harmful invaders (Moody and Mack 1988, With 2002), the range expansion of native species (Dark et al. 1998, Trakhtenbrot et al. 2005), and the effects of climate change on the maintenance of existing species and ecological communities (Dyer 1995, Higgins and Richardson 1999, Iverson et al. 2004). Finally, I am interested in performing semivariance analysis on species undergoing range contraction to see if semivariance progression patterns exhibited by such populations are reversed in comparison to patterns exhibited by the gypsy moth. 99 LITERATURE CITED Albert, D., S. Denton, and B. Barnes. 1986. Regional landscape ecosystems of Michigan. University of Michigan, Ann Arbor. Anderson, R. M., and R. M. May. 1986. The invasion, persistence and spread of infectious diseases within animal and plant communities. Philosophical Transactions of the Royal Society of London B 314:533-570. Andren, H. 1994. Effects of habitat fragmentation on birds and mammals in landscapes with different proportions of suitable habitat: a review. Oikos 71:355-366. Ashton, P. J ., and D. S. Mitchell. 1989. Aquatic plants: patterns and modes of invasion, attributes of invading species and assessment of control programmes. Pages 525 in J. A. Drake and H. A. Mooney, editors. Biological invasions: a global perspective. John Wiley & Sons, New York, New York, USA. Bailey, R. I., M. E. Lineham, C. D. Thomas, and R. K. Butlin. 2003. Measuring dispersal and detecting departures from a random walk model in a grasshopper hybrid zone. Ecological Entomology 28: 129-13 8. Bailey, T. C., and A. C. Gatrell. 1995. Interactive spatial data analysis. Prentice Hall, New York, New York, USA. Bergelson, J ., J. A. Newman, and E. M. F loresroux. 1993. Rates of weed spread in spatially heterogeneous environments. Ecology 74:999-1011. Brodie, C., G. Houle, and M. Fortin. 1995. Development of a Populus balsamifera clone in subarctic Quebec reconstructed from spatial analyses. Journal of Ecology 83:309-320. Brown, J. H. 1984. On the relationship between abundance and distribution of species. American Naturalist 124:255-279. Brown, J. H. 1995. Macroecology. The University of Chicago Press, Chicago, Illinois, USA. Brown, J. H., and A. Kodric-Brown. 1977. Turnover rates in insular biogeography: effect of immigration on extinction. Ecology 58:445-449. Brown, M. W. 1993. Population dynamics of invading pests: factors governing success. Pages 203-218 in K. C. Kim and B. A. McPheron, editors. Evolution of insect pests: patterns of variation. Wiley, New York. 100 Burgess, T. M., and R. Webster. 1980. Optimal interpolation and isarithmic mapping of soil properties I. The semi-variogram and punctual kriging. Journal of Soil Science 31:315-331. Bumham, K. P., and D. R. Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach, Second edition. Springer-Verlag, New York, New York. Carroll, S. P., and H. Dingle. 1996. The biology of post-invasion events. Biological Conservation 78:207-2 14. Case, T. J. 2000. An illustrated guide to theoretical ecology. Oxford University Press, New York. Chiles, J ., and P. Delfiner. 1999. Geostatistics: modeling spatial uncertainty. John Wiley & Sons, Inc., New York. Clark, J. S. 1998. Why trees migrate so fast: confronting theory with dispersal biology and the paleorecord. American Naturalist 152:204-224. Clark, J. S., M. Lewis, and L. Horvath. 2001. Invasion by extremes: variation in dispersal and reproduction retards population spread. American Naturalist 157:537-554. Cliff, A. D., and J. K. Ord. 1981. Spatial processes: models & applications. Pion, London. Collingham, Y. C., and B. Huntley. 2000. Impacts of habitat fragmentation and patch size upon migration rates. Ecological Applications 10:131-144. Cousens, R., and M. Mortimer. 1995. Dynamics of weed populations. Cambridge University Press, Cambridge. Cressie, N. A. C. 1993. Statistics for spatial data. John Wiley & Sons, Inc., New York. Cronk, Q. C. B., and J. L. Fuller. 1995. Plant invaders: the threat to natural ecosystems. Chapman & Hall, London. Dark, S. J ., R. J. Gutierrez, and G. I. Gould. 1998. The barred owl (Strix varia) invasion in California. Auk 115:50-56. Dennis, B., W. P. Kemp, and M. L. Taper. 1998. Joint density dependence. Ecology 79:426-441. Doak, D. F., P. C. Marino, and P. M. Kareiva. 1992. Spatial scale mediates the influence of habitat fragmentation on dispersal success: implications for conservation. Theoretical Population Biology 41:315-336. 10] Dong, P. 2000. Lacunarity for spatial heterogeneity measurement in GIS. Geographic Information Sciences 6:20-26. Donovan, M. L., G. M. Nesslage, J. J. Skillen, and B. A. Maurer. 2004. The Michigan Gap Analysis Project Final Report. Wildlife Division, Michigan Department of Natural Resources, Lansing. Dyer, J. M. 1995. Assessment of climatic warming using a model of forest species migration. Ecological Modelling 79: 199-219. Elkinton, J. S., and A. M. Liebhold. 1990. Population dynamics of gypsy moth in North America. Annual Review of Entomology 35. French, D. R., and J. Travis. 2001. Density-dependent dispersal in host-parasitoid assemblages. Oikos 95. Gage, S. H., T. M. Wirth, and G. A. Simmons. 1990. Predicting regional gypsy moth (Lymantriidae) population trends in an expanding population using pheromone trap catch and spatial analysis. Environmental Entomology 19:370-3 77. Gardner, R. H., B. T. Milne, M. G. Turner, and R. V. O'Neill. 1987. Neutral models for the analysis of broad-scale landscape pattern. Landscape Ecology 1:19-28. Garnder, R. H., R. V. O'Neill, M. G. Turner, and V. H. Dale. 1989. Quantifying scale- dependent effects of animal movement with simple percolation models. Landscape Ecology 3:217-227. Goovaerts, P. 1997. Geostatistics for natural resource evaluation. Oxford University Press, New York. Griffith, D. A. 1987. Spatial autocorrelation: a primer. Association of American Geographers, Washington, DC. Gustafson, E. J ., and R. H. Gardner. 1996. The effect of landscape heterogeneity on the probability of patch colonization. Ecology 77:94-107. Hanna, M. 1981. Gypsy moth (Lepidoptera: Lymantriidae) survey in Michigan. Great Lakes Entomologist 14:103-108. Hanski, I., and M. E. Gilpin. 1997. Metapopulation biology: ecology, genetics, and evolution. Academic Press, New York, New York, USA. Hastings, A. 19963. Models of spatial spread: a synthesis. Biological Conservation 78: 143-148. 102 Hastings, A. 1996b. Models of spatial spread: is the theory complete? Ecology 77 :1675- 1679. Hastings, A., K. Cuddington, K. F. Davies, C. J. Dugaw, S. Elmendorf, A. Freestone, S. Harrison, M. Holland, J. Lambrinos, U. Malvadkar, B. A. Melbourne, K. Moore, C. Taylor, and D. Thomson. 2005. The spatial spread of invasions: new developments in theory and evidence. Ecology Letters 8:91 -1 01. Hengeveld, R. 1989. Dynamics of biological invasions. Chapman and Hall, New York, New York, USA. Hengeveld, R., and F. Van den Bosch. 1997. Invading into an ecologically non-uniform area. Pages 217-225 in B. Huntley, W. Cramer, A. V. Morgan, H. C. Prentice, and J. R. M. Allen, editors. Past and future rapid environmental changes: the spatial and evolutionary responses of terrestrial biota. Springer-Verlag, Berlin, Germany. Higgins, S. 1., and D. M. Richardson. 1999. Predicting plant migration rates in a changing world: the role of long-distance dispersal. The American Naturalist 153:464-475. Isaaks, E. H., and R. M. Srivastava. 1989. An introduction to applied geostatistics. Oxford University Press, New York, New York. Isard, S. A., and S. H. Gage. 2001. Flow of life in the atmosphere. Michigan State University Press, East Lansing, Michigan, USA. Iverson, L. R., M. W. Schwartz, and A. M. Prasad. 2004. How fast and far might tree species migrate in the eastern United States due to climate change? Global Ecology and Biogeography 13:209-219. King, A. W., and K. A. With. 2002. Dispersal success on spatially structured landscapes: when do spatial pattern and dispersal behavior really matter? Ecological Modelling 147:23-39. Kinlan, B. P., S. D. Gaines, and S. E. Lester. 2005. Propagule dispersal and the scales of marine community process. Diversity and Distributions 11:139-148. Koenig, W. D. 1999. Spatial autocorrelation of ecological phenomena. Trends in Ecology & Evolution 14:22-26. Kot, M., M. A. Lewis, and P. van den Driessche. 1996. Dispersal data and the spread of invading organisms. Ecology 77 :2027-2042. Kowarik, I. 1995. Time lags in biological invasions with regard to the success and failure of alien species. in P. Pysek, K. Prach, M. Rejmanek, and M. Wade, editors. Plant invasions: general aspects and special problems. SPB Academic Publishing, Amsterdam. 103 Legendre, P. 1993. Spatial autocorrelation - trouble or new paradigm. Ecology 74:1659- 1673. Legendre, P., and L. Legendre. 1998. Numerical Ecology, Second edition. Elsevier, New York. Lele, S., M. L. Taper, and S. H. Gage. 1998. Statistical analysis of population dynamics in space and time using estimating functions. Ecology 79: 1489-1502. Lewis, M. A. 1997. Variability, patchiness, and jump dispersal in the spread of an invading population. Pages 46-69 in D. Tilman and P. Karieva, editors. Spatial ecology: the role of space in population dynamics and interspecific interactions. Princeton University Press, Princeton, New Jersey, USA. Li, H., and J. Wu. 2004. Use and misuse of landscape indices. Landscape Ecology 19:389-399. Liebhold, A. M., J. S. Elkinton, G. Zhou, M. E. Hohn, R. E. Rossi, G. H. Boettner, C. W. Boettner, C. Bumham, and M. L. McManus. 1995. Regional correlation of gypsy moth (Lepidoptera: Lymantriidae) defoliation with counts of egg masses, pupae, and male moths. Environmental Entomology 24:196-203. Liebhold, A. M., G. A. Halverson, and G. A. Elmes. 1992. Gypsy moth invasion in North America: a quantitative analysis. Journal of Biogeography 19:513-520. Liebhold, A. M., Z. Xu, M. E. Hohn, J. S. Elkinton, M. Ticehurst, G. L. Benzon, and R. W. Campbell. 1991. Geostatistical analysis of gypsy-moth (Lepidoptera, Lymantriidae) egg mass populations. Environmental Entomology 20:1407-1417. Mack, R. N. 1985. Invading plants: their potential contribution to population biology. in J. White, editor. Studies on plant demography: a festschrift for John. L. Harper. Academic Press, London. Malanson, G. P., and D. M. Cairns. 1997. Effects of dispersal, population delays, and forest fragmentation on tree migration rates. Plant Ecology 131:67-79. Matlack, G. R., and J. Monde. 2004. Consequences of low mobility in spatially and temporally heterogeneous ecosystems. Journal of Animal Ecology 92:1025-1035. Maurer, B. A. 1999. Untangling ecological complexity: the macroscopic perspective. University of Chicago Press, Chicago. Maurer, B. A., E. T. Linder, and D. Gammon. 2001. A geographical perspective on the biotic homogenization process: implications from macroecology of North American birds. in J. Lockwood and M. McKinney, editors. Biotic homogenization. Kluwer Academic Publishers, London. 104 May, R. M. 1974. Biological populations with non-overlapping generations: stable points, stable cycles, and chaos. Science 186:645-647. Mennechez, G., N. Schtickzelle, and M. Baguette. 2003. Metapopulation dynamics of the bog frittilary butterfly: comparison of demographic parameters and dispersal between a continuous and a highly fragmented landscape. Landscape Ecology 18:279-291. Michigan Department of Natural Resources. 1999. Michigan Resource Information System (MIRIS). MDNR 1978 Landuse/Cover. Land Cover / Use 1. Lansing, Michigan: MDNR. Land cover interpreted from aerial photography. Mladenoff, D. J ., and B. DeZonia. 2002. APACK 2.22 User's guide. Moller, H. 1996. Lessons for invasion theory from social insects. Biological Conservation 78: 1 25- 142. Moody, M. E., and R. N. Mack. 1988. Controlling the spread of plant invasions: the importance of nascent foci. Journal of Applied Ecology 25:1009-1021. Neter, J., W. Wasserman, and M. H. Kutner. 1985. Applied linear statistical models. Irwin, Homewood, Illinois. O'Dell, W. V. 1955. The gypsy moth outbreak in Michigan. Journal of Economic Entomology 48: 1 70-1 72. O'Neill, R. V., B. T. Milne, M. G. Turner, and R. H. Garnder. 1988. Resource utilization scales and landscape pattern. Landscape Ecology 2:63-69. Otter Research Limited. 2000. An introduction to AD Model Builder Version 4, Sydney, British Columbia. Pebesma, E. 1999. GSTAT User's Manual. Mathematics Department, Macquarie University, Sydney, Australia. Pimentel, D., L. Larch, R. Zuniga, and D. Morrison. 2000. Environmental and economic costs of non-indigenous species in the United States. BioScience 50:53-65. Pitelka, L. F. 1997. Plant migration and climate change. American Scientist 85:464-473. Plotnick, R. E., R. H. Gardner, W. W. Hargrove, K. Prestegaard, and M. Perlmutter. 1996. Lacunarity analysis: a general technique for the analysis of spatial patterns. Physical Review E 53:5461-5468. 105 Plotnick, R. E., R. H. Gardner, and R. V. O'Neill. 1993. Lacunarity indices as measures of landscape texture. Landscape Ecology 8:201-211. R Development Core Team. 2004. R: A language and environment for statistical computing. in R Foundation for Statistical Computing, Vienna, Austria. Reichard, S. H., and P. S. White. 2003. Invasion biology: an emerging field of study. Annals of the Missouri Botanical Garden 90:64-66. Ribiero, P. J ., Jr., and P. J. Diggle. 2001. geoR: A package for geostatistical analysis. R- NEWS 1:15-18. Richardson, D. M., N. Allsopp, C. M. D'Antonio, S. J. Milton, and M. Rejmanek. 2000a. Plant invasions - the role of mutualisms. Biological Reviews 75:65-93. Richardson, D. M., P. Pysek, M. Rejmanek, M. G. Barbour, F. D. Panetta, and C. J. West. 2000b. Naturalization and invasion of alien plants: concepts and definitions. Diversity and Distributions 6:93-107. Rossi, R. E., D. J. Mulla, A. G. Joumel, and E. H. Franz. 1992. Geostatistical tools for modeling and interpreting ecological spatial dependence. Ecological Monographs 62:277-314. Roughgarden, J. 1986. Predicting invasions and rates of spread. in H. A. Mooney and J. A. Drake, editors. Ecology of biological invasions of North America and Hawaii. Springer-Verlag, New York. Schneider, C. 2003. The influence of spatial scale on quantifying insect dispersal: an analysis of butterfly data. Ecological Entomology 28:252-256. Schumaker, N. H. 1996. Using landscape indices to predict habitat connectivity. Ecology 77:1210-1225. Schwartz, M. W. 1992. Modelling effects of habitat fragmentation on the ability of trees to respond to climatic warning. Biodiversity and Conservation 2:51-60. Sharov, A. A., and A. M. Leibhold. 1998. Model of slowing the spread of gypsy moth (Lepidoptera: Lymantriidae) with a barrier zone. Ecological Applications 8:1170- 1179. Shigesada, N., and K. Kawasaki. 1997. Biological invasions: theory and practice. Oxford University Press, Oxford. Shigesada, N., K. Kawasaki, and Y. Takeda. 1995. Modeling stratified diffusion in biological invasions. American Naturalist 146:229-251. 106 Stacey, P. B., and M. L. Taper. 1992. Environmental variation and the persistence of small populations. Ecological Applications 2:18-29. Stauffer, D. 1985. Introduction to percolation theory. Taylor and Francis, London. Taylor, P. D., L. Fahrig, K. Henein, and G. Merriam. 1993. Connectivity is a vital element of landscape structure. Oikos 68:571-573. Taylor, R. A. J ., and D. Reling. 1986. Density/height profile and long-range dispersal of first-instar gypsy moth (Lepidoptera: Lymantriidae). Environmental Entomology 15:431-435. Townsend, C. R. 1991. Exotic species management and the need for a theory of invasion ecology. New Zealand Journal of Ecology 15:1-3. Trakhtenbrot, A., R. Nathan, G. Perry, and D. M. Richardson. 2005. The importance of long-distance dispersal in biodiversity conservation. Diversity and Distributions 11:173-181. Travis, J ., and C. Dytham. 2002. Dispersal evolution. during invasions. Evolutionary Ecology Research 4:1 1 19—1 129. Van den Bosch, F., R. Hengeveld, and J. A. J. Metz. 1992. Analyzing the velocity of animal range expansion. Journal of Biogeography 19:135-150. Veit, R. R., and M. A. Lewis. 1996. Dispersal, population growth, and the Allee effect: dynamics of the house finch invasion of eastern North America. American Naturalist 148:255-274. Vermeij, G. J. 1996. An agenda for invasion biology. Biological Conservation 78:3-9. Villard, M. A., and B. A. Maurer. 1996. Geostatistics as a tool for examining hypothesized declines in migratory songbirds. Ecology 77:59-68. Wiens, J. A., and B. T. Milne. 1989. Scaling of 'landscapes' in landscape ecology, or, landscape ecology from a beetle's perspective. Landscape Ecology 3:87-96. Wiens, J. A., R. L. Schooley, and R. D. Weeks, Jr. 1997. Patchy landscapes and animal movements: do beetles percolate? Oikos 78:257-264. Williamson, M. 1999. Invasions. Ecography 2225-12. With, K. A. 2002. The landscape ecology of invasive spread. Conservation Biology 16:1192-1203. 107 With, K. A. 2004. Assessing the risk of invasive spread in fragmented landscapes. Risk Analysis 24:803-815. With, K. A., and T. O. Crist. 1995. Critical thresholds in species' response to landscape structure. Ecology 76:2446-2459. With, K. A., and A. W. King. 1999a. Dispersal success on fractal landscapes: a consequence of lacunarity thresholds. Landscape Ecology 14:73-82. With, K. A., and A. W. King. 1999b. Extinction thresholds for species in fractal landscapes. Conservation Biology 13:314-326. Yang, D., B. C. Pijanowski, and S. H. Gage. 1998. Analysis of gypsy moth (Lepidoptera: Lymantriidae) population dynamics in Michigan using geographic information systems. Environmental Entomology 27:842-852. 108 11111111111111:will: