Physiologically based demographics of Bt cotton–pest interactions

Ecological Modelling 191 (2006) 346–359

Physiologically based demographics of Bt cotton–pest interactions I. Pink bollworm resistance, refuge and risk Andrew Paul Gutierrez a,∗ , Sergine Ponsard a,b a

Division of Ecosystem Science, 137 Mulford Hall, University of California, Berkeley, CA 94720, USA b Laboratoire Dynamique de la Biodiversit´ e, University P. Sabatier, Toulouse, France Received 24 February 2004; received in revised form 26 May 2005; accepted 17 June 2005 Available online 23 September 2005

Abstract Transgenic cotton expressing the genes for the production of protoxin of the bacterium Bacillus thuringiensis (Bt) is used to control lepidopterous pests. Among the most successful applications is for control pink bollworm (Pectinophora gossypiella Saunders (i.e. PBW)) in irrigated cotton of the southwestern United States. A major threat to this technology is the development of resistance commonly assumed recessive, autosomal and controlled by a single diallelic gene. A physiologically based, distributed maturation time demographic model of Bt cotton and 10 of its major pests is developed. Here we used the model to examine the population dynamics and resistance development in pink bollworm as modified by weather and spatial and temporal refuges. The dynamics of the other pest species are reviewed in the second paper of this series. The economics of Bt cotton for control of PBW in southern California is put in the context of the historical overuse of pesticides and the alternative short season cotton technology. The analysis posits that in the short run, the Bt cotton may be risk reducing and economic, but in the longer term it may be risk increasing. © 2005 Elsevier B.V. All rights reserved. Keywords: Demographic model; Bt cotton; Pink bollworm; Resistance; Refuges; Risk

1. Introduction Susceptibility to pesticides in target agricultural pest populations is a valuable natural resource in agriculture

DOI of related article:10.1016/j.ecolmodel.2005.06.002. Corresponding author. Tel.: +1 510 642 9186; fax: +1 510 643 5098. E-mail address: [email protected] (A.P. Gutierrez). ∗

0304-3800/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2005.06.001

that must be managed carefully to keep the efficacy of a pesticide as long as possible, especially if the pesticide is environmentally friendly and if replacement options are more damaging. The incorporation of genes for the production of Bt toxin from the bacterium Bacillus thuringiensis by transgenic crops is thought to be an environmentally friendly innovation (Luttrell and Herzog, 1994). Cottons expressing the Cry1Ac ␦-endotoxin (i.e. Bt toxin) have given good control of pink bollworm (Pectinophora gossypiella

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

Saunders) and tobacco budworm (Heliothis viresens (F.)) (Williams, 2000), but variable control of other pest species is occurring (Luttrell et al., 1999). In 1999, about 8 million ha in the United States were planted to Bt cotton (USDA, 2000), and pressure is growing to introduce it worldwide (USDA-FAS, 2000; Pray et al., 2000; Qaim and Zilberman, 2003). To examine the impact of Bt toxin on population growth and resistance in pests of cotton requires the development of holistic models (Parker et al., 2000; Stark and Banks, 2003). The models presented here and in a related paper (Gutierrez et al., 2005) realistically captures the field biology of Bt cotton and of its 10 herbivore pests [pink bollworm (PBW), tobacco budworm (TBW), bollworm (BW, H. zea (Broddie)), fall armyworm (FAW, Spodoptera fugiperda (J.E. Smith)), beet armyworm (BAW, Spodoptera exigua (H¨ubner)), cabbage looper (CL, Trichoplusia ni (H¨ubner)), soybean looper (SBL, Psuedoplusia includens (Walker)), plant bugs (LY, Lygus spp., Heteroptera), whitefly (WF, Homoptera), cotton bollweevil (Anthonomus grandis Boh.)].1 Here we model the population dynamics and possible development of resistance in pink bollworm to Bt toxin, while the biology and dynamics of the other pest species in Bt cotton and the references for both papers are found in the second paper of this series (Gutierrez et al., 2005). There is a well-known trade-off between simple models amenable to mathematical analysis and more complex and realistic models that seek to capture the subtleties of the biological interactions. Mills and Getz (1996) reviewed population dynamics models and emphasized the importance of bottom-up effects, especially of plant effects on higher trophic levels, effects that are critical in the analyses of transgenic crops used in pest control. Despite knowing little about movement and migration of cotton pests (Fitt et al., 1995; Peck et al., 1999; Tang et al., 2001), simple theoretical models have been developed to assess strategies for delaying resistance to Bt toxin (e.g. Mallet and Porter, 1992; Tabashnik, 1994; Caprio, 1998; Onstad 1 The cotton bollweevil has historically been a primary pest of cotton in the southern United States (and elsewhere in the Americas), but it is currently controlled by pheromone-based technologies. It is not susceptible to the Bt toxin and hence was not included in the analysis. Whiteflies are also thought to be immune to the toxin and are not reviewed.

347

and Gould, 1998; Vacher et al., 2003, 2004; Linacre and Thompson, 2004). Although simplified and omitted plant dynamics, the models have heuristic value giving general insights into the effects of refuge size and spatial heterogeneity on Bt resistance management. Ru et al. (2002) developed an age-structure simulation model for bollworm (Helicoverpa armigera (Hubner)) that included the effect on pest dynamics and resistance of declining Bt toxin levels with plant age but did not include plant dynamics. The analysis questions the utility of the Bt refuge strategy in cotton in China. Our model of Bt cotton and its pests differs from prior models in several ways: it is a physiologically based; it has age (mass) structured distributed maturation times demographics; the same models for resource acquisition and allocation are used across species, be they plant or animal (Gutierrez, 1996); species-specific supply/demand ratios govern much of the temporal dynamics; pest dynamics modules are integrated in the system and may be implemented singly or in combination using Boolean variables; the biology of Bt cotton includes the decline in toxin levels with plant and subunit age; the lethal and sub-lethal effects of Bt toxin concentrations on pests and their vital rates are included; there is feed back between pest attack and plant compensation; the model includes pest time varying preferences for plant sub units; insecticides of varying toxicity may be included; observed weather is used to drive the dynamics of the system (see Appendix A). 1.1. Review of transgenic Bt cotton B. thuringiensis (Bt) based insecticides are examples of low-impact pesticides that have been used since the 1950s without significant development of resistance in target pests (Federici, 1999). These microbial Bt preparations are highly sensitive to light and heat and decay rapidly after being applied in the field, while Bt protoxin in transgenic cotton is continuously produced and it is quite stable. Susceptible pests feeding on Bt cotton tissues activate the protoxin in their gut killing susceptible individuals (Ostlie et al., 1997). Constant presence of Bt protoxin (hereafter toxin) in transgenic plants and the planting of Bt crops on a broad scale is thought to make the development of resistance more likely (Gould, 1997, 1998; Mellon and Rissler,

348

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

1998; Hilder and Boulter, 1999; Tabashnik et al., 2000). Although resistance has been demonstrated in several cotton pests in the laboratory (Tabashnik et al., 2000), resistance has not developed in the field (Tabashnik et al., 2003, cited in Gonzalez-Cabrera, 2003). Selection for resistance to Bt toxin is akin to the development of resistance to chemical insecticides but the process is subtle. Most models of resistance to Bt toxin assume it is recessive, autosomal and controlled by a single diallelic gene with Mendelian inheritance (Tabashnik, 1994; Gould, 1998, Tabashnik et al., 2000; Liu et al., 2001). The Bt toxin is assumed to kills susceptible heterozygous (Rr) and homozygous (rr) individuals causing rapid selection for the few homozygous resistant (RR) survivors. To combat the development of resistance, the so-called “high dose-refuge” strategy was been mandated (Gould, 1986, 1998; EPA, 2000). The tenets of the spatial refuge strategy assume that the initial frequency of the resistance allele is low, the pest is susceptible, and a proportion of the cultivated area is planted with conventional non-Bt cultivars to act as a “refuge” for the susceptibility allele. In this strategy, Bt-free refuges are expected to be sources of susceptible adults to mate with resistant ones to produce mostly susceptible heterozygous progeny. Spatial refuges for polyphagous species also occur in non-Bt host (native vegetation and other crops). Temporal refuges for pests may occur due to the time varying concentration of the toxin in the plant and differences in susceptibility of pest species (Gutierrez et al., 2005). Initially, the frequency of the resistance allele in PBW was thought to be two orders of magnitude higher than the <0.001 assumed in resistance management models (Patin et al., 1999; Tabashnik et al., 2000), but this estimate has since been discounted (R. Roush, pers. comm.). The apparent lack of resistance in field PBW populations is explored here, and cultivation of Bt cotton in the desert valleys of southern California is put in the context of the historical overuse of pesticides and evaluated against short season cotton for PBW control. 1.2. A model for Bt cotton–pest interactions Our demographic models for cotton growth and development and pest population dynamics are based on the widely used distributed maturation time model (Vansickle, 1977; see Di Cola et al., 1999 for related forms). Several models for cotton growth and develop-

ment have been developed including object-oriented models (Sequeira et al., 1991; Lemmon and Chuk, 1997). These models lack demographic properties and were not linked to pest and higher trophic level models. The cotton–pest bollworm system used here was developed, parameterized and validated in earlier work (Gutierrez et al., 1975, 1977a, 1984, 1991a,b; Stone and Gutierrez, 1986a,b; Stone et al., 1986 and others). The models incorporate weather driven meta-physiology to drive the dynamics (see Gutierrez, 1996) that has its roots in the early physiological plant modeling work of C.T. de Wit in the Netherlands (e.g. de Wit and Goudriaan, 1978). For clarity, only critical biological features of the Bt cotton–PBW model are described, while the mathematics are summarized in Appendix A and in the papers cited therein. 1.2.1. Cotton The Bt cotton–PBW dynamics model consists of eight linked “functional populations”: age structured mass population dynamics models for leaf, stem and root tissues, age-structured fruit mass and fruit numbers (Gutierrez et al., 1975, 1991a; Wang et al., 1977). The plant subunit models are linked via the metabolic pool model (Eq. (1)) that predicts photosynthetic rates under various plant states, soil moisture and nitrogen levels and weather and its allocation. GR(t) = (φ∗ (t)(D(t)β − Q10 )λ

(1)

Photosynthate for reproductive and vegetative growth (GR) is computed as the maximum demands (D) scaled by the product of all supply–demand ratios of essential resources (φ* ) such as light, soil moisture and nitrogen, and then allocated in priority order to wastage (1 − β), respiration (i.e. Q10 ), costs of conversion (λ) (see Appendix A). Photosynthate production affects subunit growth as well as birth and death rates. The plant model was modified to include the time varying expression of Bt toxin and its effects on pest dynamics and resistance. The concentration of Bt toxin varies among cotton cultivars and with plant and plant sub-units age (Finnegan et al., 1998; Fitt, 1998; Greenplate, 1999; Adamczyk et al., 2000, 2001a), but the factors controlling the expression of Bt toxin in the plant are poorly understood (Sachs et al., 1998). The highest concentrations are found in young buds and flowers and the lowest in older buds and

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

349

Fig. 1. Sub-models for pink bollworm biology: (a) percent diapause induction as a function of temperature and minutes of day length (see Gutierrez et al., 1981), (b) preference of PBW for different age fruit, (c) PBW larval developmental time in degree-days (dd > 12 ◦ C) on different age conventional cotton fruit and (d) the survival rate per day of PBW larvae on different age Bt cotton fruit.

bolls (Adamczyk et al., 2000, 2001a). Only the age dependent pattern in plant subunits was captured (see Appendix A). 1.2.2. Pink bollworm PBW population numbers are modeled by age and genotype as well as by age (toxin level) of infested plant subunit attacked. The number and mass attributes are included in the models of other pests. PBW larvae and the fruit they infest age over time, and hence the larvae experience different levels of toxin that affects their developmental time and mortality, as well as the fecundity of the adult stage. To capture this biology, a two-dimensional distributed maturation time dynamics model with ages of PBW and host fruit as dimensions was used (Stone et al., 1986). 1.2.3. Phenology Five to six PBW generations occur in Arizona and southern California with its geographic distribution limited to the relatively frost free areas (Gutierrez et al., submitted for publication). In fall, mature PBW larvae on conventional cotton enter diapause in response to decreasing temperature and photoperiod (Fig. 1a, But-

ler et al., 1978; Gutierrez et al., 1981)2 , of which less than five percent survive the winter (Bariola 1984). Survival of diapause larvae from Bt cotton is considerably lower (Carri`ere et al., 2001c). PBW adults emerge from diapause from early spring to mid July in a pattern captured by a Gompertz function (Gutierrez et al., 1977a), and are apportioned to the first adult age class of the appropriate genotype population using the frequency of resistance of diapause larvae that produced them. Modeling this biology allows linkage and simulation of consecutive season. 1.2.4. Reproduction PBW females lay eggs singly on or near buds and bolls avoiding very young buds, flowers and old bolls (Fig. 1b, Lukefahr and Griffin, 1957; Westphal et al., 1979; Henneyberry and Clayton, 1982; Stone and Gutierrez, 1986a). Each maturing boll may host as many as 15 larvae. On conventional cotton, fecundity of adults reared on buds or emerging from diapause is 2 The model for diapause initiation and termination has been used in the field in the Brazil, Egypt, The Sudan and USA (e.g. Gutierrez et al., 1977a, 1986; Gutierrez, 1996; Russell, 1995).

350

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

78% that of larvae developing on bolls (Westphal et al., 1979). At maximum toxin concentrations, fecundity of surviving adults is 50% that of adults reared on conventional cotton (Carri`ere et al., 2001b), but as toxin concentrations decline during the season, adult fecundity is assumed to increase proportionally. Assuming random mating, the eggs produced by all surviving adult genotypes are allocated in the correct number to the first age class of each genotype population. Hatching larvae attack fruit age classes at frequencies weighted for age class abundance and adult oviposition preferences (Westphal et al., 1979). 1.2.5. Developmental time Larval developmental times are longest in young buds, flowers and old bolls (Fig. 1c, cf. Lukefahr and Griffin, 1962). The shortest larval development time is 240 degree-days (dd) and the longest may be 10× that value on very young or old fruit. PBW larvae developing on Bt cotton have been shown to require 20% longer to complete development (Liu et al., 1999, 2001; Peck et al., 1999). In the model, the developmental time at maximum Bt toxin concentrations is 20% times longer in resistant genotype and 30% longer in heterozygous and homozygous recessive genotypes. 1.2.6. PBW mortality On conventional cotton, intrinsic survivorship is 0.55 (Stone et al., 1986). In addition, approximately 50–60% of the buds initiated are shed causing significant mortality to early season PBW eggs and larvae. Infested and healthy buds are shed at the same rate because the larvae feed on non-vital structures (Westphal et al., 1979). Resistant larvae are assumed to survive in Bt cotton, while survival of susceptible individuals is a function of the toxin concentration in attacked fruit. The pattern of toxin concentration in fruit and the effect on susceptible larval survivorship per day (d−1 ) (Fig. 1d) is captured as the normalized inverse of the function for PBW developmental time on fruit age (Fig. 1c). For example, a cohort of susceptible PBW larvae in young Bt squares (buds) experiences a 0.85 survivorship rate d−1 , and the average developmental time is 25 days, then the average larval stage survivorship is 0.8525 days = 0.017.

Predation in conventional cotton is an important source of mortality to PBW eggs (i.e. 0.3 d−1 ; cf. A.P. Gutierrez unpublished data; Sansone and Smith, 2001) but not larvae that penetrate fruits where they are largely invulnerable. Predation in Bt cotton is ca. 28% lower due to reduced longevity of predators feeding on Bt intoxicated prey (Ponsard et al., 2002). The combined survival from egg to pupa from all sources is 0.55 × 0.783 days × 0.8525 days = 0.0045. If we include egg–larval mortality that occurs in shed buds, the value is roughly halved. On mid-age Bt bolls, larval development is about 19 d and with a daily survivorship of about 0.94 d−1 making egg–larval survival 8%. We note that despite this high survivorship, the population on mid-age bolls is expected to be small due to the very high mortality larval populations experienced feeding on Bt buds earlier in the season. 1.2.7. Bt resistance Bt resistance in PBW has not been found in the field, but Tabashnik et al. (2000) reported 0.7–1% emergence of adults from field-collected larvae reared on late season Bt cotton bolls having low levels of toxin. In the absence of high levels of resistance, this suggests the presence, albeit very small, of a temporal refuge for susceptibility in the crop. The development of pest resistance to Bt toxin is modeled using simple Hardy–Weinberg genetics. In a single Bt toxin cultivar, three PBW genotypes arise requiring a dynamics model for each, but in a two-toxin cultivar nine genotype models are required. 1.3. Initial conditions Single and multi-year simulations were run assuming irrigated long season Bt cotton planted at 6.6 plants m−2 and non-limiting nitrogen. Meteorological data (daily max–min temperatures, solar radiation, rainfall, daily runs of wind and relative humidity) recorded at Brawley, Imperial County, California during the period 1 January 1983 to 31 December 2002 were used to run the model (UC/IPM, 2001), but weather data from other cotton growing areas of the world can be used as well. Season length was January 1 when adults could begin emerging to October 30 when crop destruction occurred.

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

The initial PBW population density was 2 diapause larvae plant−1 . In multi-year runs, predicted diapause larvae at the end of the prior season corrected for winter mortality provided the link to the next season. A high initial frequency of resistance was assumed (R = 0.15) for heuristic purposes to illustrate the direction of change of resistance. The effects of spatial refuges were included in the model as a daily dilution rate of the resistance allele by incoming susceptible migrants. Temporal refuges effects arise as part of the local population genetics and dynamics.

2. Simulation results 2.1. PBW dynamics on conventional cotton In a single season run, the model captures the first fruiting cycle and the partial second late season cycle common in long season cotton (Fig. 2a).

351

It also captures the prolonged pattern of adult suicidal emergence in spring with only a small fraction of adults emerging during the bud-boll period (Fig. 2b) (e.g. Rice and Reynolds, 1971; Gutierrez et al., 1977a; Stone and Gutierrez, 1986a). Five and a half larval generations were produced with the final one occurring after mid September and yielding mostly diapause larvae (Fig. 2c) on the partial second fruiting cycle (Fig. 2a). This illustrates that the production of diapause larvae could be avoided using short season cotton and early plowing of the crop residue (see Burrows et al., 1982; Stone and Gutierrez, 1986a,b). The dynamics of 17 years of conventional cotton and PBW are shown in Fig. 3. The fruiting patterns differ from year to year with some years producing a larger second fruiting cycle (Fig. 3a). The dynamics of PBW larvae are greatly influenced by the density of the previous season’s diapause larvae and by the bottom-up effects of within season fruit dynamics (Figs. 2a and 3a).

Fig. 2. Simulated cotton–pink bollworm dynamics on non-Bt cotton: (a) the pattern of buds (squares), immature and mature bolls, (b) the proportion (×10) of adults emerging during spring and the number of larvae entering diapause in late summer and (c) the dynamics of eggs, larvae and adults. Brawley, CA, weather data (1983) and 6.6 cotton plants m−2 were used in the simulation.

352

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

Fig. 3. Simulation of 17 years of non-Bt cotton–PBW interactions: (a) dynamics of immature cotton bolls and (b) the dynamics of PBW larvae. Brawley, CA, weather data (1983–2000) and 6.6 cotton plants m−2 were used in the simulation.

2.2. PBW dynamics on Bt cotton 2.2.1. No spatial refuge A single season simulation on Bt cotton with an initial frequency of the resistance allele R = 0.15 and no spatial refuges for susceptibility is shown in Fig. 4. PBW infestation does not alter cotton fruiting patterns (Westphal et al., 1979, Fig. 2a); hence crop dynamics are the same as in Fig. 2a and are not illustrated. Early in the season, fruit age-structure is biased toward buds with high levels of toxin, while late in the season the bias is towards bolls with low toxin levels (Fig. 2a) (Adamczyk et al., 2000, 2001a). This creates strong early season selection pressure that wanes late in the season. High Bt selection pressure returns during the late season partial fruiting cycle. The combination of prolonged PBW emergence from diapause and the decline of Bt concentration in fruit allow a small number of susceptible larvae to survive (2.7 m−1 ) during late summer (Fig. 4a) and some to enter diapause (Fig. 4b). Of these, 6.7% survive winter mortality but most experience suicidal emergence the next season, and fecundity of survivors is 50% less. In multi-year simulations of Bt cotton with no refuge, PBW summer and fall diapause larval numbers

Fig. 4. Simulated cotton–pink bollworm dynamics on Bt cotton: (a) the dynamics of eggs, larvae and adults and (b) the proportion of adults emerging during spring and the number of larvae entering diapause in late summer. Brawley, CA, weather data (1983–2000), R = 0.15 and 6.6 cotton plants m−2 were used in the simulation.

decrease to near zero after 4 years, while the frequency of resistance (R) increases to near unity after 6 years allowing PBW populations thereafter to increase to 70% of the level seen in simulations of conventional cotton (Fig. 5a versus Fig. 3b). Extrapolating back from the initial value of R = 0.15, the time to R = 0.05 is about 4–5 years, suggesting that in the absence of spatial refuges, resistance in PBW from a 0.05 level would develop in about 10 years after the introduction of Bt cotton, but 12–13 years would be required for outbreaks of PBW to begin. The lag in PBW outbreaks after resistance develops is due to low PBW density and the sub-lethal effects of the toxin on increased developmental time and reduced fecundity. In the absence of spatial refuges, increasing the toxicity of the Bt toxin decreases the time for full resistance development. For example, increasing Bt mortality 10% decreases the time to resistance breakdown to 3 years, and using higher toxicity values decrease the time further (not shown). 2.2.2. With a spatial refuge Assuming a 5% spatial refuge creates a 5% dilution d−1 of the resistance allele, outbreaks of PBW are delayed a further 9–10 years and resistance levels settle

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

353

Fig. 5. Simulation of 17 years of PBW larval dynamics and the frequency of resistance to Bt toxin: (a) no refuges, (b) 5% refuge, (c) 10% refuge and (d) a stochastic 5% ± 2.5% refuge. Brawley, CA, weather data (1983–2000), R = 0.15 and 6.6 cotton plants m−2 were used in the simulation.

near 0.95 compared to the no-refuge scenario (Fig. 5a versus Fig. 5b). This is the “high dose-refuge” strategy discussed earlier. Increasing the refuge to 10% stops the development of resistance altogether causing it to decline to near zero (Fig. 5c). Assuming an average dilution rate of 5% with a random ± 2.5% variation, resistance rises to about 0.9 and larval populations remain low initially, but increasing somewhat in later years (Fig. 5d). Assuming a 5% spatial refuge and that selection for normal developmental time and fecundity occurs with increasing resistance, the time to PBW outbreaks decreases to 3 years (not shown). We note there is no evidence that this adaptation is occurring in the field.

3. Discussion The factors affecting the development of resistance in pests to transgenic crops expressing the protoxin of the bacterium B. thuringiensis are complex and require methods of analysis that incorporate many interacting factors. The cautions in Pielke and Conant (2003) on prediction in complex systems have been taken to heart in the formulation and use of our model—especially the admonitions on complexity by N. Oreskes cited therein. Bt cotton currently provides excellent control of pink bollworm (PBW), but the development of resistance is a potential but unrealized problem (e.g. Tabash-

nik et al., 2000). Spatial refuges consisting of conventional cotton to delay or stop the development of resistance have been mandated (EPA, 2000) despite the fact that the underpinning theory for the refuges was based on simple models in the absence of sound information on migration rates from spatial refuges and other ecological factors (Gould, 1986, 1998). Vacher et al. (2003) concluded that the currently mandated spatial refuge strategy was suboptimal, and for H. virescens the optimal refuge should be ca. 25%, much higher than the mandated 5% refuge (EPA, 2000). We conclude that the failure of resistance development in Arizona after 8 years of Bt cotton is due to the massive spatial refuge provided by the roughly 50% non-adoption of Bt cotton (see Tabashnik et al., 2000, p. 12983; Carri`ere et al., 2001a). Under such conditions, only an exceedingly low level of resistance could develop, possibly below detectable levels. However, as the adoption of Bt cotton increases, resistance could become a problem despite the mandated 5% spatial refuge, and especially if fitness increases in resistant individuals. Evidence for temporal refuges for PBW is sparse. Tabashnik et al. (2000) found more than 99% mortality of PBW larvae in the field, but the same data suggest that susceptible individuals may survive in Bt bolls. For example, 19.3% of F1 progeny of individuals collected in the Eloy-Bt field survived exposure to 3.2 ␮g Cry1AC and 11.8% survived at 10 ␮g. More impor-

354

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

tant, the same pattern was found for F1 progeny from three non-Bt fields (i.e. locations Safford, Solomon and Stanfield). We posit that a temporal refuge (albeit very small) is produced by the extended emergence of PBW adults from diapause during mid summer when bolls with low Bt toxin concentrations predominate (see Greenplate, 1999). Susceptible individuals surviving in this small temporal refuge would contribute to delaying resistance, but its importance would be secondary to the effects of the larger spatial refuge created by the non-adoption of Bt cotton or even a mandated 5% refuge. We also note that the high field estimates of larval mortality observed by Tabashnik et al. (2000) could have been confounded by the intrinsic poor suitability of late season bolls (see Fig. 1c for developmental times). Our model suggests that if resistance were to develop, resurgence of PBW populations would be slowed by the sub lethal effects of the Bt toxin on developmental times, fecundity and over-winter survival. Temporal and/or spatial refuges (of various kinds) play a more important role in delaying resistance buildup in polyphagous noctuid pests such as bollworms and defoliators that are inherently more tolerant to Bt toxin than PBW, are stronger fliers, and may have alternate non Bt hosts that provide spatial and temporal refuges (Adamczyk et al., 2000a; see Gutierrez et al., 2005). Liu et al. (1999) and Peck et al. (1999) argued that increases in PBW developmental time in Bt cotton would favor homogamy and therefore accelerate resistance, but overlapping generations and distributed maturation times of cohorts would tend to decrease this effect. Cousteau et al. (2000) question the oft-made assumption of fitness cost to individuals resistant to pesticides (e.g. Gutierrez et al., 1979). Large fitness costs occur on Bt cotton as PBW larval developmental times are longer, the pupae are smaller resulting in lower adult fecundity, and the survivorship of diapause larvae is also lower (Ashfaq et al., 2000; Carri`ere et al., 2001b,c). Last, if the expression of Bt toxin in fruit were high and constant through the season as some as advocates of the technology propose, development of resistance would be faster contrary to the usual assumption (Fitt, 1998). With wider acceptance of the Bt technology (e.g. China and India) and especially if refuges are reduced or ignored, rapid buildup of resistance in PBW is more likely to occur.

Fig. 6. Summary of data on Imperial County, CA, cotton: (a) acreage and (b) yield (Imperial County, University of California Extension, County Agricultural Commissioner Office Reports, 1976–2000).

3.1. The need for Bt cotton for PBW control PBW was the major cotton pest in the desert regions of Arizona and southern California when planting longseason cotton was practiced. In the early 1980s, the cotton industry in Imperial County collapsed to less than 10% of it maximum (Fig. 6a, Imperial County Commissioner Annual Reports, 1986–2002) due to economic structural problems, insecticide costs for control of PBW and induced secondary pests, and declining cotton yield and quality. Most of the insecticide had been used to protect the 1/4 to 1/2 bale of cotton that often accrued from the second late season fruiting cycle (Figs. 2a and 3a, see Burrows et al., 1982; Stone et al., 1986a).3 The insecticides used against PBW caused massive outbreaks of secondary pests (primarily defoliators, bollworms and budworms) that became more damaging than PBW (cf. van den Bosch, 1978; Burrows et al., 1982). The problem became especially acute when insecticide resistance developed in budworm and bollworm (Gianessi and Carpenter, 1999), but surprisingly not in PBW. Several technologies were tried to make long season technology viable, among

3 Additional mitigating factors in the decline of the industry were changes in government subsidies and the insecticide induced outbreaks of whiteflies (Bemisia tabaci (Gennadius)).

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

them combinations of insecticides, pheromone for mating disruption (Stone et al., 1986), and massive releases of sterile males.4 The part PBW played in the decline was largely avoidable as early termination (mid-late September) and plow-down of the crop could reduce infestations below economic levels (see Fig. 2, Burrows et al., 1982; Stone and Gutierrez, 1986b). This technology had been available since the mid 1970s, and was made mandatory in Imperial County in 1989 after chemical control failed (Chu et al., 1996). Yields using short season cotton reached pre-crisis levels (Fig. 6b) and quality doubled, all largely without the use of insecticide (Chu et al., 1996; Eric Natwick, UC Extension, pers. comm.). Long season Bt cotton with a 5% refuge was introduced to Arizona in 1996 with good control of PBW and other lepidopterous pests (Wilson et al., 1992). USDA data for the first 5 years of Bt cotton show that insecticide use was reduced nationally an average of 15% (Benbrook, 2000, 2001), and in Arizona the reduction was much greater (Carri`ere et al., 2001a). These positive results opened up an apparent safe return to long season cotton in southern California in 1999 with resistance seen as a future, theoretical problem. However, variable control of bollworm, beet armyworm and other noctuid pests is occurring (Luttrell et al., 1999) requiring insecticides use for their control (Mahaffey et al., 1995; Lambert et al., 1996; Turnipseed and Sullivan, 1999; USDA, 2000; EPA, 2000; Gianessi and Carpenter, 1999; Benbrook, 2000; Ru et al., 2002). In California, the field evidence suggests the Bt technology is not needed for PBW control in short season cotton in southern California (Fig. 6), nor in the great Central Valley of California where PBW does not thrive (Gutierrez et al., 1977a, submitted for publication) and where cotton can be grown largely pesticide free (Falcon et al., 1968; Gutierrez et al., 1975; van den Bosch, 1978). Proponents of Bt cotton for PBW control argue persuasively that it is preferable to highly toxic pesticides, and while there is no argument on this point, there are 4 The USDA APHIS, in cooperation with California cotton growers and the California Department of Food and Agriculture (CDFA) continue to release sterile PBW males in the Great Central Valley of California to prevent the pest from establishing (USDA Bulletin, 1995), this despite the fact that the northern limits worldwide of PBW is restricted by frosts (Gutierrez et al., 1977a, in press) making it unlikely for PBW to over winters there.

355

other issues that need to be considered. First, higher yields have yet to materialize in southern California (Fig. 6b) and embedded in the technology are the seeds of its demise—i.e. resistance. Claims of 70% yield increases in India by economists Qaim and Zilberman (2003) using industry data are suspect and what gains did accrue must be viewed in the context of reduced insecticide pest disruption. Had Bt cotton been introduced into California during the 1980s at the height of the insecticide crisis, it would have also been falsely viewed as a success story for the technology (Gutierrez, 2004). Furthermore, a strategy for delaying resistance in one pest may not be suitable for delaying it in other co-occurring pests (Gianessi and Carpenter, 1999), and the Bt technology may have unexpected collateral effects including adverse effects on natural enemies that acquire the Bt toxin feeding on Bt intoxicated prey (Hilbeck et al., 1998; Ponsard et al., 2002) and increase levels of tolerant pests (see Gutierrez et al., 2005). Last, contrary to what some proponents of unfettered biotechnology assert, the market does not always chose the most suitable technology—economic welfare theory and practice has consistently shown otherwise (e.g. Regev, 1984; Zadoks and Waibel, 2000).

Appendix A. A distributed maturation time demographic model for Bt cotton and cotton pests In physiologically based plant–herbivore system models, the plant integrates both bottom-up effects of weather and edaphic factors, and the top-down effect of herbivory and intra-specific competition (Gutierrez, 1996). We illustrate the modeling approach using cotton and pink bollworm, noting that the same approach applies, with some modifications, to other pest species in the system. The cotton model used here may be viewed as a canopy model consisting of five linked demographic models for ns{n = 1, ns} linked functional Bt cotton–PBW populations: mass of leaves {n = 1}, stem {2} and root {3}, and for fruit mass and numbers {4, 5} (cf. Gutierrez et al., 1975, 1993). The age structured number dynamics of pink bollworm requires three demographic models; one for each resistance genotype {6–8} dynamically linked to the cotton fruit models {5,

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

356

6}. A time-invariant distributed-maturation time agestructure model (Eq. (A1), Vansickle, 1977) is used to model all of the populations. von Ark et al. (1984) used a similarly structured models in their study of pesticide induced outbreaks of whitefly on cotton. Legaspi et al. (1996) used the distributed maturation time model to examine the biological control of cotton boll weevil by the parasitoid Catolaccus grandis, but their model did not include plant dynamics. For cotton, the dynamics of the ith age class (Ni (t), i = 1, . . ., k) of the nth functional populations (subscript n ignored) are modeled using (A1). dNi (t) = ri−1 (t) − ri (t) − µi (t)Ni (t) dt

(A1)

Aging occurs via flow rates ri−1 (t) from Ni−1 to Ni , births enter the first age class of the population, deaths at maximum age exit the last or kth age class, and net age-specific proportional mortality (losses and gains) from all factors is included in −∞ < µi (t) < +∞. The mean developmental time of a population is v with variance V with the age width of an age class being v/k and k = v2 /V . The number of individuals (or mass units) in age class i is Ni (t)
= ri (t)v/k, and that
k in the total population is N(t) = ki=1 Ni (t) = v(t) i=1 ri (t). If k k is small, the variability of developmental times is large and vice-versa. A value of k = 45 was chosen to produce a roughly normal distribution of developmental times. The initial condition for all cotton subunit populations must be specified to assure the uniqueness of the solution. The numerical solution for (A1) is (A2). dri (t) dt k = v(t)



   v(t) dv(t) ri−1 (t) − 1 + µi (t) + ri (t) k kdt (A2)

The developmental time of PBW larvae varies with fruit host age and toxin concentration, and both the host and the larvae age on their own time scale. Hence larvae initially infesting specific age fruits at time t will in the course of their development experience changing host characteristics that affect their developmental times, mortality and potential fecundity as an adult. To handle this biology, a two-dimensional time-invariant distributed maturation time model with flows in the

fruit age and age of pest dimensions is utilized (Eq. (A3), Stone and Gutierrez, 1986a). dNi,j (t) = ri−1,j−1 (t) − ri,j (t) − µi,j (t)Ni,j (t) dt

(A3)

The mean developmental rate of a cohort of larvae (v(t, i, j)) is transient and depends on pest genotype, and the age and Bt toxin concentration of host fruit. Hence, if i is larval age and j is fruit age, the model is updated for flow first in the ith and then the jth dimension taking care to correct for differences in developmental time scales between cells. For convenience, the net proportional mortality term µnij (t)Ni,j is applied in the ith dimension and assumed zero in the jth dimension. This scheme allows mortality due to subunit shedding in some species (e.g. PBW) and to Bt toxin levels to be applied to larvae in each i,j cohort. PBW population density across genotypes (n) is: N(t) =

3  k  k  n j=1 i=1

Nnij (t) =

3  k k  v(t, j)  n j=1

k

rnij (t).

i=1

Bt toxin concentrations in food consumed by noctuid pests also varies through the season and affects their developmental time, fecundity and survival. However, because it is not possible to associate a specific noctuid cohort with a specific cohort of host, a time varying distributed maturation time model is used (Gutierrez et al., 2005). A.1. Physiological time and age Cotton and its pests are poikilotherms and hence time and age in the model are in physiological time units. The linear degree-day model (A4) or non-linear models (e.g. Lactin et al. 1995) may be used to model the temperature (T) dependent development rate (v(T )). Here we use the linear model because sufficient data for all species was not available. v(t(T )) =

1 = c1 + c2 T (t) v(t(T ))

(A4)

The lower developmental threshold for cotton and Lygus is 12 ◦ C, 10 ◦ C for boll weevil and pink bollworm, 12.2 ◦ C for all noctuid pests, and 13.5 ◦ C for whitefly. Time stepping in the model is a day of varying physiological time as appropriate for each species.

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

Development time of pest larvae (vn (aj, [Bt])) depends on larval genotype (n = RR, Rr, rr), fruit age (aj ) (e.g. Stone and Gutierrez, 1986a), and toxin concentration [Bt]. Specifically, for PBW, vn (af , [Bt]) = φn(Bt) v(af )φ(aj )

(A5)

v(aj ) is the minimum developmental time on conventional cotton. af is the age of the host fruit.      1 if [Bt] = 0 

φn(Bt) =

 

1 + 0.2[Bt]/[Bt]max

1 + 0.3[Bt]/[Bt]max

the details of α.

   −αR . S = Dh(u) = D 1 − exp D

h(u) is the proportion of the resource demanded (D) obtained. Note that intra-specific competition enters the model via  the ratio of available resource to population demand −αR in (A7). D

For cotton, R is the light energy (cal m−2 d−1 ) per unit of ground at time t multiplied by a constant that converts it to g dry matter m−2 d−1 . For a pest species, R may be the sum of all age fruit (or leaves, age = j) corrected for preference (0 ≤ ξ j ≤ 1). R=

J 

ξj R j

(A8)

j=1

A.2. Growth rates Per capita resource acquisition is allocated in priority order to respiration (i.e. Q10 ), and costs of conversion (λ) to reproductive and vegetative growth rates (GR). This is done by scaling maximum growth demands (D) by supply–demand ratios of essential resources (φ* ) (see Eq. (1) and below). GR(t) = (φ (t)D(t)β − Q10 )λ

(A7)

if n = RR is the effect of Bt toxin concentration   if n = Rr, rr

[Bt]/[Bt]max is the index of Bt toxin concentration in host fruit. φ(af ) scales PBW developmental time on fruit of age af (Fig. 1). Similar computations are made for other pest species (Gutierrez et al., 2005).

∗

357

(A6)

A.2.1. Resource acquisition Plants capture light, water and inorganic nutrients and pest larvae attack plant subunits (e.g. leaves, fruit). The biology of resource acquisition by plant and animal consumer (C(t)) involves search under conditions of time varying resource (R(t)) and temperature (T) dependent per capita demand rates (D = D* (t(T))C(t)). Resource acquisition (S(t)) is modeled using the ratiodependent Gutierrez–Baumg¨artner functional response model (Eq. (A7), Gutierrez, 1992) that is a special case of Watt’s model (1959) (see Gutierrez, 1996, for the derivation). α is the search parameter, but it may also be a concave function of C making (A7) a type III functional response (Rochat and Gutierrez, 2001). To simplify, we ignore consumer subscript n, t and T and

Note that stages with preference values equal zero are effectively removed from the calculations. The total age specific consumer demand (D) across ages (i = 1, k) may be computed in mass or number units as appropriate for the population. D=

k  i=1

Di∗ Ci

(A9)

In cotton, the demand rate (g dry matter d−1 ) is the sum of all subunit population demands and includes the costs of conversion of resource to self and respiration (see Gutierrez and Baumgartner, 1984). In pests such as PBW, the demand rate is the maximum per capita adult demand for oviposition sites computed using Eq. (A10) (cf. Bieri et al., 1983). ∗ Di=a = θφ(fruit) φ(Btn ) φT (T )f (a)

(A10)

f(a) = ca/da is the maximum age (a) specific per capita fecundity at the optimum temperature with parameters c and d. θ = 0.5 is the sex ratio. φ(fruit) is a scalar correcting for the effects of host Bt toxin concentrations on adult fecundity (see text). φ(Bt) = 0.5 is the effect of Bt toxins on surviving susceptible adults. φT (T) is the correction

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

358

for the effect of temperature dependent respiration (see (A15)). The sum of age-specific demands for hosts by PBW adults across the n genotype populations is: D=

n  k  i=1

Di∗ Ci .

(A11)

A.2.2. Supply–demand effects Consumer resource acquisition success is estimated by the supply/demand ratio (φS/D (T)) obtained by dividing both side of Eq. (A7) by the population demand D. 0 ≤ φS/D (t) =

S = h(u) < 1· D

(A12)

Some consumers may have multiple resources and they must be included in the computation. We first discuss cotton and then PBW. In cotton, success in meeting its demand is measured by the photosynthate supply/demand ratio (e.g. 0 ≤ φcot,S/D < 1), but there may be shortfall in other requisites. Shortfalls in water (w) are also computed using variants of (A7) (see Gutierrez et al., 1991b). The water 0 ≤ φw = Sw /Dw < 1 ratio is computed in three steps: (i) the potential evapo-transpiration (Dw = PET) and evaporation from the soil surface (ES) are estimated using a biophysical model (Ritchie, 1972); (ii) Dw along with available soil water in the root zone (w = Wmax − Wwp ) above the wilting point (wwp ) are substituted in Eq. (A7) to compute evapo-transpiration (Sw = ET, i.e. water use by the plant); (iii) the input–output model balances rainfall (rain, ES, ET) and runoff or flow through above maximum soil water holding capacity (Wmax ) (see Gutierrez et al., 1988). wwp ≤ w(t + 1) = w(t) + rain(t) − ES(t) − ET(t) ≤ Wmax

(A13)

Similarly for nitrogen, 0 ≤ φη = Sη /Dη < 1 is computed using analogues of Eqs. (A7), (A9) and (A13) (see Gutierrez et al., 1991b for details). The combined effect of shortfall of all essential resources is captured as the product of the independent supply–demand ratios (A14) (i.e. survivorship terms,

Gutierrez et al., 1994). 0 ≤ φ∗ = φ(S/D) φ(w) φ(η) · · · < 1.

(A14)

Eq. (A14) is functionally Liebig’s Law of the Minimum because if any component of φ* causes the supply to fall below a limiting value (e.g. respiration in cotton, see (A6)), it becomes the limiting factor. In plants, after respiration and conversion costs have been subtracted from φ* D, the remaining photosynthate is allocated in priority order to meet demands for reproduction and then vegetative growth and reserves (see Gutierrez, 1992, 1996). In addition to slowing the growth rates of subunits, φ* also reduces the production rate of new subunits, the survival of extant ones (e.g. fruit shedding), and in the extreme may causes the death of the whole plants. The computations for predators (herbivores or carnivores) are similar to those for cotton (see Gutierrez et al., 2005). PBW’s behavior is more akin to that of a parasitoid because individual larvae infest individual fruit. For this reason, the effects of temperature on respiration and hence fecundity are introduced in an approximate way. In poikilotherms, respiration increases with temperature and a plot of the net assimilation rate (supply–respiration) on temperature typically yields a humped or concave function over the range favorable for development with zero values occurring at the lower and upper thermal thresholds and the maximum assimilation rate occurs at Topt . This occurs naturally in models for cotton and predators where mass dynamics are modeled directly. However, to capture the effect of temperature on PBW fecundity, we assume that the normalized effect of temperature on fecundity (i.e. the physiological index for temperature, 0 ≤ φT (t) ≤ 1) is similarly concave (see Gutierrez et al., 1994, Rochat and Gutierrez, 2001). The simplest form for φT is symmetrical (Eq. (A15)).     (T − Tmin ) − γ 2  if Tmin ≤ T ≤ Tmax 1− φT = γ   otherwise 0 (A15) The temperature thresholds for development are Tmin and Tmax , and γ = (Tmax − Tmin )/2 is half the favorable range and Topt = (Tmax + Tmin )/2 is the midpoint. The effect of temperature on PBW oviposition demand (A11) is include as a scalar on total demand

A.P. Gutierrez, S. Ponsard / Ecological Modelling 191 (2006) 346–359

(φT D) which is substituted in random search model (A7) to compute the number of eggs deposited (S) at any level of weighted available fruit (R). A.3. Mortality Mortality in cotton may accrue to whole plants or plant subunits from resource shortages, herbivores, temperature and other factors. Here we focus on PBW. Sources of mortality to PBW are summarized in Eq. (A17), and include temperature, predation on the egg stage, Bt toxin effects on pest larvae, and emigration of adults due to sites for oviposition supply/demand shortfalls. The effect of temperature on PBW survival is also concave and summarized by (A15). Only PBW eggs are exposed to attack by generalist predators, hence we use the 30% d−1 predation rate (µ(pred),i (t)) for noctuid eggs estimated from field data (A.P. Gutierrez, unpublished data), but in Bt cotton, the predation rate is 72% of that rate. Little is known about predation of pupae and adults, and it is ignored in the model. Bt toxin concentrations in fruit kill susceptible genotypes PBW

359

larvae (rr and Rr) at a rate µ(Bt) (af ) per day (see Fig. 1). Susceptible age cotton fruit are shed at a rate proportional to resource shortfalls, and PBW eggs and larvae in them also die (µ(cot S/D),i ). The rate is easy to compute because cohorts of larvae are associated with specific age fruit. PBW females failing to find hosts are assumed to emigrate from the field at a rate proportional to the shortfall (µ(PBW S/D),i ) and hence can be viewed as mortality to the population. Total PBW age specific (i) survivorship (lxi (t)) at time t due to all factors is the product of all age specific survivorship rates (A17). lxi (t)) = 1 − µi (t) = (1 − µ(T ),i (t))(1 − µ(Bt),i (a, t)) × (1 − µ(Pred),i (t))(1 − µ(cot S/D),i (t)) × (1 − µ(PBW S/D),i (t))

(A17)

Note that sources of mortality for some stages may be zero because the stage is not affected. Mortality is incorporated in PBW model (A1) as µi (t).