EPA/600/R-10/055 May 2010 Simulation Models Evaluation of Pest Resistance Development to Refuge in the Bag Concepts Related to Pioneer Submission Letter Report By Michael Caprio1, John A. Glaser2 ^isssissippi State University, Department of Entomology and Plant Pathology, Mississippi State, MS 39762-9775, 2US Environmental Protection Agency, National Risk Management Research Laboratory, Sustainable Technology Division, Cincinnati, OH 45268 January 22, 2010 Version 1.0 Sustainable Technology Division National Risk Management Research Laboratory US Environmental Protection Agency Cincinnati, OH 45268 Keywords: PIP corn variety, PIP crop registration, resistance management, simulation modeling, refuge in bag A portion of the work outlined in this report was accomplished under contract with Dr Michael Caprio ------- Introduction The USEPA, under its administration of the Federal Insecticide Fungicide and Rodenticide Act (FIFRA), requires the registration of all pesticides and pesticidal materials. The GM crops containing pesticidal traits are subject to FIFRA registration requirements. Concerns relating to future environmental effects of these crops can be investigated only using simulation models. Model application in this instance offers an analysis of the useful lifetime for the pesticidal traits. USEPA has clearly supported the use of these traits as replacements of broad spectrum pesticides that have significant environmental footprints. Background Resistance management (RM) simulation models are designed as deterministic, stochastic, and spatially explicit analytical methodologies. Simulation models can provide a realistic assessment of the risk of resistance evolution given the allele frequencies used as initial conditions. New concepts of PIP crop deployment such as refuge-in-the-bag require a detailed evaluation to ensure the claims associated with each concept offer the sustainable protection of the crop. Simulation models offer a means to evaluate the threat of crop loss through the development of resistance in the near future. Purpose This collaborative research effort was designed to evaluate the relative merit of refuge-in-bag concepts as submitted by Pioneer Hi-Bred International, Inc. This research specifically analyzed the relative risks related to the control of western corn rootworm in block refuge deployment and 5% seed- mixture scenarios. Research Design To understand the relative risks related to the evolution of resistance in western corn rootworm between the currently mandated block refuge deployment and a 5% seed-mixture deployment as submitted Pioneer Hi-Bred International, Inc., two different models were developed and evaluated. The two models consisted of a modification of a spatially-explicit stochastic model (Caprio et al. 2006) and a simpler, frequency-based deterministic model. The latter could be run in a single simulation mode with a graphical interface to enter parameters or in a risk-assessment mode capable of running thousands of simulations to estimate the effects of parameter uncertainty. Realistic estimates for ------- several parameters using in the deterministic model were developed through use of the stochastic model. Impacts and Outcomes This research was designed to provide substantive information about the operation and capabilities of different resistance management models to assist the regulatory expert in its proper use and interpretation of results. Models Employed EPAPOPGEN (Population Genetics Simulation Models) are a series of generalized, flexible population genetics modeling systems emphasizing the evolution of resistance to conventional pesticides as well as transgenic crops. The programs are capable of simulating complex landscapes that vary spatially as well as temporally over the course of a growing season. POPGEN-D The deterministic model POPGEN-D is a simple frequency-based model that incorporates non- compliance, non-random mating in compliant fields and larval movement. The model includes a switch to vary between a block-refuge design and a seed-mixture design. Biological parameters were primarily derived from that model with the exception of mortality rates, block size, dispersal rates and spatial aspects of refuge placement and movement. POPGEN-S2 The stochastic model POPGEN-S2 is a modification of a previous written stochastic, spatially explicit, dual locus simulation model designed to evaluate the population dynamics of crop pest insects relating to the development of resistance to pesticidal controls installed in the crop. (Caprio et al 2006). Materials and Methods Stochastic Model The stochastic model was a modification of a previous western corn rootworm model (Caprio et al 2006). Biological parameters were primarily derived from that model with the exception of mortality rates, developmental delays, block size and dispersal rates. These changes are discussed below. The stochastic model did not incorporate larval movement and was used only for evaluations of ------- different block refuge scenarios. Mortality rates: Mortality rates on transgenic corn were adopted from Pioneer's submission. We used two different dominance values (//), 0.05 and 0.2. If wss is the likelihood that a susceptible individual will emerge from transgenic corn relative to non-Bt corn, and we standardize Wn- as 1, then the fitness of heterozygotes (wrs) is estimated as (WIT- wss)*/7+wss. When /z=0.05, heterozygotes were ca. 5x more likely to survive than susceptibles, and were ca. 20x more likely to survive when h=0.2 when wss = 0.01. Those rates decrease as the fitness of susceptible homozygotes increase (a reflection of the difficulty in using dominance rather than the wrs:wss ratio). Developmental delays: To replicate the developmental delay of larvae on transgenic plants in the Pioneer model, we added the ability to delay insects by not moving a portion of any age group by a fraction. Normally, all the individuals in an age group would age one day in the model per time step. In the current model, 5% were held back from each age group (though they did experience mortality) and did not age. This resulted in a mean developmental delay of 2.3 days. Block size: The simulated universe consisted of a 40x40 matrix of 30m x 30m blocks, or 144 ha. Four different types of simulations were run with different field sizes. Because the current refuge strategy requires that refuges be associated with each field, the four different scenarios were field sizes of 5x5 blocks (each block is 30m x 30m, for a total of 2.25 ha), 10x10 blocks (9 ha field size), 20x20 blocks (36 ha field size) and finally a single 40x40 block field of 144 ha. Mean field sizes in the U.S. range from as low as 2.5 ha/field to as high as 150 ha/field, with most fields in the 30-50 ha/field range. In each case, a 20% refuge was planted along one edge of the field. With the 5x5 block fields, this refuge consisted of a single row of 5 blocks along one edge. In the 10x10 block fields, the refuge was 2 rows wide and 10 blocks long. The refuge widths for the remaining scenarios were, of course, 4 and 8 rows wide, maintaining in each case a 20% refuge (by area). Pioneer maintained that a square field is a best-case scenario (Flexner 2009), but we feel the opposite is true. It is difficult to determine from the Pioneer submission, but it appears that they simulated rectangular fields with refuges planted at the short ends of these fields. Most growers, however, plant lengthwise along such fields to minimize the number of turns they must make. It is most common to plant the refuge along one long edge of the field. If dispersal of adults is in question, then as fields become more rectangular (length:width ratio increased), refuges become thinner but the maximum distance between any Bt-plant and the refuge decreases. Under these conditions a square field (as we used) maximizes the distance between refuge and Bt plants and can be seen as a worst case. Pioneer's model (or at least my interpretation of their 4 ------- model), maximizes the distance between refuge and Bt plants, but is not consistent with typical grower practices. Adult dispersal: Dispersal was modeled based on Pioneer's dispersal data. We combined the data from Pioneer Fig. 3 and built probability distributions based on the frequency of beetles moving different distances. In short, if the total height of the bars in the above graphs summed to one, then the individual bar heights would represent the likelihood that a beetle would move that distance. One change was made. The Pioneer data includes beetles that did not disperse. The zero movement class consisted of 45% of the beetles, so the stochastic model assumes that 45% beetles do not disperse (on a daily basis) while the remainder move according to the Nowatski data. The terminology difference may be important. Pioneer's interpretation of the data is that beetles move an average of 14.5m/day. Using the same data, but only looking at beetles that move, the average movement rate is ca. 29m/day. This is purely a difference based on definition. In the our stochastic model, each day a beetle that was dispersing would be assigned a random position in the 30x30m block, chose a random direction, determine the distance it moved (based on the Nowatski data) and move it to a new block. As beetles disperse, they may move only a short distance and not change the 30x30m block where they are currently located. To determine how well the stochastic model dispersal kernal matched the Nowatski data, a simulation was run in which one block was populated with adults and the remaining blocks were empty. The number of beetles along a transect through the fields was counted every 3 days (Table 1, Fig. 1), corrected for the number of cells at a distance and multiplied by the distance they had dispersed to estimate the mean dispersal rate (per Nowatski) at different times (Table 2). The dispersal rate estimated after 3 days is 17.6m/day, which approximates the Nowatski estimate. Given the model uses actual field data to parametrize the dispersal, it was felt that this was a suitable match. Interestingly, the data from Table 2 demonstrate that the estimate of the mean dispersal rate drops over time. This occurs because beetles don't move in a straight line day after day. Brownian diffusion (based on a random walk) is usually assumed for insect dispersal. By confounding measurements over 3 days and using that data in a daily time-step model, Pioneer may have underestimated the daily dispersal rate of western corn rootworm. To demonstrate this process, a simple spreadsheet was constructed in which 10 individuals are simulated that all move 14.5m/day. At the start of each day a new direction is randomly drawn and the beetles move 14.5m in that direction. The right-hand column shows estimates of the dispersal rate using the Nowatski method. Only the estimates from the first day accurately ------- estimate dispersal rate/day (because all the individuals are moving away from the release point). Estimates based on beetles displacement on subsequent days tend to underestimate daily movement. After three days of movement the estimated dispersal rate/day is less then half the actual value (because the spreadsheet uses random numbers, values may differ). Unfortunately, because we do not have access to the original data and do not know the capture distributions for beetles of different ages, we cannot correct for this problem. In retrospect, analysis of the data using a diffusion model may have been more appropriate. If we assume that Table 3 can show the extent of underestimation of the dispersal, and that equal numbers of beetles were caught on days 1,2 and 3, a "ballpark" figure is indicates Nowatski underestimated the daily dispersal rate by a factor of 0.68. Using a definition of estimating the dispersal as the distance a beetle moves when it moves (this is purely semantics but critical to define), and correct by the above factor, the resultant dispersal rate is ca. 41m/day. It should be noted that, because a different definition of dispersal was used, this rate is not the same as the 45m/day dispersal rate used in the sensitivity analysis in the Pioneer submission. Refuge movement: Several scenarios were explored for yearly refuge deployment. One scenario, used for the transect experiments presented below, used refuges that were moved every generation. All growers planted their refuge on the west side of their fields in one generation and moved them to the east side in the next. Because the east side of one field is adjacent to the west side of a neighboring field, this scenario resulted in refuges always being planted adjacent to a refuge from the previous year. Ultimately, random placement of refuges each year was chosen as the most likely grower behavior. In this case, at the start of year, the simulation first determines if a field will be compliant, and if so, randomly chooses which side (east or west) of the field to plant the refuge. This means that occasionally refuges are not moved, while in other cases they are moved to relatively isolated areas. Deterministic model The deterministic model is a relatively simple frequency-based model that incorporates non- compliance, non-random mating in compliant fields and larval movement. It includes a switch to vary between a block-refuge design and a seed-mixture design. Appendix A includes descriptions of the different parameters. A schematic of the model is shown in fig. 2. The primary difference between the Bt population in the non-compliant population and the Bt-population in the non-random mating population is that no refuge was planted for the non-compliant fields. Non-compliant populations reduce both non-random mating and refuge size. As noted in the Pioneer submission, a critical ------- component of their modeling is that refuges are moved every generation. As mentioned before, while we believe that this may overestimate movement of refuges by growers, this was incorporated into the model by defining the refuge size not as the physical area devoted to the refuge, but rather the proportion of the eggs in a field that start out underneath the new refuge. We refer to this as the effective refuge size to distinguish it from physical refuge size (the proportion of an area planted to a refuge). These small refuge sizes are independent of non-compliance and result from growers moving their refuge to areas of the field were few eggs were deposited the previous year. To estimate this important parameter, we used the stochastic model to estimate egg densities at the start of a year across our four different field sizes (Fig 3a-e). In each case, we estimated a best and worst scenario. For the best case scenario, the grower moved their refuge to the opposite side of the field, next to where the refuge would have been planted in the neighboring field the previous year (these simulations were run with all growers planting refuges to the same sides of their fields and swapping them every year). The worst case scenario was based on planting the refuge in the area of the field where the egg density was the lowest. Because the egg density distribution across fields has never been measured, it is unlikely that growers could accomplish this, but it does represent and theoretical lower bound for effective refuge size. The worst case effective refuge sizes for the four different field sizes were (from smallest to largest field size): 15%, 11%, 3% and 0.055%. As noted by Pioneer, few insects are moving into the center of large fields as reflected by the effective refuge size of 0.055% for 144 ha fields. Other estimates for dispersal are higher and may depend on the region of the country. To test the effect of higher dispersal rates we doubled the vector of daily movement (to 29m/day) and repeated the large field size simulation (fig. le). In this case the density of eggs in the center of the field was fifty times greater under the worst case scenario and the effective refuge size was 1%. Most fields in the corn belt, however, are smaller, in the range of 30-50 ha, and we decided to use a lower limit of 1% for this parameter. Asymmetrical movement was included in the model, though the amount movement of susceptible larvae from transgenic corn is uncertain (Hibbard et al 2005). Asymmetrical movement was defined as the probability that larvae that have survived exposure on a Bt-plant will move. When larvae move, they move onto plant genotypes according to the ratio of those genotypes in the field. In this case, a larvae leaving a transgenic plant had a 95% probability of moving onto another transgenic plant. As with other rates in the model, this movement rate is converted into a daily rate (depending on the days of larval movement parameter). ------- To incorporate our uncertainty in the parameters used in the deterministic model, we developed distributions based on data for the lowest and maximum possible values we felt the parameter might take, as well as a most likely value. These values were used in PERT distributions based on a beta distribution with the most likely value weighted by 4 (Vose 2001). This was done for both the block refuge model (Table 4) and larval movement model (Table 5). Each model was run 2000 times, each time drawing parameter values from the relevant PERT distributions and using the sampled parameters as inputs to the model. The time it took for the resistance allele to exceed 50% was recorded for each simulation. We measured the mean and 25th and 75th percentiles of the time to resistance. We determined the proportion of simulations in which the time to resistance exceeded 10 years (i.e., 11 years or longer). At a cutoff point of 11 years, the gene has already been in use for 4 years and previous EPA-SAPs have suggested 15 years as a reasonable minimal goal for the lifetime of single gene products. Uncertainty distributions for the deterministic model. The uncertainty distributions are divided into those that are common to both models (Tables 4- 5), those specific to the block refuge model (Table 4) and those specific to the seed-blend model (Table 5). Common parameters 1. Fitness of susceptible homozygotes on transgenic plants (wss). The Pioneer estimate for this parameter was 0.0125. In the past, however, as transgenic events have been exposed to a wide variety of environments, wss has often been found to exceed initial estimates. The minimum value used was 0.0125, the maximum value was 0.08 and the most likely value was 0.04. 2. Dominance (h). Occasionally, dominance associated with events that are not a true high-dose is greater than the Pioneer estimate. Our estimate was centered on the Pioneer estimate but included the possibility that it could be 5x less or 5x greater. 3. Initial resistance allele frequency (IRF). While many estimates of the initial resistance allele frequency for loci providing resistance to Bt toxins are in the range of 0.001, recent studies with western corn rootworm have suggested they may be more common, perhaps in the 0.02 range. These two values were used as our maximum and minimum values, with 0.005 used as the most likely value. Block refuge parameters 1. Effective refuge size. Without non-compliance, the physical refuge size was always assumed to be 20% of the area. Non-random oviposition and refuge placement could, however, lead to that 20% of the area having more or less than 20% of eggs initially (see Figs 3a-e). This is a complex parameter that depends on the size and shape of fields, how growers physically place the refuges in the fields, how they move the refuges ------- between years and on corn rootworm dispersal rates. It is likely that this parameter is responsible for many of the differences between the results of this model and the model used by Pioneer. The effective refuge size varied from 1% to 20%, with 15% as the most likely value. 2. Non-compliance. Non-compliance must be measured on a field by field basis, and EPA's interpretation that a non-compliance rate of 30% is an overestimate seems reasonable. There may be localized areas, however, where non-compliance may be greater than 30%, and to incorporate the effects of these areas the maximum value of non-compliance was set to 60%. The minimum value was 5% with a most likely value of 10%. 3. Non-compliance mixing. This parameter is related to the rate of mixing per generation between compliant and non-compliant fields and may be seen as a measure of persistent non-compliance versus growers that an accidentally non-compliant in one or more fields. Different monitoring groups suggested 5-7% of growers were persistently non- compliant, so 93% was chosen as a most likely value, 95% as the maximum and 85% as a conservative minimum value. 4. Non-random mating. This parameter is an estimate of the proportion of the population emerging in the transgenic fields that mates among themselves with no contribution from refuge insects. Based on the simulated dispersal experiment, it appears as though as many as 50% of these individuals would be unlikely to mate with refuge individuals in the largest field sizes. Because refuges are moved each year, these non-random populations are merged with the random mating populations for reproduction. The maximum value used was 50%, with a minimum of 5% and a most likely value of 30%. Seed-blend model 1. Base larval movement. Data suggest there is limited larval movement among non- transgenic plants in the absence on high larval densities (which we find unlikely in a 5% seed-mixture). The range for this parameter was from 1%, to 5%, with a most likely value of 3% 2. Days of larval movement. Based on Mallet and Porter (1992), larval movement effects the evolution of resistance by shifting dominance. This shift occurs gradually over a number of days. The range of time over which larval movement occurred ranged from 1 to 18 days, with a most likely value of 10 days. The mortality rate was assumed to be constant over this period, with the overall mortality rate equal to 1 - w, the relative fitness of the genotype. 3. Asymmetrical movement. Food choice bioassys in many species suggest that larvae will tend to taste and then move from transgenic tissue at a greater rate than non-Bt tissue. Hibbard et al (2005) found, in a two year study, that this rate varied in those two seasons between 0 and 30%. Because of the limited data on this parameter, the range of this parameter was 0 to 50%, with a most likely value of 30%. Two locus deterministic model One final set of resistance simulations was conducted for two locus block refuge and seed- mixture models. In this case the fitness of the susceptibles and the dominance of resistance at each of two loci were independently determined according to the PERT distributions from the single locus model. This distribution is informative because it tells us the relative longevity of two genes stacked ------- into the same plant versus using the same two genes sequentially in single gene plants. This is relevant because introducing a single gene event into a landscape with dual gene events may lead to rapid resistance at the first locus and essentially turn all dual gene events that share the toxin into single gene events. Results and Discussion Deterministic Model The distributions of the time to resistance for the block refuge strategy and the larval movement strategy with and without asymmetrical movement are presented in figure 4. Because the bars each represent a year, the density of the bars is the proportion of simulations that resulted in resistance evolving in that generation. The mean time to resistance for the block refuge strategy was 11.95 years, and 56.55% of the simulations exceeded 10 years. For the larval movement model, the mean time to resistance was 8.1 years and 13.5% of the simulations exceeded 10 years. The results from the two locus block refuge deterministic model (Fig. 5a) show that, though the second toxin has identical uncertainty distributions as the first toxin, the mean expected lifetime of the dual gene product is 132 generations in a block refuge deployment, 6 fold longer than if the two toxins had been used sequentially in a block refuge design. Similarly, for the two locus seed-mixture model, the dual gene product lasted a mean of 52.6 years (Fig. 5b), 3.25 fold longer than if the two toxins had been used sequentially. Single gene products that share toxins with dual gene events can lead to a pseudo sequential use because the single gene event selects rapidly for resistance to one gene, and once this occurs the dual gene events are effectively rendered single gene events. Stochastic Model The stochastic model was used to evaluate 4 different sizes and two different dispersal rates (Fig. 6). Because the largest field size took the entire simulated universe, it was not possible to incorporate non-compliance. To address this issue, the model was run in an 80x80 cell mode (6400 30x30m cells) for the largest field size, allowing simulation of four 144 ha fields. These larger simulations for the 144 ha field size were used in Fig. 6 and in the data analysis. Analysis of variance showed no significant impact of dispersal rate, but field size was highly significant (p « 0.001). When the simulation was run with a single 144 ha field and no non-compliance, the times to resistance were much longer (mean = 3292 d). It was not expected that doubling the dispersal rate would have so little 10 ------- impact on the time to resistance, and it suggests the model is relatively insensitive to this parameter. Using Pioneer's 40 ha field size as a benchmark, the stochastic model predicted an increase in durability of 3.94 fold compared to simulations with no refuge (see discussion). This is much higher than the 1.4 fold reported for this field size in the Pioneer simulations. The simulations presented predict that a block refuge system with 36 ha fields would last 13.9 y, almost double the time predicted by Pioneer. The 144 ha simulations showed much lower durability (8.4 y, or 2.4 fold longer than without a refuge). It may be relevant that the dimensions of these fields had a 240m wide refuge and a 960 meter wide strip of Bt-corn. Pioneer's simulations may have been run with a 200m wide refuge and an 800 wide strip of Bt-corn (this occurs in a 40 ha field size in the Pioneer model but a 144 ha field in the current model). To determine further the sensitivity of the model to the rate of dispersal, the 36 ha field size was run with 2x, 4x and 8x dispersal. Dispersal was again not significant in an ANOVA (p=0.058). The 36 ha field size has only 4 refuges in the default 144 ha system, and this increases variance when in any given year each field has a 30% of being non-compliant. There is a 24% chance that no fields will be non-compliant in any year and a 0.8% that all four fields might be non-compliant. In a larger system the mean number of non-compliant fields would remain more constant over time (Law of Larger Numbers) and there might be a better opportunity to investigate how dispersal interacts with refuge size and placement. It is possible that dispersal has less effect than might have been expected because there are several different situations that arise in the model because of yearly random refuge placement and non-compliance. When a refuge is not moved, lower dispersal rates result in a larger effective refuge. Simulations suggest this frequently delays resistance evolution even though it reduces the probability that an individual beetle will move and mate at random with beetles emerging in Bt-fields (Caprio 2001). Similarly, if the refuge is moved next to the location of where a refuge had been planted in a neighboring field in the previous year, lower rates of dispersal might also maximize effective refuge size. In contrast, when refuges are moved to areas relatively distant from the locations of the previous years refuges, higher rates of dispersal would maximize the effective refuge size. It is clear, however, that these results suggest that the model is less sensitive than expected with respect to dispersal rates. Discussions over whether the Nowatski data over or underestimate movement rates for the entire corn belt are, perhaps, less critical than expected. Discussion 11 ------- Pioneer presents in Table 1 (MRID PHI-2009-197) the results from their model using the benchmark parameters listed in Table 2. Using the deterministic model presented here with similar parameter values, similar but not identical durabilities were found. In the case of no refuge, Pioneer's model suggested resistance would evolve in 5 years, while this model suggested durability would be 3 years. The stochastic model predicted a durability of 3.5y in the absence of refuges. For comparison's sake, we measure other rates relative to these baseline numbers. For the block refuge model, Pioneer suggested resistance might evolve in 7 years, or 1.4x the time without a refuge. Assuming our most likely value of a 15% effective refuge (and otherwise using Pioneer's benchmark values), we obtain a predicted durability of 14 years, 4.7 fold longer than the simulation without a refuge. One possible reason for the large discrepancy is the effective refuge size. To duplicate the 1.4x durability increase noted by Pioneer, the effective refuge size would need to be decreased to 1%. This was the minimum value in our uncertainty distribution for this parameter. This parameter was set as the minimum because it would require large fields (on the order of 150 ha), that the dispersal rates used by Pioneer were not underestimated and are typical of rootworm dispersal rates throughout the corn belt (there is some data to suggest otherwise), and that growers consistently move refuges to regions of their fields with the lowest egg density. The overall risk profile from our uncertainty simulations for block refuges is impacted by other parameters that deviated from the Pioneer benchmark values as well. Although relative fitness, initial frequency and dominance used the same uncertainty distributions for both the block and seed-blend simulations are were likely to affect both models in a similar manner, these parameters could shift both distribution in time. The inclusion of the possibility of higher levels of dominance, in particular, could lead to shorter predicted durability estimates. Other parameters specific to the block refuge model could lead to the longer durability predicted in these simulations compared to Pioneer's simulations. The distribution for non-compliance was weighted more heavily towards lower values (with a most likely value of 10% non-compliance) than reported by CSPI (Center for Science in the Public Interest), based on EPA's analysis that report overestimated non-compliance when measured on a field by field basis. We assumed that most non-compliant growers were "accidentally" non- compliant, perhaps due to miscalculations or mistakes during planting. Only 5% of non-compliant growers were assumed to be persistently non-compliant year after year. For the larval movement model, the current model predicted a durability of 9 years, or a 3-fold increase in durability compared to no refuge. This compares reasonably to the Pioneer estimate of 10 years, or a 2-fold benefit over their results run without a refuge. The values both compare quite well to the mean of our risk profile, 8.1 years. 12 ------- These data suggest there are minor differences between the larval movement results generated by the Pioneer model and the model presented here, but the major differences lie in the block refuge model. Although Pioneer suggests that non-random mating is important in their results, we did not find that to be the case in our simulations (Table 6). The parameter most likely to be responsible for the differences between the two models is the effective refuge size. We agree with the Pioneer assertion that it is possible that areas of large fields could have very low egg densities, and should growers manage to locate and plant refuges over those areas, the block refuge strategy could be severely compromised. Understanding of the differences in the results is limited by details missing in the Pioneer data submission. They simulated a size of 4x10 cells (each cell as 100m xlOOm), or 40 ha. It was unclear if this is a single field or represented multiple smaller fields, but given their results, a single field is likely. The exact placement of the refuge in the 4x10 cell universe is important and not stated. Based on an ESA presentation (Flexner 2009), it appears as though the refuge was located to one "thin" edge of this system (perhaps 2 4 cell rows to one edge), which would maximize the furthest distance between transgenic plants and the refuge, but would seem to be at odds with common grower behavior (this is supported by the fact that 20% of the cells by planting along the long edge in this size is not possible in this geometry). Again, based on the same ESA presentation (Flexner 2009), it appears as though the only way for insects to reach the plants to the right of the system (when the refuge was planted on the left side) was to move across the intervening 800m. In actual field systems there are other fields surrounding the simulated field with refuges that, on average, would be adjacent to the current field, and the mean distance between the furthest from a refuge and the refuge would be 400m. Most simulations are designed as toruses (an insect moving off the left side of the system arrives at the right side, as if the system were a donut) to overcome this problem, but the animation (Flexner 2009) did not appear to be designed in this manner. This model minimizes the effects of dispersal and would tend to lower Pioneer's estimate of the durability of a block refuge system. In the current stochastic model, we simulated multiple fields arranged in a large torus. Refuges were placed randomly along one long edge of fields, and non-compliant fields were also chosen randomly each year. In these simulations, occasionally refuges in adjacent fields would be side by side, while in other cases they would lie at opposite ends. In the latter case there would clearly be an area with low egg densities between the two refuges, but that is apparently offset by refuges that must logically be more closely spaced in other areas. The Pioneer model and our dispersal kernel based on the Pioneer dispersal data ignore the potential of longer distance dispersal (Storer et al 2003). Although the contribution from 13 ------- such moths is likely to be small, in areas of fields were there are essentially no eggs (see Fig 3d), the contribution of these few individuals could be significant. A number of explanations for the differences in the predicted longevity of the block refuge strategy between the Pioneer simulations and those presented here are possible. Important details of the Pioneer model are missing and it is difficult to determine how certain key elements were modeled. The Pioneer model used a single 40 ha block of lOOxlOOm cells arranged in a rectangular fashion (4x10). In discussions, corn entomologists expected most growers would plant lengthwise in such a field, and the refuge would be planted along either one of the long edges (here termed the east and west edges). Simulations of refuges planted to the short ends of these fields (the north and south edges) would not seem to be the normal grower practice and would tend to increase the mean distance between refuges and Bt plants. This choice, if made, would increase the likelihood of the holes in random mating mentioned by Pioneer in their submission and reduce predicted durability of the block refuge design. Pioneer provided little detail in the spatial aspects of their simulation. Was a single field simulated with all the refuge gathered in one place, or were smaller fields simulated with multiple refuges? When the refuge was moved each, where was it moved to? Was it moved to a random location, or was it always moved to the opposite end of the field from it's previous location? The latter choice might lead to a significant decrease in the effective refuge size and more rapid resistance development. How were boundary conditions treated? When an insect moved from the field to the west side, for example, what happened to it? Did it die, was it reflected back into the system, or was a torus simulated? The choice here may be important as it can again impact the mean distance between refuge and Bt plants and ultimately effective refuge size. As an example, consider a string of 4x10 fields, each 1000m long, each with it's own 200m wide refuge. The maximum distance between any two refuges would be 1600m, and the maximum distance between a refuge and a Bt plant would be 800m. The average distance would be 800m between refuges and 400m between Bt plant furthest from a refuge and the nearest refuge. If Pioneer chose not to use a torus design, and simulated a single field with a refuge at either the west or east end, the maximum distance between a Bt plant and the nearest refuge would be doubled (800m). In a simulation of a torus, insects leaving the east edge of the field end up in the west edge, in essence creating a long string of identical fields. This design correctly models the 14 ------- mean distance, but there is no variance. The simulations described here included multiple fields with random placement of refuges on a yearly basis to either the east or west edge of a field, which seems to best mimic grower behavior. This design meant that in some cases the refuges would be at opposite ends of adjacent fields, while in other cases they might be side by side. This design choice allowed simulation of both the correct mean distance between refuges, but the variance in this distribution. Our interpretation of the data indicates that rather than non-random mating being the threat to a block refuge system, it is non-random oviposition and refuge placement that are most important. The conclusion is supported the results of the stochastic model run with just one 144 ha field. The time to resistance in this case was over five years longer than any other scenario, despite the fact that substantial areas of non-random mating must have existed in the center of these blocks (see Fig. 3d). What did occur is the field was always moved adjacent to a previous years refuge, so there was a relatively high density of eggs under the new refuge and the effective refuge size was high. This suggests a potentially plausible method that Pioneer might use to bolster their claims that failure of the block refuge system is imminent. If the relative density of eggs could be estimated in transects across fields of different sizes (similar to Fig. 3a-e), and if more accurate estimates of how growers plant and move their refuges over years, then major differences between the two models might be resolved. Such a study might provide an estimate of the contribution of beetles dispersing long distance distances. 15 ------- References Caprio, M A, T Nowatzki, BD Siegfried, LJ Meinke, RJ Wright, and LD Chandler, 2006, Assessing Risk of Resistance to Aerial Applications of Methyl-Parathion in Western Corn Rootworm (Coleoptera: Chrysomelidae), J. Econ. Entomol. 99(2): 483-493. Caprio, M. A., D. Sumerford & S. Sims. 2006b. Evaluating transgenic plants for suitability in pest and resistance management programs, 2nd edition. L. Lacey & H. Kaya [eds]. Field Manual of Techniques for the Application and Evaluation of Entomopathogens. Hibbard, BE, MI Higdon, DP Duran, YM Schweikert, and MR Ellersieck. 2004. Role of egg density on establishment and plant-to-plant movement by western corn rootworm larvae (Coleoptera: Chrysomelidae). J. Econ. Entomol. 97(3): 871-882 Hibbard, BE, TT Vaughn, IO Oyediran, TL Clark, MR Ellersieck. 2005. Effect of Cry3Bbl expressing transgenic corn on plant-to-plant movement by western corn rootworm larvae (Coleoptera: Chrysomeldiae). J. Econ. Entomol. 98(4): 1126-1138 Mallet, J, P. Porter. 1992 Preventing Insect Adaptation to Insect-Resistant Crops: Are Seed Mixtures or Refugia the Best Strategy?, Proc R Soc London B 250, 165-169. Meihls, LN, ML Higdon, BD Siegfried, TA Spencer, NJ Miller, TW Sappington, MR Ellersieck, and BE Hibbard. 2008. Increased survival of western corn rootworm on transgenic corn within three generations of on-plant greenhouse selection. Proc.Nat. Acad. Sci. 105: 19177-19182. Pan, Z., DW Onstad, BH Stanley, JL Flexner, Q Chen, and SALefko, Addendum to: A Spatially Explicit Computer Simulation Study of Corn Rootworm, Diabrotica spp. [Coleoptera: Chrysomelidae], Pioneer Hi-Bred International, Inc., Study ID: PHI-2009-125 Rood, T, Z Pan, D Onstad, B Stanley, JL Flexner, Q Chen, T Nowatzki, M Ziegler, and S Lefko, Additional Addendum to MRID 473567-08: A Spatially Explicit Computer Simulation Study of Corn Rootworm, Diabrotica spp. [Coleoptera: Chrysomelidae], Pioneer Hi-Bred International, Inc., Study Number: PHI-2009-197 Vose, 2000. Risk Analysis: A quantitative guide. Wiley, New York. 16 ------- Beetle Position By Day After Release JD M— O JD E 5000 500 50 0.5 0 30 60 90 120 150 180 210 240 270 300 330 360 Distance (m) Figure 1 Beetle Dispersal Time Series 17 ------- Non-compliant Bt Bt Compliant Refu< Figure 2. Schematic of the deterministic frequency based model for both seed-blend and block refuge deployments. 18 ------- 6000 Figure 3a Small 2.25 hectare field size with a 30m row of refuge along the edge of a 120m wide field oftransgenic corn. The refuge is always planted along an edge and is switched every year. Previous year's refuge was in strip 5, while the current one is in strip 1. Growers always plant to the same side of their fields, so the reason there is a larger egg population under strip 1 than in strip two is because there was a refuge to the left in another field. In this case, if growers planted in strip one, the new refuge would get 19.8% of the overwintering eggs. However, if they planted in the area with the lowest density, the refuge would 15% of the overwintering eggs 25000 7 8 10 Figure 3b 9 hectare field. The highest density is under the previous years refuge, but continued "drift" from neighboring fields is observed. In this case 18.6% of the eggs would be under the refuge in the default planting pattern, but the refuge would only get 11% if growers planted the refuge in the area with the lowest density of eggs. The refuge in this case was 2 strips (60m) wide. 19 ------- 40000 35000 30000 25000 20000 15000 10000 5000 10 11 12 13 14 15 16 17 18 19 20 Figure 3c 36 hectare field (about an average field for many parts of the U.S.). With the default planting pattern, the refuge would receive 20% of the eggs, but only 3% if growers planted over the area with the lowest density. The refuge in this case was 4 strips wide (120m). 25000 20000 15000 10000 5000 111 .ll 1 23456789 10111213141516171819202122232425262728293031323334353637383940 Figure 3d 144 hectare field size. At this size, the default planting would result in 23% of the eggs in the new refuge, but the worst case drops to 0.055% (if growers plant in the middle). The refuge is 8 strips wide (240m). 20 ------- 40000 35000 30000 25000 20000 15000 10000 5000 40 112 184 256 328 400 472 544 616 688 760 832 904 976 10481120119212641336140814801552 4 76 148 220 292 364 436 508 580 652 724 796 868 940 101210841156122813001372144415161588 Figure 3e Doubled dispersal distance in a 144 hectare field. This figure contains a column for each 30x30m block rather than an average of 3 transects as in the other four figures. Figure 3. Simulated egg densities in transects across Bt-corn and associated refuges. 21 ------- S - Generations Figure 4a. Block refuge model ~ -, O.I ' TE !^r OOOOOOOOOO 20 ne ratio ns Figure 4b. Seed-mixture model Figure 4. The distributions of the time to resistance 22 ------- 100 200 Dataset$time 300 400 Figure 5. The distribution of time to resistance from 2000 simulations of a block refuge deployment with two resistance loci. The mean time is 132 generations and 100% of the simulations persisted beyond 26 years. 23 ------- Plot of Means o o 10 LU CO F O o in o o o 225 9 FieldSizec-$ D ispe rsa 36 Field Size 144 Figure 6. Means of four simulations at each field size (given in hectares). Simulations used for the dark line used the Nowatski data, simulations used in the red line doubled that dispersal rate 24 ------- Days cells 0 4604 675 224 110 58 44 20 1 30 542 386 169 103 61 31 26 8 60 119 182 99 60 39 31 20 16 90 20 64 72 57 35 26 13 24 120 2 19 37 30 25 28 13 32 150 3 12 15 17 9 18 9 40 Distance 180 0 4 9 13 7 5 4 48 210 0 3 2 7 5 0 3 56 240 0 3 0 0 2 1 1 64 270 0 2 0 0 2 4 1 72 300 0 0 4 1 0 3 1 80 330 0 0 0 0 0 1 0 88 360 0 0 0 1 1 1 0 96 Table 1. Beetle counts over time in a simulated dispersal experiment. Distance is the distance in meters from the initial release point. Days 0 0 0 0 0 0 0 0 30 130080 92640 40560 24720 14640 7440 6240 60 114240 174720 95040 57600 37440 29760 19200 90 43200 138240 155520 123120 75600 56160 28080 120 7680 72960 142080 115200 96000 107520 49920 Distance 150 180 18000 72000 90000 102000 54000 108000 54000 0 34560 77760 112320 60480 43200 34560 210 0 35280 23520 82320 58800 0 35280 240 0 46080 0 0 30720 15360 15360 270 0 38880 0 0 38880 77760 19440 300 0 0 96000 24000 0 72000 24000 330 0 0 0 0 0 29040 0 360 0 0 0 34560 34560 34560 0 Sum 313200H 705360 720480 675840 50112oM 580800 286080 n/day 27.22 17.64 13.66 11.09 9.28 8.98 6.81 Table 2. Estimate of dispersal rates over time in a simulated dispersal experiment. Each cell is the product of the corrected number of beetles at a distance in time multiplied by the distance. These are summed, averaged by beetle numbers and divided by the number of days for the sample to obtain the mean displacement in m/day. 25 ------- Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 m/day X-position Y-position direction displacement x displacement y total displacement direction displacement x displacement y total displacement direction displacement x displacement y direction displacement x displacement y total displacement 0 0 Day 1 2.43 -10.94 9.51 Distance 14.5 Day 2 1.95 -16.38 22.96 Distance 28.2 Day 3 0.99 -8.48 35.12 Distance 36.13 Day 4 3.05 -22.92 36.38 Distance 43 0 0 Day 1 1.06 7.03 12.68 Distance 14.5 Day 2 5.11 12.7 -0.67 Distance 12.71 Day 3 1.43 14.71 13.69 Distance 20.1 Day 4 1.13 20.84 26.83 Distance 33.98 0 0 Day 1 3.26 -14.39 -1.77 Distance 14.5 Day 2 0.28 -0.45 2.23 Distance 2.28 Day 3 6.06 13.67 -1.03 Distance 13.71 Day 4 5.24 20.92 -13.59 Distance 24.95 0 0 Day 1 3.6 -13.03 -6.36 Distance 14.5 Day 2 1.48 -11.77 8.09 Distance 14.28 Day 3 1.88 -16.25 21.88 Distance 27.25 Day 4 1.71 -18.27 36.24 Distance 40.58 0 0 Day 1 2.08 -7.03 12.68 Distance 14.5 Day 2 2.13 -14.71 24.98 Distance 28.99 Day 3 5.86 -1.47 19.08 Distance 19.14 Day 4 1.73 -3.74 33.4 Distance 33.61 0 0 Day 1 3.12 -14.5 0.25 Distance 14.5 Day 2 4.59 -16.26 -14.14 Distance 21.55 Day 3 5.41 -6.94 -25.25 Distance 26.18 Day 4 5.9 6.5 -30.68 Distance 31.36 0 0 Day 1 5.85 13.14 -6.13 Distance 14.5 Day 2 5.52 23.57 -16.2 Distance 28.6 Day 3 6.25 38.06 -16.71 Distance 41.57 Day 4 0.45 51.1 -10.35 Distance 52.13 0 0 Day 1 0.52 12.56 7.25 Distance 14.5 Day 2 1.24 17.28 20.96 Distance 27.16 Day 3 1.48 18.54 35.4 Distance 39.97 Day 4 0.26 32.55 39.16 Distance 50.92 0 0 Day 1 1.08 6.81 12.8 Distance 14.5 Day 2 4.42 2.57 -1.06 Distance 2.78 Day 3 0.91 11.5 10.36 Distance 15.48 Day 4 5.69 23.52 2.25 Distance 23.62 0 0 Day 1 2.01 -6.13 13.14 Distance 14.5 Day 2 3.89 -16.73 3.25 Distance 17.05 Day 3 4.29 -22.63 -9.99 Distance 24.74 Day 4 1.85 -26.63 3.94 Distance 26.92 14.5 9.18 9.03 Table 3. Example of estimate daily dispersal rate over time from insects moving in a random walk at 14.5m/day. Direction is given in radians. The mean displacement (m/day) is an average from 10 individuals and will chance depending on the random directions drawn. 26 ------- Min Most Likely Max Wss 0.01250 0.04000 0.08000 Dominance 0.01000 0.05000 0.20000 IRF 0.00100 0.00500 0.02000 Effective refuge 0.01000 0.15000 0.20000 Non-compliance 0.05000 0.10000 0.60000 Non-compliance mixing 0.85000 0.93000 0.95000 Non-random mating 0.05000 0.30000 0.50000 Table 4. Uncertainty parameters used for the block refuge simulations Min Most Likely Max Wss 0.01250 0.04000 0.08000 Dominance 0.01000 0.05000 0.20000 IRF 0.00100 0.00500 0.02000 Refuge 5.00000 5.00000 5.00000 Base larval movement 0.01000 0.03000 0.05000 Days of larval movement 1.00000 10.00000 18.00000 Asymmetrical movement 0.00000 0.30000 0.50000 Table 5. Uncertainty distributions used for the seed-mixture simulations 27 ------- Coefficients: Estimate Std. Error lvalue Pr(>|t|) (Intercept) 11.6958 2.2321 5.240 1.78e-07 *** dominance -84.9810 1.3097 -64.887 <2e-16*** IRF -629.7121 14.7731 -42.626 <2e-16*** Non. compliance -8.8690 0.4426 -20.037 <2e-16*** non.compliance.mixing 1.5415 2.4028 0.642 0.521 non.random.mating -0.7173 0.5336 -1.344 0.179 Refuge 54.6024 1.2589 43.374 <2e-16*** wss 61.5433 3.4862 17.653 <2e-16*** Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1'' 1 Residual standard error: 1.797 on 1992 degrees of freedom Multiple R-squared: 0.8095, Adjusted R-squared: 0.8088 F-statistic: 1209 on 7 and 1992 DF, p-value: < 2.2e-16 Table 6. Regression coefficients for the regression on the results from the block refuge model. Coefficients: Estimate Std. Error lvalue (Intercept) 13.25428 0.20471 64.748 <2e-16 *** BaseLarvalMovement 7.43119 4.18984 1.774 0.0763 . dominance -63.84439 0.91981 -69.410 <2e-16 *** IRF -466.03593 10.13300 -45.992 <2e-16 *** SelectionDays 0.09145 0.01017 8.992 <2e-16 *** wss 57.75906 2.49046 23.192 <2e-16 *** Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1'' 1 Table 7. Regression coefficients for the results from the seed-mixture model with asymmetrical movement. 28 ------- Appendix A. Deterministic parameter description Wss - The fitness of homozygous susceptible genotypes on the transgenic plants. Pioneer used a value of 0.0125. Without better knowledge of the Pioneer field tests it's difficult to develop a range, but a specific factor over and under was used to deve;op our estimate. The sorts of questions to be asked might be what is the lowest/highest performance under field conditions. Instead of 98.75% mortality, could it be 95% (as an example)? Dominance - The relative fitness of heterozygotes. I use the traditional population genetics equation: Wrs = (Wrr - Wss) * h + Wss where h is the dominance and Wrr = 1.0. Our distribution is derived from Pioneer's estimate of h=0.05, but that seems low compared to many estimates based on LCSOs. LCSOs, however, are based on a log scale, however this research has used a linear scale. In general, for products with lower effectiveness, dominance tends to be higher (see Caprio et al. 2006b). IRF - The initial resistance allele frequency. My impression from dealing with CRW folks is that this value has tended to be higher in empirical estimates than we have seen for other insects. I believe I used Pioneer's default value as my most likely value. Block refuge model parameterization: Values cited above are shared with the larval movement model and changes in would effect both the block and larval movement models to a similar degree. Refuge proportion - Initially set this to 20%. However, it can be viewed as the proportion of eggs in a compliant field that end up being in a refuge (perhaps after being moved), and the non-random oviposition patterns can be addressed using this parameter. If more individuals survive in refuges, and they tend to lay eggs near that area, then the number of eggs in the soil to start the new growing season will be greatest near the refuge from the previous year. If the refuge is moved, it stands to reason the new refuge, even if it covers 20% of the area, will receive less than 20% of the eggs oviposited in the field. Based on this modeling research, this fact and the likelihood that refuges are moved are identified as key assumptions of the Pioneer model. Unfortunately, there is little empirical data on the non-random oviposition patterns of CRW. Non-compliance - This is the proportion of fields that are planted without a refuge. Individuals emerging from this area mate as a separate population. The most likely value is a literature value (CSPI). This parameter should be based on a field by field basis and incorporate some estimate of amount of non-compliance. A person planting a 3% refuge in a field is not the same as someone not planting any refuge. Alternatively, it shouldn't measure if a field is compliant (as yes./no answer), but how close to compliant a non-compliant field is. General comments have been encountered which indicated that 30% non-compliance was too low. It is not obvious how non-compliance was estimated. This parameter may be important because it determines how much Bt-corn is planted without any refuge. Non-compliance mixing - This determines how much movement there is between compliant and non- compliant fields. If non-compliant growers are accidentally non-compliant (i.e. they never got around to planting the refuge because of weather), then they will generally be compliant in the next year, and this parameter should be close to 1.0. On the other hand, if non-compliant growers are growers who 29 ------- purposefully don't plant refuges and continue this practice year-by-year, then this parameter should be low. How low would depend on how much movement you think naturally occurs between non- compliant and compliant fields. Inevitably this leads to discussion of the clustering of non-compliant fields. Matt Carrol has done some work on this with the Monsanto model. Non-random mating - This parameter addresses the proportion of the population in the Bt-portion of compliant fields that is far enough away from the refuges that they do not mate at random. These differ from non-compliant fields in that refuge is planted, but it was too far away. The non-random oviposition that results from this non-random mating (and limited dispersal capability) is covered by the reductions in refuge size discussed above. Larval movement model parameters The refuge proportion is fixed at 0.05 in the larval movement model. This was based on Pioneer's last data that indicated they had reduced the variance considerably. Base larval movement - This is the amount of movement that might occur in a non-transgenic field with populations low enough so that significant density-dependent dispersal is not a factor. In the stochastic models, population decline is found to be a greater problem than density-dependent factors such as mortality or dispersal. The model does include asymmetrical movement, with the movement of individuals on Bt plants increased by an amount related to their daily mortality rates. Susceptibles therefore move more than heterozygotes, and resistant homozygotes would only move at the base larval movement rate. Of course, an individual leaving a transgenic plant only has a 5% change of moving onto a non-Bt plant. Various estimates of larval dispersal rates have been encountered but mostly on the low side. Generally, this research assumes that populations will be low in seed-mix fields and there will be little to no density-dependent larval dispersal. This could be an area we need help with, or we might find by selecting alternate numbers is relatively unimportant compared to other parameters and a precise estimate would be a misapplication of resources. Days of larval movement - Movement and larval development are assumed to occur over this period of days. Movement occurs once per day. Daily survivorship is constant and is determined as the rate that raised to the days power results in the overall fitness listed for the genotype in the original table. In the case of Pioneer, that was 0.0125 for susceptible larvae. If the days of larval movement is 10, then the daily survivorship rate is 0.645. This means that on day 1, 36.5% of the SS larvae on Bt die (I assume conservatively that mortality occurs before dispersal) , and of the remaining 64.5%, 36.5% will disperse (see above). Of these dispersants, only 5% will move onto non-Bt plants, increasing the refuge size from 5% to ca. 6%. If the base larval movement is low (remember it is spread over days as well), then this process can lead to an increase in the refuge population and actually decrease the rate of resistance evolution. I put this into the model just to see what effect it might have. The value was varied between 5, 10 and 15 and it had very little effect on the results. 30 ------- |