Multigenerational migration of fall armyworm, a pest insect

Multigenerational insect migration commonly expands poleward, but meteorological influences are not clearly understood. We coupled biological and physical processes for the agricultural and invasive pest Spodoptera frugiperda (fall armyworm), by modeling its seasonal migration, and comparing simulated migrations to observed captures, and population genetic markers, at a continental scale. Simulations corroborated the spatial distribution and mixing of Texas and Florida source populations defined by genetic haplotypes. Positive relationships were found between first weeks of simulated and observed immigration, and between genetic and simulated metrics. The capacity to project biotic-, migratory-, and meteorology-induced shifts in insect distributions will aid strategic implementation of crop protection measures and economic analyses of host-resistant germplasm deployment in response to a warming climate.


INTRODUCTION
Migratory flight allows many insect species to rapidly expand their access to resources, which is critically important when local resources are inadequate or depleted (Drake andGatehouse 1995, Chapman andDrake 2019). Migratory pest insects can abruptly invade and exploit new crop production regions before natural enemy populations increase or crop protection activities are performed. However, the timing and intensity of generational insect migrations and the associated influence of prevailing atmospheric factors are poorly understood. Such information will become increasingly relevant as changes in climate alter air transport systems and destabilize ecosystems, with unknown consequences to pest migration patterns although poleward shifts of pest distributions have already been reported (Bebber et al. 2013).
Numerous species of insects engage in longdistance migrations that may require multiple generations (Drake 1985, Stefanescu et al. 2012. The aerial flux of migrating insects is profound, reaching 3200 tons of biomass over the southern UK annually (Hu et al. 2016), and displacements are linked to seasonally beneficial wind directions (Chapman et al. 2008, Wotton et al. 2019. Migrants find refuge, food resources, or reproductive habitat and may benefit by escaping predators and inadequate local resources (Drake and Gatehouse 1995). Further, migrant pests may convey genetic traits that provide tolerance to insecticides (Metcalf 1989) and plant-incorporated toxins (Storer et al. 2010, Tabashnik et al. 2013; may carry pathogens capable of infecting plants (Thresh 1989); and may directly injure plants by feeding (Dhaliwal et al. 2010). Coincidentally, beneficial arthropods (Weyman 1993), birds (Kelly et al. 2013), and bats (McCracken et al. 2008) are also known to migrate at high altitudes within the atmospheric boundary layer (Westbrook 2008) and may track the migration of pest insects to gain favorable prey resources.
Herbivorous migratory insects represent a type of migration that allows migrants to track available food resources (Chapman et al. 2004). Often migrations extend the insect population to geographic locations that are not sustainable for year-round survival, requiring migration back to refuge areas (e.g., winter-breeding regions; Westbrook 2008, Hu et al. 2016, Wotton et al. 2019. A critical difference between migrant and resident insects is in the pulsed nature of migrations and the timing of migrant interactions relative to host productivity (Bauer and Hoye 2014). Pulses of migration allow migrants to quickly exploit ephemeral resources; irrigation and other crop production activities generally help migrants accomplish long-distance seasonal movements for which they might otherwise encounter large areas of unsuitable host resources. Productivity of both crops and non-cultivated plants is subject to climatic variability, which can impact the availability, growth stage, suitability, and geographic distribution of hosts on which migratory insects can feed and reproduce (Arag on and Lobo 2012, Shaw 2016).
Environmental conditions largely regulate the initiation and displacement of insect migrations. Populations of migratory insects may delay emigration flights if the wind speed is excessive. Further, precipitation or cool air temperature may inhibit emigration flights. Vertical stratification of air temperature, wind speed, or wind direction can lead to concentration of migrating insects into relatively thin altitudinal layers (Drake 1985, Drake andReynolds 2012), in which the insects may be immersed in a low-level wind jet of high speed and favorable heading (for arriving at locations with advantageous resources). The development of low-level wind jets commonly occurs during the nocturnal period, when radiative cooling leads to development of a temperature inversion in the atmospheric boundary layer that reduces the frictional effects of surface roughness on the wind velocity at the top of the atmospheric boundary layer (Bonner 1968, Westbrook andEyster 2017).
Simulation trajectory and dispersion models have identified likely patterns of migration, but validation has been rare and typically limited to acute migration events for a single source or generation. Increased availability of atmospheric models, global atmospheric data, monitoring networks, and population genetics techniques (Kim et al. 2010, Nagoshi et al. 2012) now offer the means to both model and validate estimates of insect migration based on serial observations (e.g., trap collections) of insect distributions and genetic markers defining natal sources. New efforts are needed to simulate and validate insect migration models that span multiple sources and generations relative to a dynamic and heterogeneous host crop distribution (Stefanescu et al. 2012, Wang et al. 2019, dynamic meteorological processes, and to assess the relative local immigrant contribution of multiple source populations (Nagoshi et al. 2009).
We selected fall armyworm, Spodoptera frugiperda (J. E. Smith), as the migrant insect species for this study because it is a crop pest of tropical origin that cannot survive extended periods of freezing temperature (Sparks 1979). This means that winter-breeding populations of fall armyworm in the United States are historically found below 28°N latitude in southern portions of Texas and Florida, yet infestations during the subsequent local crop production period routinely extend as far north as Canada, a range approximating 3000 km (Mitchell et al. 1991). These observations reflect an annually recurring multigenerational migration behavior that we previously used to identify long-distance population movements (Westbrook et al. 2016). Expanding from this, multi-year simulations of multigenerational fall armyworm moth migration in North America were generated and tested by genetic and field-monitoring methods with the specific objectives to (1) identify seasonal and annual differences in migration patterns; (2) quantify the accuracy in estimating initial immigration; and (3) quantify the accuracy in estimating the mixing ratio of immigrants from two source populations.

METHODS
We compiled crop and moth development, moth migration simulation, trap capture, and genetic patterns to develop and validate a biophysical process defining multigenerational migration of corn-strain fall armyworm. Because of the lack of haplotypes distinguishing between source locations in the rice strain (RS), this paper is limited to corn-strain populations (the primary strain in corn [Zea mays L.] and sorghum [Sorghum bicolor (L.) Moench]). Migration was modeled in three phases, as previously detailed for one year (Westbrook et al. 2016), which are expanded here for four years and summarized below. The first phase created areas of identified field corn available for oviposition and development. We excluded sweet corn which is frequently treated with insecticides and would not produce significant numbers of fall armyworm moths, relative to field corn, at the scale of this study. A second phase accumulated degree-day values derived from air temperatures which were used to drive phenological development of corn, fall armyworm development and population dynamics, and plant-insect interactions used to define migratory propensity. In the third phase, we executed a migration simulation model by incorporating the distribution and growth stage of corn and the abundance of local and migratory fall armyworm moths within 1600-km 2 grid cells on a daily time step from 10 February to 31 December each year. Each year, simulation results of immigration timing were compared to independently observed trap captures, and simulation results expressing the dynamic geographic pattern of moth migration from two sources were compared to independently measured sources defined by haplotyping specimens from the trap captures.

Corn distribution
We downloaded corn planting distribution for 2011-2014 from the U.S. Department of Agriculture-National Agricultural Statistics Service CropScape-Cropland Data Layer (CDL; http:// nassgeodata.gmu.edu/CropScape/; Han et al. 2012Han et al. , 2014. Within an area of interest spanning 75°W to 105°W longitude, used to define corn production areas for potential infestation by immigrant fall armyworms ( Fig. 1), CDL data were extracted in 40 km latitudinal bands. A total of 66 bands of the 30-m resolution CDL data were compiled, and the planted areas summed and upscaled into 1600-km 2 blocks for spatial compatibility with air temperature data used in simulations of corn growth, insect population development, and migration flights.

Development of corn plants and fall armyworms
We assumed that corn fields were planted on 15 February in the winter-breeding source areas of southern Texas (including Tamaulipas state in northeastern Mexico) and southern Florida and that initial infestation occurred at the 3-4 leaf stage, or 167 DD 10°C after planting, which is corroborated by the timing of larvae in two-to fourleaf corn in northern Tamaulipas, Mexico (Blanco et al. 2014). Corn plant and moth development calculations were driven by air temperature. Archived (40-km resolution) meteorological data from the Eta Data Assimilation System (EDAS) were downloaded from the USDC-NOAA Air Resources Laboratory (http://www.ready.noaa. gov/archives.php; Rolph 2013) to calculate cumulative degree-days that indicate the phenological growth stage of corn plants and fall armyworms. We projected the EDAS data (and CDL data) to latitude and longitude. We extracted (2 m above ground layer) air temperature data at 0600 Universal Coordinated Time (UTC) and 2100 UTC, representing approximate daily minima and maxima, respectively. Daily and cumulative degree-days (DD 10°C and DD 13.8°C ), restricted to a maximum of 30°C for estimating corn plant growth, were calculated to simulate growth of field corn (Neild and Newman 1990) and fall armyworms (Hogg et al. 1982), respectively.

Simulating migration of fall armyworm moths
We simulated a total initial infestation of 1 9 10 6 fall armyworms in corn field source areas each year, comprised of 9 9 10 5 and 1 9 10 5 fall armyworms from southern Texas and southern Florida, respectively, accounting for differences in the respective total area of planted corn at similar infestation densities.
The propensity of migratory flight of S. frugiperda is poorly understood, but based on a related noctuid, corn earworm (Helicoverpa zea), we estimated newly eclosed fall armyworms to be as much as 30 times more likely to migrate from the silking stage than from the whorl stage ❖ www.esajournals.org of corn (Westbrook and Lopez 2010). We assumed a ninefold variation in generational emigration of fall armyworm moths for all locations using an estimated proportion of 0.10 emigrants in moths that developed as larvae in whorl stage corn, and 0.90 emigrants in moths eclosing at the time of silking or later stage corn. The remnant (non-emigrant) proportion of moths (0.90 in whorl stage corn and 0.10 in silking stage corn) was identified as resident moths, which was restricted to local reproduction.
The EDAS meteorological data were input to the NOAA Air Resources Laboratory (ARL) HYbrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (HYSPLIT PC version 4.9; Draxler and Hess 1998). Nightly emigration flights of fall armyworms were simulated from the HYSPLIT concentration model. We defined Pollutant Species in HYSPLIT as TEX (Specie 1) and FLA (Specie 2) for the Texas and Florida source populations, respectively. Although atmospheric deposition values were not defined per se, we assigned dispersed moth values at the end of each nightly flight to respective grid cells for updating moth abundance and oviposition into susceptible corn fields. A point source emission file was defined to run emigration flights at 0000 UTC for a duration equal to the number of hours between sunset and sunrise (scotophase) minus 1 h and based on the nightly availability of migrant moths at each source location (i.e., grid cell). The emigrant portion of fall armyworm populations was translated to the emissions file exactly for each nightly flight. However, for the first emergence of the generation, we distributed the population of newly emerged moths over one week using an exponential decay function to simulate uncertainty in this event where E is the number of emerged moths, and the subscripts i and t indicate day-of-week and total per weekly cohort, respectively. Further, we assigned peak adult moth emergence of the first (F 1 ) generation on dates associated with the end of the week following a cumulative value of 346.2 DD 13.8°C (Hogg et al. 1982) after infestation. Peak moth emergence of the second (F 2 ) and third (F 3 ) generations was estimated to occur on dates associated with cumulative values of 692.4 DD 13.8°C and 1038.6 DD 13.8°C after initial infestation, respectively.
After nightly flights, moths were run through a biological model that linked to corn plant distribution and phenology. Migrants that terminated flight over large bodies of water or areas without corn were removed. We split the remaining population into resident moths and migratory moths for subsequent flights. We estimated the proportion of migrants depending on the suitability of corn plants in the area, which was associated with the crop growth stage based on cumulative DD 10°C within respective grid cells. Resident and migrant moths contributed to infestation contingent on this crop growth stage, but only moths that had not reached their ovipositional capacity contributed to the infestation. We pooled the number of fall armyworm larvae that developed from eggs during a specific week into a single cohort at the end of that week. Cohorts of immatures eclosed as adults (moths) after accumulating 346.2 DD 13.8°C . Subsequently, weekly patterns of simulated abundance of fall armyworm moths were estimated, retaining respective contributions from the TEX and FLA source areas.
Programs that executed the biological model and the HYSPLIT migration model were written in R (R x64 version 3.1.1; R Foundation for Statistical Computing, Vienna, Austria), Surfer (version 12.4.784; Golden Software, Golden, Colorado, USA), and MeteoInfo (version 1.0.9.0; Wang 2014). Collection of one year of data required 70 h of computer processing time, 67 h of which were used to upscale the CDL data from 900 m 2 to 1600 km 2 . The main program that executed the biological model and HYSPLIT model required 6 h to run a one-year simulation. We executed the moth migration model on an Intel Core 2 Quad CPU Q9650 3.00 GHz computer. Initial conditions for the biological model are the same as those noted in Westbrook et al. (2016). Meteorological model parameters have been revised slightly from those included in Westbrook et al. (2016): (1) Nocturnal migration duration was changed from a fixed duration of 12 h to a variable duration (i.e., number of hours between sunrise and sunset minus 1 h) that accommodated seasonal and latitudinal effects; (2) TEX winter-breeding source population remained the same quantity but was distributed over a larger area (i.e., grid cells <28.5°N); and (3) FLA winter-breeding source population remained the same quantity but was distributed over a larger area (i.e., grid cells <28.5°N).

Comparing captured moths and simulated immigration
We compared weeks of first moth capture in pheromone traps with weeks of simulated first moth immigration. We selected pheromone trap data from the PestWatch database for locations with annual total capture ≥10 from UniTraps baited with sex pheromone (Appendix S1: Table S1) serviced once or twice weekly. We summarized, plotted, and analyzed weekly accumulations of simulated fall armyworm migration and total fall armyworm trap captures at each selected trap location using Grapher 11.4.770 (Golden Software, Golden, Colorado, USA), MapViewer 8.0.212 (Golden Software, Golden, Colorado, USA), and JMP 12.1.0 software (SAS Institute, Cary, North Carolina, USA). Due to the large data range and numerous zero values, data were logarithm-transformed (i.e., log 10 (x + 1)). We examined the ability of the model to explain variation in initial immigration using linear regression and a mixed model with year as a random variable and used paired t test to compare the first weeks of simulated immigration to first weeks of capture in traps. Further, natal source (TEX or FLA) attribution at sites with trap capture was estimated using the simulated moth migration data and parameterized as a migration index [((log 10 (FAW FL + 1) À log 10 (FAW TX + 1))/((log 10 (FAW FL + 1) + log 10 (FAW TX + 1))]. The simulated migration index varies from À1 to +1, where negative values indicate higher contributions of the TEX natal origin and positive values indicate higher contributions of the FLA natal origin.

Determination of fall armyworm strain and haplotypes
Specimens were initially identified for strain by analysis of the mitochondrial gene cytochrome oxidase I (COI) COI-891F/COI-1472R ❖ www.esajournals.org PCR amplification product. After polymerase chain reaction (PCR) amplification was completed, 5 units of the restriction enzyme EcoRV (New England Biolabs, Beverly, Massachusetts, USA) was added to each 20-lL PCR mix along with 1 lL of the manufacturer recommended 109 restriction enzyme buffer (final volume taken to 30 lL with water). Restriction digests were incubated at 37°C 1-3 h. For each reaction, 6 lL of 69 gel loading buffer was added and the entire sample run on a 1.8% agarose horizontal gel containing GelRed (per manufacturer's instructions; Biotium, Hayward, California, USA) in 0.59 Tris-borate buffer (45 mmol/L Tris base, 45 mmol/L boric acid, 1 mmol/L EDTA pH 8.0). Fragments were visualized on a long-wave UV light box. Only the RS-associated COI allele has an EcoRV site in the amplified region. Uncut fragments were preliminarily identified as corn strain (CS) and were cut out from the gel. Fragment isolation was performed using Zymo-Spin I columns (Zymo Research, Orange, California, USA) according to manufacturer's instructions. The isolated fragments were sequenced using primer COI-891F by the University of Florida Interdisciplinary Center for Biotechnology Research (UF-ICBR). DNA sequence comparisons were performed using Geneious version 5.6.2 created by Biomatters (http://www.geneious.c om/). The h1-4 haplotypes were identified by the nucleotides found at COI sites 1164 and 1287 (Nagoshi et al. 2009). Each site was associated with two alternative bases to produce four possible haplotypes: h1 (A 1164 A 1287 ), h2 (A 1164 G 1287 ), h3 (G 1164 A 1287 ), and h4 (G 1164 G 1287 ). The frequencies of the four haplotypes were calculated for each collection, as was the quotient of the CS-h4 to CS-h2 frequencies, designated as the h4/h2 ratio.
Haplotype h4/h2 ratio data were calculated from collections of ≥15 captured fall armyworms during a week or pooled over several weeks for numerous locations (Appendix S1: Table S1). Haplotype ratio values defined the estimated source attribution of the local fall armyworm moth population. A haplotype ratio ≤0.5 defined the Texas source population; a value >1.5 defined the Florida source population; and 0.5 < value ≤ 1.5 defined a Mixed source population (Nagoshi et al. 2012). To validate the migration model, the ability of simulated migration index values to categorize source attributions defined by haplotype ratios was tested with logistic regression.

Insect and host plant growth
Corn planting dates (Appendix S1: Fig. S1) based primarily on long-term climatic data followed a northward gradient, but short-term, regional weather further modified geographic patterns of subsequent growth of corn plants within and between years. Climatic effects were reflected in the planting dates and acreage within a given year. Short-term, regional air temperature patterns directly influenced the pattern and growth rate of corn plants and fall armyworms. Further, short-term weather influenced the timing of fall armyworm eclosion (escape of the moth from the puparium), insect-plant interactions driven by physiological development of both crop and insect, propensity for migratory behavior, and ovipositional acceptance behavior, thus creating a dynamic, heterogeneous pattern of emigrant sources and immigrant establishment during the growing season (Fig. 1).

Insect migration flights
Migration patterns were simulated for four years, 2011-2014. Initial flight simulations began with the first projected emergence of adult fall armyworms in the winter-breeding areas, and subsequent migratory projections occurred following the emergence of successive generations of adult progeny in the winter-breeding areas and in areas of successful immigration and establishment. Moths were simulated to migrate at night at an initial flight altitude of 500 m above ground level.
Simulated migration of fall armyworms initially progressed northward from Texas and subsequently veered eastward into the Ohio River Valley as early as 30 June in 2011, 30 July in 2012 and 2014, or 10 September in 2013 (Fig. 2). Fall armyworms were simulated to migrate to the Canadian border of the United States by 10 September in 2011-2014.
The Florida winter-breeding population of fall armyworms was simulated to migrate by 10 May in 2011, 2013, and 2014 but not until 30 May in 2012 (Fig. 2). Simulated migration of the Florida population of fall armyworms advanced northward and along the mid-Atlantic coast of the United States. Aside from 2013, the Florida population of fall armyworms was simulated to migrate westward into the Mississippi River Valley by 30 July in 2011 and 2012 or 10 September in 2014.
The simulated first week of fall armyworm immigration increased from about Week 14 in Texas to greater than Week 34 in the northeastern United States in 2011-2014 (Fig. 3). Simulated first weeks of FAW immigration were regressed with observed weeks of first capture of ≥2 FAW in traps (county level), where ≥10 FAW were captured in a given year and no FAW were collected on the week of first trap inspection (Fig. 4). There was a significant linear regression fit FC = 18.93 + 0.31 (SI), where FC is the week of first captured FAW and SI is the week of first simulated immigration (F = 17.5; df = 1, 104; P < 0.0001). The mixed model also showed a significant effect for trap capture (P < 0.0001) and an interaction with year (P < 0.004).
Simulated migration of moths from source populations in southern Texas and Florida resulted in varying ratios of the respective cohorts (Fig. 5). Simulated fall armyworms in the Plains states from Texas to North Dakota were exclusively from the Texas source population in 2011-2014. Although simulated fall armyworms from the Florida source population were predominant along the Atlantic coast of the UDS in Genetic patterns supported the dynamic geographic patterns that emerged from the migration simulations (Fig. 6). Logistic regression showed that the probability of assignment of natal origin assignment based on haplotype ratio (Nagoshi et al. 2012) was strongly associated with a migration index derived from the moth migration simulations (see Appendix S1), which ranged from À1 to +1, where negative values indicate higher contributions of the TEX natal origin and positive values indicate higher contributions of the FLA natal origin (log-likelihood v 2 = 60.2; df = 2; P = 0.001).

DISCUSSION
Simulated migratory flights were generated from short-term weather data, but resulting synoptic distributions revealed association with interannual climatic variability. Migratory flights were principally influenced by wind velocity (speed and direction) but were also directly linked with dynamic heterogeneity of emigrant sources driven by biological, meteorological, and behavioral factors. Incorporation of insect flight orientation could profoundly alter model simulations of individual flights (e.g., along an oceanic coastline) and generational migrations because flight orientation has been found to optimize seasonal migration displacements in several diurnal and nocturnal migratory insect species (Chapman et al. 2010).
There are several possible reasons for the limited fit between simulated immigration and capture of FAW in pheromone traps. One likely factor was practical limitations in monitoring ❖ www.esajournals.org range as years with more values across the full range of dates (2011 and 2014) showed a consistent positive relationship, whereas years with limited sampling (2012 and 2013) did not. Furthermore, pheromone traps attract only males at a stage receptive to a sex pheromone, which may not accurately reflect the timing, number, and geographic distribution of the total population. Methods that are not sex-biased and data with strong temporal and spatial resolutions would be ideal. There are also simplifying assumptions made in the model that could lead to inaccuracies. The current model did not incorporate flight speed or flight heading (crab angle relative to wind heading) and only incorporated field corn as a potential host plant. While corn is the primary host for the population of fall armyworm examined, significant contributions from other hosts are possible. It is also possible that  winter-breeding FAW populations along mild Gulf coastal areas may have occurred farther north than the assumed northern limit (i.e., 28.5°N latitude; Garcia et al. 2018). In these locations, initial FAW captures would be from locally developed populations and not be indicative of a migration event.
Coupled biophysical air transport-based migration models have the potential to identify regions at risk for pests capable of long-distance flight and, in combination with future climate scenarios, project the consequences of climate change on pest migration. This will be particularly relevant to risk assessments involving the entry and establishment of invasive pest flying insects and the dispersal of newly arising alleles in endogenous pest populations. This is particularly relevant given the recent establishment of fall armyworm in Africa and its rapid spread throughout that continent and in India (Day et al. 2017, Ganiger et al. 2018 and China (Li et al. 2019). Li et al. (2019) recently produced air transport-based model projections that supported the plausibility of annual fall armyworm migrations in Asia analogous to what is observed in North America. In this case, regions in Indochina and South China are projected to serve as overwintering sources for progressive waves of migration into eastern and northern Asia with the potential for substantial economic losses along the pathways. One notable distinction in their modeling is the assumption that fall armyworm will undergo three successive nocturnal flights, which differs from the variable number of flights, ranging from one to three, that are sufficient to describe the multigenerational migration observed in the United States. It will be interesting to see whether this difference is significant as it may indicate the conditions in Asia are less favorable for long-range migration than North America, which could mitigate the infestation risk. Efforts are underway in China to identify and introduce or augment natural enemy populations to reduce fall armyworm infestations (Silver 2019). Under such circumstances, migration models can inform conceptual designs for areawide or regional pest management programs by enabling estimates of the geographic potential of pest suppression that could be achieved with earlier pest suppression in smaller source regions. Specifically, interannual variability of flights and relative contributions of migrants from respective source areas may significantly impact the severity of pest infestations, and the mixing of genetic traits such as resistance to synthetic insecticides and insecticidal transgenic crops. T. Bailey assisted in trapping moths and reporting weekly moth capture data. USDC-NOAA provided access to HYSPLIT PC version 4.9 atmospheric transport modeling software. Mention of trade names or commercial products in this article is solely for the purposes of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture. USDA is an equal opportunity provider and employer. This material is based upon work that was supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture, under award number 2011-67003-30209. JW, SF, RM, and RN conceived the ideas and designed methodology; JW, SF, RM, and RN collected the data; JW, SF, SJ, and RN analyzed the data; SJ developed and implemented the model; and JW, SF, SJ, RM, and RN led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication. The authors declare no competing interests.