Modeling fall armyworm resistance in Bt-maize areas during crop and off-seasons

Entomologists have often used computational modeling to study the dynamics of insects in agricultural landscapes. Recently, important issues such as the movement of adults and immatures associated with insect resistance to GMO (genetically modified organism) crops have been addressed using computational models. We developed an individual-based model using the cellular automata approach (CA) to investigate how an intercropping system composed of maize engineered with Bacillus thuringiensis (Bt) gene, refuge areas (non-Bt maize), and grasses combined with off-season periods might influence the evolution of resistance in Spodoptera frugiperda (Lepidoptera: Noctuidae), one of the leading agricultural pests targeted by GMOs. We designed the Bt and non-Bt plants in two different arrangements: (a) a seed mixture and (b) strips rows, adding grasses in areas adjacent to the field. We added the seasonal planting dynamics (crop season and off-season), to evaluate a total of six agricultural scenarios. We followed a crop calendar from the United States to create simulations close to agricultural practice. The results showed that the frequency of the resistance allele was strongly related to the landscape arrangements and their dynamics. Since the adult insects are mobile, the seed-mixture scenario increased the frequency of the resistance the most (95.86%), followed by strips (82.10%), without grass fields. The maize harvest made it possible to reduce the frequency of resistance allele below 1%. Based on our results, we can expect that the maintenance of pasture areas, for instance next to the corn crops, will act as a reservoir of susceptible insects during off-season periods.


Introduction
GMO (genetically modified organism) plants expressing genes from Bacillus thuringiensis (Bt) have been used to control insect pests worldwide. However, the resulting selection pressure on target insects may result in evolution of resistance, compromising the benefits of this technology (Lu and Perkins 2013). To overcome this problem, a high-dose strategy has been used, i.e., Bt toxins that kill 95% of heterozygote insects with one resistance gene R. This strategy requires farmers to plant some areas with non-Bt plants to serve as a refuge for susceptible insects (Huang et al. 2011). Crossing individuals from refuge and Bt areas results in heterozygotes that this technology can control, avoiding the rise of new resistant individuals (Carrière and Tabashnik 2001;Sisterson et al. 2005;Tabashnik et al. 2008).
One of the most important agricultural pests is the fall armyworm, Spodoptera frugiperda (Lepidoptera: Noctuidae). Fall armyworms were the first insect pests to develop Communicated by Salvatore Arpaia. resistance to Bt crops in different countries around the world (Molina-Ochoa et al. 2003;Yang et al. 2016;Sisay et al. 2018). These insects have two genetic strains, rice, and corn, which have different host-plant preferences (Piggott et al. 2021;Unbehend et al. 2014;Pashley 1988).
Many studies have used theoretical models to address issues related to resistance to Bt crops and the importance of adopting refuge areas as a strategy for managing S. frugiperda resistance (Vacher et al. 2003;Cerda and Wright 2004;Carroll et al. 2012). These models can potentially provide insights into the evolution of insect resistance in heterogeneous environments composed of Bt and non-Bt plants, considering the fitness of the insect pest on each host type . Including ecological information for different pest, strains can improve these models, which are essential tools in the study of insect resistance and support efficient programs of resistance management, especially in tropical countries (Malaquias et al. 2017). However, information about how resistance evolves in mosaics of host crops, complex landscapes, and buffer zones (e.g., non-cultivated land such as pastures) is still lacking. For instance, pastures can be favorable for fall armyworm survival, sheltering adults, and helping to maintain populations year-round in tropical regions, even when the preferred host is not present (Múrua and Virla 2004;Kanya et al. 2004;Dias et al. 2016).
The role of grass areas in insect fitness was studied by Hay- Roe et al. (2011), who showed that the rice strain is more adapted to Cynodon spp. (Poaceae) than the corn strain. However, the authors did not investigate the role of alternative grass hosts in complex landscapes, such as in a mosaic combining Bt crops with refuge areas. Grass areas can be important in avoiding resistance evolution in the fall armyworm, since they can provide shelter for susceptible insects (Bates et al. 2005;Juárez et al. 2012), especially in the absence of maize fields. In Bt crops with adjacent refuge areas, we assume that susceptible insects have high mortality in Bt with high survival in non-Bt plants. On the other hand, resistant insects can have high survival rates in either Bt or non-Bt areas, regardless of the fitness cost associated with resistance (Caprio 2001).
In order to study the dynamics of S. frugiperda in a heterogeneous landscape composed of host plants (Bt maize, non-Bt maize, and grass), we used an individual-based model appropriate for this purpose (Garcia et al. 2016), based on transition rules and representation of biological traits of the study organisms . We developed a predictive model to describe the dynamics of a rice strain population in an agricultural landscape composed of maize crops surrounded by Cynodon nlemfuensis, a grass commonly found in Florida, USA. This approach was chosen because S. frugiperda is a highly polyphagous species and has previously been recorded from spontaneous and cultivated grasses (Montezano et al. 2018). Using this model, we investigated the role of grass, as an alternative host, in the evolution of insect resistance. We also tested different values of oviposition probability for fall armyworms in grass fields, in order to evaluate whether a S. frugiperda population more highly adapted to this host could affect the frequency of the resistance allele to Bt maize.

Modeling process
A Monte Carlo simulation model was developed to study the impact of grass areas (C. nlemfuensis) adjacent to Bt maize fields. Grass can occur on field edges during or after cultivation, and also when it is part of a crop-livestock system (Ceccon et al. 2015). In this study, we assumed that C. nlemfuensis occurred naturally in the simulated landscape.
We created a bi-dimensional grid with 506 x 506 cells to represent a 25.6-ha agricultural landscape, large enough so that insects did not reach the fixed grid borders during the simulations. The grid was subdivided into four areas containing 250 x 250 cells representing four different patches (arranged as a Cartesian plane), separated by two rows of empty cells (without plants and insects), plus two rows of empty cells at the grid edges. We also used four empty rows and columns as a fixed border during all simulations. We simulated six scenarios that varied in the composition and arrangement of the patches in space and/or time. All scenarios started with Bt fields and refuge zones arranged in one of two configurations: strips or randomly distributed, representing the seed-mixture refuge. Two of the four patches of each scenario were composed exclusively of grass (C. nlemfuensis), simulating intercropping. We added temporal dynamics for two scenarios in which the maize patches were harvested at time step 240, simulating the crop calendar in the USA (USDA 2020), and became empty, representing the off-season, while the grass patches remained available for the insect population. Two scenarios were designed to represent only Bt and refuge zones during all simulations. We did not include the dynamics of these plants in the system. Each time step corresponded to one day, and each simulation had 365 time steps. Below is a detailed description of each scenario and its respective name ( Fig. 1): Constant strip refuge + grass fields (CSG); Constant seed-mixture refuge + grass fields (CSMG); Intermittent strip refuge + grass fields (ISG); Intermittent seed-mixture refuge + grass fields (ISMG); Constant strip refuge (CS); or Constant seed-mixture refuge (CSM).
Each cell could be occupied by immatures, adults, or both insect stages, limited by the carrying capacity (k) of the cell for each stage. One cell could contain up to two immature individuals (k i ) of S. frugiperda, as described by Farias et al. (2001), for each maize plant. For adults, the carrying capacity (k a ) was assumed to be 15 for each cell, that is, the total k of cells was equal to 17. We assumed that both maize and grass crops had the same carrying capacity for both immature and adults, since our study focused on the nutritional role of grasses for the fall armyworm, and did not attempt to treat intra-specific competition on different plants, which is not clearly understood.

Biological parameters of S. frugiperda
We divided the fall armyworm life cycle into two phases: immatures, which included eggs, larvae, and pupae; and adults. We modeled only the dynamics of females since they are responsible for laying eggs and population increase. We also assumed that only they could move through the grid of cells since larvae move only a short distance (Garcia et al. 2019). The biological parameters of the S. frugiperda rice strain, including development time and viability of larval and pupal stages when feeding on maize (non-Bt) and grass, were obtained from experiments carried out in Florida, USA, on grass and by Garcia et al. (2017) on corn. We decided to use an American crop calendar to standardize our study. Fall armyworms were raised under controlled temperature (22 °C), humidity (80%), and photoperiod (14:10 light:dark), and subjected to feeding experiments with C. nlemfuensis. The biological parameters obtained are reported in Table 1. The oviposition period was obtained from the study by Pashley et al. (1995).

Dynamics of immatures
Each cell covered with grass or maize plants could be occupied by additional immatures (i) if an adult lays an egg in it. The number of immature insects in each cell could not reach its carrying capacity (k i ). A cell could be without immature if the insects die during the immature phase, emerge as adults, or if adults did not lay eggs in it. The daily natural mortality probability μ i (d) of each immature i is given by (Eq. 1): where d is the immature developmental time i; ϑ l (s) and ϑ p (s) are larval and pupal viability, respectively, which depend on the crop species (s) on which the immature is feeding; and D l (s) and D p (s) are the durations of the larval and pupal stages, respectively. The mathematical explanation is given in Supplemental Material 1.
The first four days of the immatures' lives corresponded to the embryonic stage. The larva hatched on day five and could feed on the plant on which it was deposited until age 4 + D l (s) + 1 when it pupated. At age 4 + D l (s) + D p (s) + 1, the pupa became an adult if there was enough space in the cell, i.e., with fewer adults than the adult carrying capacity (k a ); otherwise, the pupa was assumed to die and was not included in the simulation.

Adult dynamics
An empty cell could be occupied by an adult (a) if it emerged in this cell or if some adult moved to it. Females were assumed to reach the reproductive period on the third day after emergence. The daily probability of natural adult mortality μ a (d) was set at 0.01 during the fertile period (f s ), which depends on the plant species (s) on which the individual fed as a larva. When a female reached day f s + 1 it was assumed to no longer lay eggs; and was removed from the grid.

Adult movement
Females could move in the grid as soon as they emerged, with no preferred direction, to the surrounding cells, according to the following equation (Garcia et al. 2016): where P is the probability that an adult will travel over a distance S per time step if the number of adults in the cell destination is less than the k a , regardless of the plant species; otherwise, the number of individuals would exceed the assumed carrying capacity. Females could not move to empty cells, except for those that were already in those cells before the harvest period.

Oviposition
After moving, each female randomly chose one of its neighboring cells of radius two (Garcia et al. 2020), and oviposit one egg per time step, with a probability of ovipositing σ(s), which depends on the crop s growing in the chosen cell.
Oviposition was possible only in cells covered with grass or maize and if the cell had not reached its carrying capacity (i < k i ). Otherwise, the female had to search for another cell in its neighborhood of radius two until a cell met these conditions. If all the neighboring cells of radius two had reached their carrying capacity for S. frugiperda immatures (i = k i ), this female did not oviposit in that time step. The probability of ovipositing σ(s) is defined by: where T r is the ratio between the probability of ovipositing on grass and the probability of ovipositing on maize. We tested different values of T to investigate the effect of different acceptability of grass and maize by females for ovipositing, on the dynamics of the resistance allele in S. frugiperda populations. These differences in grass acceptability could

Genetic component
We assumed that resistance to Bt maize was recessive and autosomal (R gene), determined by a single locus and that three different genotypes could occur: susceptible homozygote (SS), susceptible heterozygote (SR), and resistant homozygote (RR) (Huang et al. 2014;Vélez et al. 2014;Garcia et al. 2016). We assumed a 1% initial frequency of the resistance allele R, and the initial genotype frequencies were determined by the Hardy-Weinberg equation: p 2 + 2pq + q 2 = 1, where p is the frequency of allele S and q is the frequency of allele R (p = 1 -q). The components p 2 , 2pq, and q 2 determine the frequency of the genotypes SS, SR, and RR, respectively. For each egg laid, the genotype of the offspring was defined according to the conditional probability functions defined by Garcia et al. (2016) from the Hardy-Weinberg equation (Table 2). Since the paternal genotype was unknown, the offspring genotype was defined based on the maternal genotype and on the frequency of adult alleles p and q within a radius of 35 m from the oviposited cell (Table 2). We consider this distance as the maximum distance reached by adults in one day within the simulated area. Thus, we assumed that the paternal genotype of each offspring was represented by the frequency of the alleles in the adult population within the daily dispersal radius of each female, corresponding to the area where the females would probably have access to the potential fathers of their offspring.
As a fitness cost, we used viability reduced by 20% and a delay of four days in the duration of the larval stage for individuals with at least one copy of the resistance allele (SR or RR) in the absence of selection pressure (Table 3) (Dangal and Huang 2015). We assumed complete resistance, in which the survival and duration of resistant larvae (RR) in Bt maize were the same as those of susceptible homozygotes in non-Bt maize. We assumed that Bt maize produced a high-dose event, which should kill 99.99% of the susceptible larvae, with minimum mortality of 95% for heterozygotes (Environmental Protection Agency 1998). Then, we assumed that all susceptible larvae (SS and SR) that fed on Bt maize were killed (Table 3).

Simulation and data analysis
After the maize harvest, in time step 240, all immatures on Bt or non-Bt crops died, and the adults were able to search for grassy cells and move toward them if these cells were within a radius of 35 m and if there was enough space in these cells (a ij < k a ). Otherwise, the female moved randomly to any harvested cell as described above.
The simulations started with 2,500 sexually mature adult females, randomly released in the 200 × 200 cells in the central grid. The genotype frequency of these initial adults was defined according to the Hardy-Weinberg equation as described above. Fifty repetitions were run for each combination of the scenarios and for the probability of oviposition on grass. The following information was recorded per time step: the total number of insects and the frequency of genotypes SS, SR, and RR in the entire grid.
The mean frequency of resistance and susceptible alleles in the final time step (365) of the 50 repetitions was used to compare the different scenarios, using a χ 2 test. We use the tendency χ 2 test (Liu et al. 2005) to investigate if there was a significant tendency (α ≤ 0.05) toward an increase or decrease in the resistance allele frequency with increasing probability of oviposition on grass for scenarios CGS, CSMG, ISG, and ISMG. All analyses were performed using the R software (R Core Team 2020). Table 2 Conditional probabilities of the occurrence of each offspring genotype, based on the maternal genotype and on the frequency of adult alleles p and q within a radius of 35 m from the oviposited cell, according to Garcia et al. (2016)

Total population dynamics
The relative frequency of the resistance allele differed significantly among the populations, according to the landscape configuration (Fig. 2). The scenarios that maintained the same configuration throughout the simulation (CSG, CSMG, CS, and CSM) showed an increase in the frequency of the resistance allele in the S. frugiperda population, compared to the scenarios that did not. The highest mean frequency was obtained with a seed mixture without grass (CSM; 95.86%), followed by the strip configuration without grass (CS; 82.10%). The presence of grass in the landscape (CSG, CSMG, ISG, and ISMG) contributed to reducing the frequency of the resistance allele compared to the landscapes with only maize (CSM and CS). The intercropping scenarios in which the maize was harvested (ISG and ISMG) were the only landscape scenarios that contributed to reducing the frequency of the resistance allele below the initial frequency (1%), for both refuge configurations. For all oviposition probabilities, the highest frequencies of the resistance allele at the end of the simulation were observed for scenario CSMG, followed by CSG, ISMG, and ISG (Fig. 3). The frequency of the resistance allele increased continuously in the scenarios that maintained maize in the system throughout the simulations (CS and CSM), especially when the alternative host was absent to reduce the selection pressure. The same pattern was observed when grass was included in the system, but the insects showed the lowest probabilities of oviposition on this host (CSG and CSMG). The frequencies of the resistance allele for the different oviposition probabilities in scenarios ISG and ISMG did not differ significantly, except when the oviposition probability on grass was equal to (grass) = 0.99. For the lower probabilities, the effects of the different refuge configurations on the frequency of the resistance allele were dissipated by the absence of selection pressure. After 130 time steps, the frequency of the resistance allele increased substantially until the end of the simulation or until the maize was harvested. The frequency of the resistance allele decreased rapidly after the maize harvest in time step 240 for the scenarios where maize fields were removed from the landscape (ISG and ISMG) and remained at a low level in the subsequent time steps.
There is a direct relationship between the oviposition probability in alternative host and the total number of insects in the scenarios shown in Fig. 4. The crop system simulated in scenario CSG showed the highest insect population, regardless of the probability of oviposition on grass. As expected, the overall total population decreased significantly in the crop systems in which maize was harvested (scenarios ISG and ISMG), which did not differ from each other and were almost half as large as those in landscape scenarios CSG and CSMG.
In Fig. 5 it is apparent that the best strategies to reduce resistance evolution of fall armyworm in maize crops are those closest to the origin of the graph or in the lower-left quadrant, which corresponds to the lowest probabilities of oviposition on grass for scenarios ISG and ISMG.

Relative frequency of genotypes
The higher the relative frequency of resistant (RR) insects (Fig. 6a) through the simulation time, the lower the frequency of homozygous susceptible (SS) insects in the Fig. 2 Relative frequency of the resistance allele (± SE) in Spodoptera frugiperda populations in the different landscape scenarios, in the last simulation time step. Bars with different letters differed significantly (p ≤ 0.05) from each other based on a 2-by-2 χ 2 -test of independence Fig. 3 Relative frequency of the resistance allele in Spodoptera frugiperda populations for the simulated landscape scenarios and different ratios between the probability of oviposition on grass and on maize (T r ) for scenarios CSG, CSMG, ISG, and ISMG, through 365 time steps. Colored areas around the curves correspond to one SE. Arrows indicate the time of maize harvest in scenarios ISG and ISMG, in time-step 240, representing the off-season period. Note different scales of y-axes Fig. 4 Total population of Spodoptera frugiperda in the last time step as a function of the ratio between the probability of oviposition on grass and on maize (T r ) and landscape scenarios CSG, CSMG, ISG, and ISMG (colors). Shaded areas around the curves correspond to one SE population (Fig. 6b). After approximately time step 200, the relative frequency of heterozygous susceptible (SR) individuals increased slightly (Fig. 6c) in scenarios CSG and CSMG, and then returned to the initial frequency level. Comparison of scenarios ISG and ISMG shows that the SR frequency decreased for all T r probability values on grass.
In the landscapes with only Bt and non-Bt maize without grass, the population of resistant insects rose sharply, in contrast to the susceptible homozygous populations (Fig. 7), which rose to a plateau or remained stable over time (scenarios CS and CSM, respectively).

Discussion
Our simulation model showed that different refuge configurations, combined with alternative plant hosts such as C. nlemfuensis, could influence the frequency of the Bt-resistant allele of the fall armyworm in maize. The seed-mixture landscape, in the absence of grass (CSM), increased the frequency of the resistant allele more than when grass was present (CSMG and ISMG). This was expected, since this strategy consists of planting Bt and non-Bt seeds next to each other, allowing faster elimination of susceptible insects (Tabashnik and Carrière 2017;Garcia et al. 2016), especially due to the high mobility of S. frugiperda in the field. The negative effect of the seed-mixture configuration, i.e., its contribution to accelerating the evolution of resistance in the population , was reduced when grass fields were present, because the grass area serves as a refuge for susceptible individuals, thus slowing the rate of increase in the frequency of the resistance allele.
In the fixed landscape scenarios CSG and CSMG, both the total insect population and the relative frequency of the resistance allele were higher than in the two dynamic landscape scenarios ISG and ISMG (which had an off-season period after the maize harvest), even with the higher probability of oviposition on grass. However, the total insect population did not reach the same proportion as observed for ity of oviposition on grass and on maize (T r ). In scenarios ISG and ISMG, arrows indicate the maize harvest in time step 240, beginning the off-season period. Colored areas around the curves correspond to one SE the seed mixture and strip scenarios without grass (CS and CSM). The lower frequency of the resistance allele observed after the maize harvest indicated that a significant part of the pest population with this allele was occupying maize fields, and evidently, when the maize was harvested, the insect population and the frequency of the resistance allele both decreased markedly. However, the total population continued to increase in the grass fields after the harvest, since S. frugiperda is able to reproduce on this alternative host.
Our results also indicate that refuges structured with 20% non-Bt crops were not sufficient to reduce the frequency of the resistance allele under the test conditions and that crop management, such as periods of absence (off-season) must be implemented for optimum results. Furthermore, an ideal agricultural landscape should be heterogeneous in both space and time. This landscape should be composed of Bt and non-Bt plants, which could be maize or any plant species that allows susceptible insects to survive and would delay any increase in the frequency of the resistance allele during the crop season, as well as a period without Bt plants (Cerda and Wright 2004).
In Brazil, forage crops are cultivated for other purposes. Species such as Bermuda grass (Cynodon dactylon) are commonly used as a "green bridge" in the off-season (Ribeiro et al. 2020), or in integrated crop livestock systems where green forage is used for animal feed (Auad et al. 2016). As far as we are aware, published studies have addressed only the contribution of alternative host species to the maintenance of pest populations in the field during the crop season (Ribeiro et al. 2020). The role of these species in maintaining susceptible individuals in the field, as well as their impact on the evolution of resistant insects, remains unclear ).
The population of resistant insects in scenario CS grew nearly as much as in scenario CSM, but the population of susceptible insects was lower in the seed-mixture scenario. This indicates that a strip configuration can maintain susceptible individuals (SR and SS genotypes) in the field for a longer period. Combined with the fitness cost to individuals carrying the resistance allele in refuge zones (in the absence of Bt maize), this configuration can generally delay resistance evolution, even in situations where the initial frequency of the resistance allele is high, as assumed here (Visser and Van den Berg 2020).
In the field, heterozygous individuals are more abundant than RR insects (Dangal and Huang 2015), although our results indicated the opposite (the spatial distribution is given in Supplemental Material 2). Even though they have only one resistance allele, we assumed that heterozygous insects could not survive in Bt fields (high dose) and would have higher mortality on non-Bt plants (maize and grass) due to the fitness cost, compared to homozygous susceptible insects, leading to fewer heterozygotes than the other genotypes in both Bt and non-Bt plants.
To achieve the goals of resistance management, refuge zones must be abundant near Bt fields and contain host plants that are suitable and attractive for target pests (Visser and Van den Berg 2020). Suitability and availability during seasonal changes can increase the potential of a crop to become a source of insects, depending on the adaptability of the insect to this host (Kennedy and Storer 2000). In the case of S. frugiperda, its polyphagy must be exploited from the point of view of IRM, as well as its oviposition preference (Téllez-Rodríguez et al. 2014;Nascimento et al. 2020). There are indications that female S. frugiperda chooses host plants that are less damaged by larvae to deposit their eggs, which would have direct implications for the evolution of resistance. Using computational models, Gustafson et al. (2006) concluded that alternative hosts can act as refuges for Helicoverpa zea (Boddie) (Lepidoptera: Noctuidae) and are important in delaying resistance. They also indicated that pests such as heliothines and Spodoptera spp. can use a wide variety of hosts, and, therefore, modeling studies of insect resistance should include alternative hosts and off-season periods, although the latter is not part of the standard requirements for implantation of Bt crops. Simulations such as these should be conducted for other types of agricultural landscapes and for other pests, in order to provide insights into how resistance evolves in farm fields.