The invasive Asian citrus psyllid (Diaphorina citri Kuwayama) is the vector of the bacterial pathogen (‘Candidatus Liberibacter asiaticus’) that is the putative causal agent of citrus greening disease (Huanglongbing disease) in citrus in many areas of the world. The capacity to predict the potential geographic distribution, phenology and relative abundance of the pest and disease is pivotal to developing sound policy for their management. A weather-driven physiologically based demographic model (PBDM) system was developed to summarize the available data in the literature, and used to assess prospectively the geographic distribution and relative yield of citrus, the relative densities of the psyllid, its parasitoid (Tamarixia radiata Waterston), and the potential severity of citrus greening disease in North America and the Mediterranean Basin. The potential for natural and biological control of citrus psyllid was examined prospectively.
The Asian citrus psyllid, Diaphorina citri Kuwayama (Hemiptera: Liviidae) is a destructive invasive species causing direct feeding damage to species of citrus and other species in 25 genera of Rutaceae (Shivankar et al. 2000; Halbert & Manjunath 2004). The psyllid is a vector of the bacterium, ‘Candidatus Liberibacter asiaticus’ (and other species of the genus) believed to cause greening disease (Huanglongbing) in citrus (Catling 1970; Bové 2006). Halbert & Manjunath (2004) reviewed the literature on this diseasevector problem.
The psyllid is reported from tropical and subtropical Asia (see Halbert & Manjunath 2004), tropical islands (Reunion, Guadaloupe, Mauritius), parts of South and Central America, Mexico, and the Caribbean (Da Costa Lima 1942; Shivankar et al. 2000; Halbert & Núnez 2004, Étienne & Aubert 1980). Diaphorina citri was first detected in the USA in Florida in 1998 (Halbert 1998; Etienne et al. 2001), Texas in 2000 (French et al. 2001), and southern California in 2008 where there is concern it will spread northward along the coast and into the Central Valley (Grafton-Cardwell 2012). Citrus greening disease was first discovered in the Americas in Brazil in 2004 (Coletta Filho et al. 2004), in south Florida in 2005 (Halbert 2005), and in southern California in 2008 ( http://cisr.ucr.edu/Asian_citrus_psyllid.html).
Successful natural and biological control of heteropterans has occurred (e.g., Pilkinton et al. 2005) and is being pursued for control of D. citri (Hoy et al. 1999; Hoddle & Hoddle 2013). Coccinellid beetles (Harmonia axyridis Pallas, Olla v-nigrum Mulsant, Cycloneda sanguinea L. and Exochomus childreni Mulsant), lacewings (Ceraeochrysa sp. and Chrysoperla rufilabris Burmeister), and a spider (Hibana velox (Becker)) (Araneae: Anyphaenidae) are reported attacking the psyllid in Central Florida (Michaud 2004; Qureshi & Stansly 2009) with O. v-nigrum exhibiting a strong numerical response to psyllid density (Michaud 2001).
Hall et al. (2008) and de León & Sétamou (2010) reviewed the history and effectiveness of the host specific idiobiont ectoparasitoid Tamarixia radiata (Waterston, 1922) (Hymenoptera: Eulophidae) in controlling D. citri. This parasitoid has provided significant control of D. citri on the islands of Reunion and Guadaloupe (Aubert & Quilici 1984; Etienne et al. 2001) with parasitism in Puerto Rico ranging from 79 to 88% (Pluke et al. 2008). In Brazil, levels of parasitism reach 60% (GÓmez-Torres 2009). First releases of T. radiata in California occurred in the southern part of the state in December 2011 (Hoddle & Hoddle 2013).
The goal of this study was to bring the various elements of the biology of this system together in weather driven physiologically based demographic models (PBDM) in a GIS context, and to use the models to assess prospectively the geographic distribution and abundance of the pest in North America and the Mediterranean Basin.
Predicting Geographic Distribution and Relative Abundance
Pivotal for developing sound policy for eradication or management of pests is the capacity to predict their potential geographic range and the relative favorability of sub regions for population growth. Although not established in Australia, Aurambout et al. (2009) developed a stock flow diagram model of citrus flushing (i.e., newly developing clusters of young leaves on the expanding shoot terminals) and D. citri dynamics. They used the model to explore prospectively the potential geographic distribution of the pest in the face of extant weather and projected climate warming due to increasing greenhouse gases (Fischlin et al. 2007). Using simplifying assumptions for mathematical tractability (see Wang & Gutierrez 1980), Chiyaka et al. (2012) developed a mathematical model to characterize the dynamics of the vector and the spread of the pathogen between new terminal shoots (e.g., flushes) within a tree.
Ecological niche models (ENMs) are commonly used to predict the geographic distribution of species, and may be physiological index models (Fitzpatrick & Nix 1968; Gutierrez et al. 1974; Sutherst & Maywald 1985; Sutherst & Bourne 2009), statistical models (see Estrada-Peña 2008; Mitikka et al. 2008), or models based on concepts of information theory (e.g., MaxEnt; Phillips et al. 2006). The MaxEnt approach is gaining considerable attention (Elith et al. 2011) because of its apparent ease of use and its better prediction relative to other ENMs. All ENM approaches attempt to characterize the ecological niche of a species based on aggregate weather and other factors in the recorded range of the species, and then to estimate the potential limits in its recorded range, and the prospective range in new areas where it might invade. ENMs have implicit ecological and mathematical assumptions and lack mechanistic biological underpinnings (Soberón & Nakamura 2009). Similarly, the Fourth Assessment Report of the Intergovernmental Panel on Climate Change Working Group 2 (Fischlin et al. 2007) concluded that ENMs are unable to account for species interactions, and lack physiological mechanism and population processes.
Barlow (1999) reviewed demographic modeling approaches, including the PBDM approach used here to summarize available data and to assess prospectively the distribution and relative abundance of citrus (e.g., yield of navel orange, Citrus sinensis L. Osbeck), D. citri, T. radiata, and greening disease across North America and the Mediterranean Basin. The PBDM system model includes the bottom-up effects of plant growth and development, and the top-down action of natural enemies (Gutierrez et al. 2006, 2007, 2009, 2010a; Ponti et al. 2009). PBDMs capture mechanistically the weather-driven biology of the interacting species to predict their population dynamics and geographic distribution and relative abundance across time and space independent of field distribution records. The PBDM approach resolves many deficiencies of ENMs and other demographic modeling approaches. We note that the PBDM approach has some relationship to the physiological index ENM approach (see Gutierrez et al. 1974, 2010b).
Methods
The PBD System
Our PBDM system is modular consisting of 11 linked functional age-structured population dynamics models {n = 1,…,10} that may be in units of numbers, mass (dry matter) or both. Each of the model is based on the distributed maturation-time dynamics model proposed by Vansickle (1977). The citrus plant model is a canopy model consisting of subunit models for the mass of leaves {n = 1}, stem {2}, shoots {3} and root {4}, and fruit mass and numbers {5, 6}. The models for D. citri consist of age-structured models for eggs {7}, unparasitized nymphs {8}, and the adult stage {10}. The parasitoid model consists of a model for parasitized psyllid nymphs {9} containing parasitoid immature stages, and a model for the freeliving adults {10}. The model for the pathogen is a simple temperature-dependent growth rate model {11} that measures the favorability of temperature for its development. The model for the effects of predation is summarized below.
The underpinning PBDM modeling concepts are reported in Gutierrez & Baumgärtner (1984) and Gutierrez (1992), while the mathematics of the time-invariant distributed maturationtime dynamics model used is found in Vansickle (1977), Gutierrez (1996), DiCola et al. (1999) and Severini et al. (2005) (see appendix). Data in the literature used to parameterize the models are of varying degrees of completeness and quality, and care was required in their interpretation and use. All of the parameters for the model are reported in the text.
Citrus
With variation, the physiological bases and methods for modeling plant growth and development are well established in the literature (see Marcelis & Heuvelink 2007). The PBDM model for olive tree phenology and growth (see Gutierrez et al. 2009) was used as a template for the development of our citrus model designed to capture the phenology and relative yield of navel orange (see initial conditions in Table 1). Our goal was to estimate citrus phenology and quantity of flushing important to psyllid dynamics, and only secondly fruiting and relative yield. The flow of dry matter within the system is illustrated in Fig. 1a–c.
Citrus Developmental Rates. There are many species and varieties of citrus having varying developmental phenology. Data on citrus development across the complete range of temperatures (T) are not available, and estimates of the thermal threshold for citrus vary within a narrow range. For example, Bellows et al. (1985) developed a phenology model for citrus flowering using a lower threshold of 12.5 °C (θL). Hardy & Khurshid (2007) posited that θL = 13 °C and that the upper limit (θU) favorable for growth is 35 °C. Colauto-Stenzel et al. (2006) proposed for ‘Folha Murcha’ sweet orange that θL = 12.8 °C with an upper threshold of θU = 36 °C. Based on this literature, we formulated a normalized developmental rate function using the equation proposed by Brière et al. (1999, eqn. 1) with θL = 12.8 °C. The developmental rate increases to a peak at about 35 °C and then declines to zero near 39 °C.
The maximum degree-days per day (Δtcitrus = dd) above and below θL is computed as, Δtcitrus(T) = 22.2dd × Rcitrus(T), where 22.2 dd is the maximum at 35 °C.
Vegetative Growth, Flowering and Fruit Development. Despite the importance of citrus, the extensive literature is scattered and concerns many species and varieties. Iglesias et al. (2007) reviewed research progress on the physiology of citrus fruiting, and cautioned that information from modeling other plant species cannot always be applied to citrus because of its unusual reproductive biology. Mattos et al. (2003) estimated biomass partitioning and nitrogen allocation in 6-yr old ‘Hamlin’ orange trees. Lovatt (1999) describes the phenology of fruiting in ‘Washington’ navel orange, while García-Luis et al. (2002) estimated the developmental time and biomass accumulation in clementine orange fruit. The data from these sources were used as rough guides for developing the PBDM for citrus.
Leaf flushing and flowering within a tree vary among citrus species and varieties, but are especially prominent during late winter and spring at rates dependent on the availability of new shoots. The presence of fruit reduces summer/fall shoot growth and hence the number of nodes that can develop as spring shoots and bear inflorescences (Valiente & Albrigo 2004; Verreynne & Lovatt 2009). In sweet orange trees, 2.52 to 3.59 flower buds are produced per summer shoot, while only one flower develops per shoot produced during spring (Valiente & Albrigo 2004). We assume complete harvesting, and hence the number of flower buds per tree in year y (i.e., flbuds(y)) is a function of the previous years’ shoot dry matter (gshoot(y-1)) on which buds are formed (see Cohen et al. 2005).
Fig. 1.
The citrus/citrus psyllid/parasitoid/disease system: (a) trophic relationships, (b) energy flow within citrus and to the psyllid, and (c) the biology of the psyllid and its natural enemies (the parasitoid T. radiata and coccinellid beetles).
Flowering is assumed to begin 250 dd after 1 January and last approximately 250 dd. Colauto-Stenzel et al. (2006) estimated that 2500–3600 dd > 12.8 °C are required for fruit maturation in ‘Folha Murcha’ sweet orange. In our model, flowers developing into fruit initially grow slowly for 200 dd and then begin rapid growth that terminates at 2,200 dd, with harvest occurring at age 2,450 dd, usually in late fall. Shedding of flowers and young fruits occurs in response to shortfalls of photosynthate (Iglesias et al. 2003, 2007) (i.e., the ratio of 0 < photosynthate supply/demand<1 in our model). Further, photosynthate production may not be synchronized within the tree allowing vegetative growth (flushing) and the production of new plant subunits until limited by fruit demands. The timing of leaf flushing is important because it affects psyllid reproduction and survival (see below).
Citrus is susceptible to freezing damage (Found 1965) that kills shoots, leaves and fruit, and if severe enough may kill the whole tree. The effect of freezing on citrus at time t (i.e., day) is captured as a survivorship term (lxfrost) by eqn. 3 modified from Gutierrez et al. (2009).
where Tmin is the minimum daily temperature at time t.
Yield in the model was calibrated roughly to data reported by Wheaton et al. (1995) for navel orange in Florida.
Asian Citrus Psyllid
The biology of citrus psyllid and its ectoparasitoid was reviewed by Parra et al. (2010, Fig. 1c). Among the factors estimated from the literature for these species are their temperature thresholds, non-linear temperature-dependent developmental rates and average duration of life stages, maximum per capita age-specific fecundity, temperature-dependent scalars for reproduction and longevity, temperature-dependent mortality rate, and the sex ratio.
Rate of Development of D. citri. The data and function for the psyllid developmental rate for the egg to adult period on temperature (T) are illustrated in Fig. 2a (Liu & Tsai 2000 (
); Nakata 2006 (
)) (eqn. 4, see Brière et al. 1999).
A lower thermal threshold of 12.85 °C was estimated using both sets of data. Data from Parra et al. (2010) suggests a higher threshold of 13.53 °C, and shows the influence of the host on the developmental times. Using the 12.85 °C threshold, average developmental times for eggs, nymphs and adult longevity in dd are 50.4, 155.7 and 435 dd respectively computed in the midrange of favorable temperatures (dd = days × (T - 12.85)). The daily change in dd for development at time t is computed using the egg to adult period as Δte - a(T(t)) = 206.15 Re - a(T(t)).
Reproduction. Adult females prefer to oviposit on new foliage (flushes) where progeny develop to the adult stage (Shivankar et al. 2000). Although adults prefer young leaves, adults may survive for extended periods on mature leaves (Michaud 2004, Pluke et al. 2008). Data on average age-specific oviposition (f(x) = eggs / ♀) at 28 °C (Liu & Tsai 2000, Fig. 2b) was fit using the model proposed by Bieri et al. (1983) (eqn. 5). The pre-oviposition period is 75 dd > 12.85 °C after which f (x) increases nearly linearly to 15 eggs d-1 at age x = 200 dd > 12.85 °C, and then declines to zero (Fig. 2b).
Fig. 2.
The biology of citrus psyllid (a–d) and greening disease (e): (a) the developmental rate per day of Diaphorina citri on temperature [(Liu & Tsai 2000 (); Nakata 2006 ()], (b) the age specific fecundity profile at 28 °C (Liu & Tsai 2000), (c) the effect of temperature on psyllid egg production (Hall et al. 2011), (d) psyllid mortality rate per hour at low temperatures (Tsukuba 2007; Hall et al. 2011; see text), and (e) the reciprocal of time to plant death in potato due ‘Ca. Liberibacter solanacearum’ ( , 
), Munyaneza et al. 2012), and estimates for ‘Ca. Liberibacter asiaticus’ development in citrus (,
) based on observations from Bové (2006) and Lopes et al. (2009).
The total number of eggs (F(t,T)) produced at time t (i.e., day) by the age structured adult psyllid population (N(x,t)) is affected by age (x), sex ratio (sr = 0.5), availability of new leaves (θNL(t)), and temperature (θT(T(t))).
ϕNL(t) is computed as a function of the ratio of new foliage (NL, age < 700 dd) to total leaves (TL).
The effect of temperature on oviposition reported by Liu & Tsai (2000) has a narrower range (i.e., 15–33 °C) than that reported by Hall et al. (2011) for oviposition during a 48 h period (e.g., 16.25 to 40 °C with a peak near 32 °C) (Fig. 2c). We used the Hall et al. (2011) data to estimate the temperature scalar (θT(T(t)) (eqn. 6iii, Fig. 2c).
Temperature Dependent Mortality Rate. The mortality rate of D. citri life stages at low average daily temperatures (Fig. 2d) was estimated from hourly data (µpsyllid (T)h-1) in Tsukuba (2007) and Hall et al. (2011). Converting the mortality rates to a daily basis by multiplying by 24, we get eqn. 7 that declines from unity at approximately -3.2°C.
Mortality rates at temperatures up to 30 °C are low (Liu & Tsai 2000), but further work is required at higher temperatures. (µpsyllid(T)d-1) was also used to estimate mortality rates in parasitized psyllids.
Greening Disease
Greening disease causes yellow mottling of leaves and green banding along the major veins, as well as small-misshaped fruit with bitter taste (Halbert & Manjunath 2004; Chung & Brlansky 2005). The disease is associated with three species of Liberibacter, namely ‘Candidatus Liberibacter asiaticus’, ‘Ca. L. africanus’ and ‘Ca. L. americanus’. Of these,‘Ca. L. asiaticus’ and ‘Ca. L. americanus’ are vectored by D. citri (Munyaneza et al. 2012). ‘Ca. L. africanus’ is heat sensitive, whereas ‘Ca. L. asiaticus’ is heat tolerant. ‘Ca. L. asiaticus’ vectored by D. citri is associated with greening disease in Asia, North and South America, and other parts of the world, with disease symptoms developing under humid conditions at temperatures up to 35 °C (Bové 2006). After 90 days, ‘Ca. L. asiaticus’ infected trees have high disease titers at 32 and 35 °C, but not at 38 °C (Lopes et al. 2009). In a related pathogen, ‘Ca. L. solanacearum’ in potato, Munyaneza et al. (2012) estimated the times to symptom expression and plant death, with plant death at 33.5 and 37.5 °C (symbol ) occurring before symptom expression could develop (Fig. 2e, solid line).
Time to death of citrus from ‘Ca. L. asiaticus’ infection has not been well documented. The model for greening disease progression (eqn. 8) was developed for heuristic purposes and should be viewed as a qualitative index of favorability of temperatures for greening disease development. To develop the model, we extrapolate the disease progression function for ‘Ca. L. solanacearum’ using the upper thermal limits for ‘Ca. L. asiaticus’ estimated from Bové (2006) and Lopes et al. (2009) (Fig. 2e, dashed line).
Parasitoid - Tamarixia radiate
The life history parameters for the T. radiata model were estimated from Gómez Torres (2009) and Gómez-Torres et al. (2012).
Developmental Rates. The lower thermal threshold for egg-pupal stage development is 7.13 °C, and the temperature where the rate of development begins to decline to zero is 36 °C (eqn. 9, Fig. 3a).
Mean egg-pupal (subscript I) developmental time and adult longevity in dd in the favorable range of temperature for egg-pupae development are X1 = 196.6dd>7.13 °C and XA = 212dd>7.13 °C . The daily increment of dd for the parasitoid is computed as ΔtTr(T(t)) = 196.6RTr,I(T(t)).
Reproduction. Female parasitoids oviposit on late instar D. citri nymphs and host feed on early instars. At high host density, a single female may lay 300 eggs. At 25 °C, the pre-oviposition period is approximately 2 days, after which the per capita age-specific fecundity profile declines linearly with age (x in days at 25 °C; Gómez-Torres 2009, Fig. 3b).
Temperature affects fecundity and is captured by a scalar functions θTr(T) in the temperature range 13.5 °C < T < 39.5 °C (eqn. 10ii, Fig. 3c) (Gómez-Torres 2009).
Fig. 3.
Biology of Tamarixia radiata: (a) rate of development per day on temperature for the egg-pupal period, (b) per capita fecundity profile for females of age (x) in days at 25 °C, (c) normalized fecundity on temperature, and (d) the proportion females at 5 day average temperature (data from Gómez-Torres 2009; Gómez-Torres et al. 2012, see also Parra et al. 2010).
The sex ratio (srTr(
), eqn. 10iii) is the proportion of ♂♂ that changes with temperature, and in the model is a function of the running 5-day average temperature (Fig. 3d, eqn. 10iii; Gómez- Torres 2009, Gómez-Torres et al. 2012).
If PTr is parasitoid density, the combined effects of the different factors on reproduction (demand for hosts) by all parasitoid females is FTr(t,T).
Parasitoid Functional Response. A type II parasitoid-form functional response model (Fraser & Gilbert 1976; Gutierrez & Baumgärtner 1984; Gutierrez 1996) that incorporates intraspecific competition in the exponent was used to estimate the total number of hosts attacked (Na) (eqn. 11).
FTr (eqn. 10iv) is the total demand for hosts by all parasitoid females (eqn. 10iv), N = N4 + N5 is the number of available unparasitized 4th and 5th instar hosts, and α (= 0.1) is the assumed area of search calibrated roughly to psyllid densities reported in Florida (e.g., Fig. 4e, Hall et al. 2008). The attacks are apportioned among 4th and 5th instar psyllid nymphs proportional to their density. Skelley and Hoy (2004) found under quarantine conditions that T. radiata parasitized ca. 36% of psyllids and killed 57% by host feeding (ca. 1:1.6 ratio). We assume host feeding on instars 1 – 3 at a rate one per parasitized host with the mortality apportioned by instar according to their density.
Predation
Using exclusion cage experiments, Michaud (2004) found that H. axyridis and O. v-nigrum were the most important biological control agents at high psyllid densities, and that intraguild predation caused > 95% mortality of T. radiata immature stages. The ladybeetle O. v-nigrum exhibits a strong numerical response to psyllid densities (Michaud 2001), and for heuristic purposes, we use a composite mortality rate model to capture the effect of predation on both healthy and parasitized psyllid nymphs (eqn. 12, cf. Gutierrez et al. 2010a).
Predation (0 < µ < 1) increases with both temperature via Δt = Δdd>10°C and the density of healthy and parasitized nymphs, (Nnymphs, N parasitized), but at a decreasing rate. The constant (0.00005) was selected to give roughly the population levels observed by Michaud (2004) and Hall et al. (2008).
Weather Data
Daily weather data for the US and Mexico (1983–2003) and for the Mediterranean Basin (1990–2010) were obtained from the following sources: the global surface summary of daily weather (GSOD) from the National Climatic Data Center, ( http://www.ncdc.noaa.gov/); HelioClim solar radiation data provided by SoDa ( http://www.soda-is.com/); E-OBS gridded observed temperature and precipitation data obtained from the ENSEMBLES ( http://ensembles-eu.metoffice.com) and ECA&D ( http://eca.knmi.nl/) projects; Gauge-Based Analysis of Global Daily Precipitation data sourced from the Climate Prediction Center (CPC) ( http://ftp.cpc.ncep.noaa.gov/); surface meteorology data obtained from the NASA Langley Research Center POWER Project ( http://power.larc.nasa.gov/) funded through the NASA Earth Science Directorate Applied Science Program; Daymet data ( http://www.daymet.org/). The few missing records in the data set (<0.0005% of daily records) were estimated by linear interpolation within the data.
Simulation and GIS Mapping
The PBDM system model is modular, and Boolean variables in a setup file were used to determine the combinations of species included in the different runs and to run the model across locations for the period of available weather data. The same initial conditions for citrus (Table 1), and psyllid population densities (3 eggs, 8 nymphs, 0.6 adults) and parasitoid (0.25 immature, 0.25 adults) were used for all locations with weather determining the resulting dynamics.
Numerous life-history variables are computed daily for each species, and hence we must chose standard metrics across locations (e.g., total nymph-days y-1 for D. citri). Output variables from the model were geo-referenced and written by year to batch files. Means, standard deviations and coefficients of variation for all output variables were computed across yrs for each location, and the results were written to batch files for mapping. The system was assumed equilibrating to local conditions during the first year, and hence these data were not used in the summary analyses. A 20 yr run at one location requires about 10–15s on a laptop computer.
The geo-referenced simulation data were written to files by yr and mapped for elevations below 1500 m using the open source Geographic Resources Analysis Support System (GRASS) GIS software (GRASS Development Team 2012). Inverse distance weighting interpolation of the data via the v.surf.idw GRASS module was used in mapping. The patterns in the maps reflect not only the site-specific effects of weather on the dynamics of the system, but also the location and distances between weather stations. Predicted densities should be viewed as indices of relative favorability and not precise predictions of population density.
Results
We first review the dynamics of the system to illustrate some of the richness of the model output, and then examine the summary dynamics and distribution prospectively across the US and Mexico. We use Visalia, California as a reference location.
System Dynamics at Visalia, California
Citrus. Fig. 4a illustrates the simulated dynamics of fruit and leaf mass for the period 1983–2003, while Fig. 4b shows the cumulative dd below -3.5 °C that cause mortality to fruit, shoot and leaf mass, as well as the psyllid and its parasitoid. For example, the winters of 1987, 1989, 1990 and 1999 had freezing periods that affected spring growth and flowering. The severe 1990 cold period affected the citrus growth dynamics that year and during 1991 with some carryover into 1992 (i.e., the stippled vertical bar).
Psyllid Dynamics. Oviposition occurs during periods of leaf flushes in citrus that usually occur during spring. Flushing may also occur when the fruit demand is reduced due to harvest or maturation of fruit, and at other times in response to rain following drought. Psyllid populations decline as foliage matures and becomes less suitable for reproduction (eqn. 6ii).
In the absence of natural enemies, simulated D. citri populations cycle with flushes reaching very high levels (not illustrated). Including the action of predation by coccinellid beetles reduces psyllid densities 8–10 fold yielding the dynamics patterns in Fig. 4c. Cold weather, especially frosts, cause steep declines in psyllid densities such as occurred during late 1990 (Fig. 4c, see also Tsukuba 2007) with populations rebounding in 1991 in response to tree regrowth and favorable weather.
Fig. 4.
The simulated dynamics of the citrus/Diaphorina citri system at Visalia, California for the period 1983 to 2003: (a) leaf and fruit dry matter growth dynamics (g), (b) cumulative degree days below -3.5C, and D. citri dynamics: (c) with only predation, (d) with only parasitism, and (e) with both predation and parasitism. The stippled vertical area indicates a period following freezing temperature that affected all species in the system (see text).
Including the parasitoid T. radiata alone results in periodic outbreaks of the psyllid with increased numbers of adults relative to nymphs (Fig. 4d). In contrast, the combined action of both coccinellid predation and parasitism is predicted to cause a further 8 to 10 fold decrease in psyllid density compared to predation alone (Fig. 4c vs. Fig. 4e). The difference between the simulated effects of predation plus parasitism, and parasitism alone is a measure of the indispensable mortality and stabilizing effect of ladybeetle predation (Fig. 4d vs. Fig. 4e; see Michaud 2004; Qureshi & Stansly 2009).
Regional Analyses - North America.
Fig. 5 shows regional maps for relative citrus yield (g dry matter per tree, Fig. 5a), the distribution of favorability for greening disease as an index (0 ≤DI ≤1, Fig. 5b), and D. citri nymphal density (total nymph-days y-1) with predation (Fig. 5c) and plus the added effects of T. radiata (Fig. 5d). The product of the disease index (DI, Fig. 5b) and normalized psyllid density (0 ≤CPI ≤ 1) computed from data in Fig. 5d is illustrated in Fig. 5e (i.e., DI x CPI). The results shown in Fig. 5 are also displayed with greater clarity in color as Suppl. Fig. 5 in Supplementary Material for Florida Entomologist 96(4) (2013) online at http://purl.fcla.edu/fcla/entomologist/browse.
Fig. 5.
The simulated average dynamics of the citrus/citrus psyllid/greening disease system across North America below 1500m: (a) simulated yield in g dry matter per tree, (b) the normalized growth rate of greening disease (DI), (c) psyllid nymph days per tree given the effects of coccinellid predation, (d) psyllid nymph days per tree given the added effects of Tamarixia radiata parasitism, and (e) the product of DI and the normalized psyllid density (i.e., CPI) computed from the data used in Fig. 5d
Citrus Yield. Highest average yields are predicted in tropical areas of Mexico, Arizona-California and south Florida. Commercial citrus production in California occurs in the southern coastal region, the lower eastern areas of the Central Valley and desert areas of southern California, but citrus is also widely grown in private gardens across California (e.g., Berkeley, California). Citrus in arid areas of Arizona - California and Mexico requires irrigation that in some cases is not available.
Citrus Greening Disease. Favorability for citrus greening disease is measured as the normalized index of the cumulative daily pathogen developmental rates (eqn. 8) per year (0 ≤ DI ≤ 1). The results suggest that the disease organism has a wider prospective distribution than citrus (comparing areas in the upper 2 quartiles, Fig. 5a vs. 5b). The disease is limited northward by cold temperatures that also limit citrus and the psyllid vector.
Psyllid Dynamics. With predation by native coccinellid beetle adults and larvae, mean D. citri densities range from 80 to 14,520 nymphal-days y-1 with intermediate densities predicted in marginal areas of citrus production across the southeastern USA and Texas (Fig. 5c). The warmer mostly frost-free areas of Arizona, California, Texas, Florida and Mexico have higher densities.
With the addition of parasitism by T. radiata, psyllid densities are reduced approximately 8–10 fold across the region. Areas marginal for citrus in the SE USA and some areas of south Texas, south Florida, and tropical Mexico have predicted high psyllid densities due to adverse effects of high temperatures on the parasitoid (fig. 5c vs. 5d).
Effects of Natural Enemies on Disease Distribution. The prospective distribution of greening disease given the action of predation and parasitism on the vector is computed as the product of the disease index (0 ≤ DI ≤ 1, Fig. 5b) and the normalized cumulative nymphal-days computed from the data in Fig. 5d (0 ≤ CPI ≤ 1). Prospectively, in order of decreasing favorability (Fig. 5e), locations in the upper risk quartile (DI x CPI = 0.75-1) include the Yucatan and south Mexico which is favorable for both (fig. 5c, d). Quartile 0.5 – 0.74 includes south Texas, eastern coastal Mexico, parts of Baja California, and central and south Florida. The 0.25 – 0.49 quartile encompasses much of central California, Arizona and northern areas of western Mexico across to Texas and the SE USA, while very marginal areas (0–0.24) include much of near coastal California. Color figures are provided as supplemental materials.
Regional Analyses - The Mediterranean Basin
The model captures the potential citrus growing region in the Mediterranean Basin in the upper half of the distribution of yield quite well (Fig. 6a). The prospective distribution of favorability for citrus greening disease is greatest in North Africa, southern Spain and the Middle East (Fig. 6b), with the countries in southern Europe, including Turkey being intermediate in favorability. The prospective distribution of citrus psyllid is more restricted (Fig. 6c) than that of greening disease, and to estimate the prospective joint distribution we again multiply normalized psyllid days (i.e., from Fig. 6c) and the disease index (DI x CPI, Fig. 6d). The results shown in Fig. 6 are also displayed with greater clarity in color as Suppl. Fig. 6 in Supplementary Material for Florida Entomologist 96(4) (2013) online at http://purl.fcla.edu/fcla/entomologist/browse.
Fig. 6.
The simulated average dynamics of the citrus/citrus psyllid/greening disease system across the Mediterranean Basin below 1000m: (a) simulated yield in g dry matter per tree, (b) the normalized cumulative growth rate of greening disease (DI), (c) psyllid nymph days per tree given the effects of coccinellid predation and Tamarixia radiata parasitism, and (d) the product of DI and normalized psyllid density (i.e., CPI) computed from data used in Fig. 6c.
Prospectively, the joint favorability suggests the eastern Mediterranean region is at greatest risk with only Sicily and small areas of southern Spain included in the upper half of the range.
Discussion
Asian citrus psyllid and greening disease occur across a wide area of Asia, the Caribbean region and South America, and both have the capacity to invade all of the citrus growing regions of North America and the Mediterranean Basin where the psyllid and the disease are not yet reported. Eradication of the psyllid pest is unlikely, and hence the goal must be to reduce or disrupt the tripartite interactions of diseased trees availability, psyllid densities and disease transmission (sensu Nakazawa et al. 2012). Like citrus, the geographic range of the psyllid is limited by cold winter temperatures with both having similar lower thresholds (12.8 and 12.85 °C respectively), with seasonal patterns of vegetative growth (flushing) being an important driver of D. citri phenology and dynamics (Hall et al. 2008, Chong et al. 2010). As modeled, the greening disease bacterium appears to have a wider tolerance to temperature (lower and upper thermal thresholds, 10.1 °C and ∼38 °C respectively) than its psyllid vector (12.85 °C – 37 °C). Full data on disease progression in citrus need to be developed to enable a more accurate assessment of favorability for ‘Ca. L. asiaticus’ across areas of citrus culture.
Laboratory studies by Pelz-Stelinski et al. (2010) provide valuable insights on greening disease transmission by Asian citrus psyllid, and are summarized below. For example, they found that nymphs reared on ‘Ca. L. asiaticus’ infected plants were more likely to acquire the bacterium than adults, with disease acquisition ranging from 60 to 100% during nymphal development, but only 40% in adults after 5 wk of feeding. Similar rates of pathogen acquisition were observed in the field. Transovarial transmission of the bacteria from parent to offspring occurs at a rate of 2–6%. Adults acquiring the pathogen as nymphs are much better vectors than those acquiring the pathogen as adults. Transmission of the disease increases with vector density being 4 to 10% by individuals whereas groups of 100 or more transmit the pathogen at a rate of ∼88%. The epidemiology of the disease also depends on temperature (Bové 2006; Lopes et al. 2009) and vector density, especially of adults that move between hosts.
Gilligan & van den Bosch (2008) reviewed progress on the analysis of plant epidemiology from a theoretical perspective. Nakazawa et al. (2012) developed a theoretical model integrating the tripartite interactions among host plant, plant pathogen and herbivore vector in a community context, and posited that as well as efficiency of disease transmission, mass effects (i.e., abundance of infected trees) are crucial for disease persistence.
Prior studies have not included the effects of natural enemies in reducing vector density. Here, the product of normalized cumulative nymphal-days given the action of predation and parasitism (0 ≤CPI ≤ 1; Fig. 5e) and the normalized annual cumulative disease growth rate (0 ≤ DI ≤ 1, Fig. 5d) were used as a measure of tripartite favorability for the disease. In the interaction term (0 ≤ DI x CPI ≤ 1), lower favorability for either or both factors decreases favorability for disease prevalence. For example, central California and the upper margins of the Gulf States are favorable for the vector but less favorable for the disease resulting in a low joint favorability. Conversely, Baja California is quite suitable for the disease but less favorable for the vector, again lowering the joint favorability for the disease. In the Mediterranean Basin, the region with highest joint favorability is the eastern Mediterranean (Fig. 6d).
The most obvious way to control greening disease is to lower vector density, and the only viable option is effective biological control as occurred for glassy winged sharpshooter (Homalodisca coagulata (Germar)) (Pilkington et al. 2005) that is the vector of Pierce's disease. In that case, high joint favorability for the disease and vector occurs only in the hot desert regions of southern California (Gutierrez et al. 2011).
In the absence of natural control, very high average psyllid densities are predicted that maximize the likelihood for disease transmission. Adding coccinellid predation, psyllid densities decline 8–10-fold (e.g., Fig. 5c) across the region. Reported levels of parasitism are quite variable with good control reported in Reunion and Guadalupe (Aubert & Quilici 1984; Etienne et al. 2001), levels of ∼56% in Florida (Qureshi et al. 2009, Chong et al. 2010), 79 to 88% in Puerto Rico (Pluke et al. 2008), and ∼60% in Brazil (Gómez-Torres 2009). The parasitism rate reported by Pluke et al. (2008) is similar to that predicted by our model for Visalia, CA using parameters estimates for the ecotype of T. radiata in Brazil (Gómez Torres 2009), and the search parameter (α = 0.1, eqn. 11) calibrated roughly to psyllid densities reported in Florida (e.g., Fig. 4e, Hall et al. 2008). Adding T. radiata to the system model reduces psyllid densities a further 8–10 fold, but the results are not uniform across the region. Highest psyllid densities are predicted in tropical Mexico and south Florida (Fig. 5c, d) where relatively high DI x CPI are also predicted.
Natural enemy introduction remains an inexact science, though modeling can provide insights as to why some species are efficacious and other are not, say due to differences in search related behavior (see Gutierrez et al. 2011). For example, a 10-fold reduction in the parasitoid search parameter α predicts increases in the average of peak psyllid density of ∼3-fold at our reference location of Visalia, California. However, other factors may also influence natural enemy efficacy. For example, McFarland & Hoy (2001) found differences in survival of 2 ecotypes of T. radiata and the endoparasitoid Diaphorencyrtus aligarhensis (Shafee, Alam & Agarwal) under different humidity (and temperature) conditions, and this suggests that a priori selection for the best parasitoid or ecotype seems unlikely. Hence, multiple introductions may be required to allow self-selection to occur across the ecological regions of citrus culture.
Acknowledgments
Special thanks are due Drs. M. L. Gómez-Torres, D. E. Nava and J. R. P. Parra for providing a copy of the Doctoral Thesis on the biology of T. radiata. J. P. Michaud and P. A. Stansly provided early reviews and helpful comment. We continue to be grateful to Dr. M. Neteler (Fondazione Edmund Mach - Centro Ricerca e Innovazione, Trento, Italy; http://gis.fem-environment.eu/) and an international network of co-developers for maintaining the Geographic Resources Analysis Support System (GRASS) software, and making it available to the scientific community.
No extramural funding for this work came from US sources. This research was supported entirely by the Center for the Analysis of Sustainable Agricultural Systems, Kensington, California ( http://cnr.berkeley.edu/casas/), and by a Marie Curie International Reintegration Grant within the7th European Community Framework Program (project number: 224091; acronym: GlobalChangeBiology).
References Cited
Appendices
APPENDIX
OVERVIEW OF THE DISTRIBUTED-MATURATION-TIME DYNAMICS MODEL
The physiologically based demographic model (PBDM) approach builds on the idea that all organisms are consumers and all have similar problems of resource acquisition (inputs) and allocation (outputs) (Gutierrez 1992). This notion allows the use of the same resource acquisition model and birth-death dynamics models to describe all trophic levels including the economic one (Regev et al. 1998). For consumers such as citrus and citrus psyllid, the model for per capita resource acquisition (i.e. the supply, S) is driven by total consumer demands for dry matter (D) in priority order for egestion, conversion costs, respiration (i.e., the Q10 rule in ectotherms), and reproduction, growth, and reserves. The ratio 0 ≤ S ≤ D ≤ 1 is always less than unity due to imperfect consumer search, and is used in the model to scale maximal vital rates of species. Parasitoid seeks hosts (e.g., nymphs) as outlined in the text.
The Erlang distributed-maturation-time demographic model is now widely used to simulate the age structured population dynamics of species (Vansickle 1977 and DiCola et al. 1999, p. 523–524). The numerical solution for the time varying form is found in Severini et al. (2005). The numerical solution for the time invariant form for ith age class of a life stage with i = 1, 2,…, k age classes is equation A1. The time step is a day (t) with forcing variable being temperature that from the perspective of the poikilotherm organisms is of varying length in physiological time units (degree-days, dd).
The state variable ri(t) is the number density (it could also be mass) as a rate in age class i,
is the flux rate between age classes, X is the mean developmental time in dd, and Δx(t(T)) is the dd increment of physiological age (x) at time t and temperature T (see below). µ1 is the proportional net loss rate that includes the rich biology affecting age class rates of deaths, growth, predation, and net immigration. Births (r0) flow into the first age class (k=1) as A · (r0(t) - ri(t)), and some individuals exit at maximum age (k) as A · rk(t). The number in an age class is and the total population is
. Because field data on an area basis were unavailable, density was modeled on a per tree basis.
The distribution of stage maturation times is determined by the number of age classes k and the variance of maturation times (var) (e.g., k = X2/var). The larger the value of k, the narrower is the Erlang distribution of developmental times of cohort members. A value of k =25 was selected that assumes an intermediate level of variability in developmental times for all species and substages.
Table 1.
Initial Parameters For The Citrus Model
Notes
[1] Supplementary material including color photos for this article in Florida Entomologist 96(4) (2013) is online at http://purl.fcla.edu/fcla/entomologist/browse
