The potential geographical distribution and phenology of Bemisia tabaci Middle East/Asia Minor 1, considering irrigation and glasshouse production

The Bemisia tabaci species complex is one of the most important pests of open field and protected cropping globally. Within this complex, one species (Middle East Asia Minor 1, B. tabaci MEAM1, formerly biotype B) has been especially problematic, invading widely and spreading a large variety of plant pathogens, and developing broad spectrum pesticide resistance. Here, we fit a CLIMEX model to the distribution records of B. tabaci MEAM1, using experimental observations to calibrate its temperature responses. In fitting the model, we consider the effects of irrigation and glasshouses in extending its potential range. The validated niche model estimates its potential distribution as being considerably broader than its present known distribution, especially in the Americas, Africa and Asia. The potential distribution of the fitted model encompasses the known distribution of B. tabaci sensu lato, highlighting the magnitude of the threat posed globally by this invasive pest species complex and the viruses it vectors to open field and protected agriculture.

Species distribution modelling methods employed by Herrera Campo et al. (2011) resulted in some very poor projected distributions. Herrera Campo et al. (2011) modelled the potential distribution of B. tabaci sensu lato based solely on distribution data, using an ensemble method training five models (climate space, environmental distance, GARP, MaxEnt and support vector machines) with four sets of covariates. The results were highly erratic, ranging from excessively conservative to excessively liberal, and often heavily biased. Even the models that were very liberal had poor sensitivity. The environmental distance trained with all covariates had the highest sensitivity score of 0.941, meaning that it was unable to account for 6% of records, and yet indicated (implausibly) that Central Australia, the Sahara Desert and the entire wet tropics were all modelled as highly suitable. Furthermore, six of the 16 model results indicated that Australia was completely unsuitable for B. tabaci, whilst the remaining ten models indicated suitability in the central desert regions of Australia where cropping could only be undertaken where cropping does not currently occur, but could be undertaken if irrigation was possible. The results of an elaborate weighting scheme in Herrera Campo et al. (2011) indicate climate suitability in areas that are too dry (e.g., Sahara Desert, Namibian Desert and Central Australia) and excessively cold (Southern Argentina and mid-west USA). Clearly, ensemble methods do not, as is often claimed (e.g., Araujo and New 2006), necessarily correct for poor performing models (Hannemann et al., 2015).
Whilst the present known distribution of B. tabaci sensu lato is broad, the known distribution of B. tabaci MEAM1 is more confined. For example, in the USA, it has a northern limit at Tulare in the San Joaquin Valley in California, and greater numbers at Bakersfield at the south end of the valley (Gruenhagen et al., 1993). It is abundant at Brawley in the Imperial Valley (Zalom et al., 1985), and also occurs in smaller numbers at Palm Springs and Phoenix (Arizona), but is absent from New Mexico and Los Angeles (California). In the Texas, it is abundant around Brownsville and the Lower Rio Grande Valley and occurs as far north as Del Rio, but cannot survive in San Antonio or along the northern rim of the Gulf of Mexico. It occurs in southern Florida, especially on crops around Immokales near Fort Myers; further north it is absent from areas north of Orlando, except for a pocket at Charleston (South Carolina) (McKenzie et al., 2004).
In this paper, we develop a climate suitability model to estimate the potential distribution of B. tabaci MEAM1 (De Barro et al., 2011). We compare the fitted model to the known distribution of another invasive member of the complex, B. tabaci MED, which has been noted invading Canada, China, Guatemala, Japan, Mexico, the Netherlands, New Zealand, South Korea, Uruguay and the United States (McKenzie et al., 2009;Qiu et al., 2009).
The Australian distribution of B. tabaci was determined by a series of surveys across Australia between July 1996 and June 1997. Surveys involved searching available known hosts in a particular location for a maximum of 20 min. The pest status of B. tabaci was determined on a location-specific basis through grower surveys and consultation with key grower representatives. B. tabaci MEAM1 was considered to hold pest status in a particular location if growers were taking management decisions to minimize numbers. Identifications were carried out using randomamplified polymorphic DNA-polymerase chain reaction (De Barro and Driver, 1997).

Development rate studies
The development and fecundity of B. tabaci MEAM1 and two populations of Australia (B. tabaci EAN and WAN, formerly B. tabaci biotype AN) were tested experimentally. The B. tabaci MEAM1 culture was established from individuals collected from tomato in Bundaberg, Queensland; B. tabaci EAN was established from collections made on Euphorbia cyathophora in Bundaberg, Queensland and WAN from Emilia sonchifolia in Kununurra, Western Australia. Single-mated females of the B. tabaci MEAM1, Eastern and Western accessions of the Australia biotypes were allowed to oviposit for 24 h onto E. cyathophora seedlings kept in cylinder cages at 15, 20, 25 or 30°C. Each treatment was replicated ten times per temperature. After 24 h, the insects were removed. The time taken for development from newly laid egg to emerging adult was noted.

Modelling
The CLIMEX modelling package (Sutherst and Maywald, 1985;Kriticos et al., 2015) was used to fit an ecological niche model to the known distribution of B. tabaci records identified as B. tabaci MEAM1 (De Barro et al., 2011). CLIMEX is a climate-response model that produces indices derived from the responses of a nominated species to each of the components of climate. Rainfall and relative humidity or evaporation are combined into a single soil moisture index (MI), using a hydrological model to represent the availability of moisture. The effects of extreme climatic values are taken into account in a series of stress indices that estimate the threat to B. tabaci posed by prolonged adverse periods of excessively cold, hot, dry or wet conditions. Finally, the growth and stress indices are combined into an ecoclimatic index, scaled from 0 to 100, to represent the overall favourability of the given geographical location for the permanent survival and propagation of populations of B. tabaci. The annual minimum heat sum for development (in degree-days) is used to estimate the expected number of generations in a year.
B. tabaci MEAM1's climatic requirements were inferred from information on its known global geographical distribution, relative abundance and seasonal phenology. Some life cycle data, such as developmental threshold temperatures, were used to finetune or interpret the CLIMEX parameter values and to provide more confidence in these values (table 1). However, published reports on development rates and population dynamics of B. tabaci MEAM1 provide a wide range of parameter values for lifecycle processes, largely dependent on the host plant (Zalom et al., 1985;Zalom and Natwick, 1987;Van Giessen et al., 1995;Tsai and Wang, 1996;Liu and Stansly, 1998;Muniz and Nombela, 2001;Lin and Ren, 2005;Musa and Ren, 2005). CLIMEX parameters were fitted using three assumptions about the conditions under which B. tabaci can exist: (i) rainfed field environment, (ii) irrigated crop environment (top-up irrigation to produce 25 mm week −1 in summer, 15 mm week −1 in winter) and (iii) glasshouses enabling unrestricted over-wintering (cold stress removed).
Published data on geographical distribution and abundance do not clearly define the role of glasshouses in maintaining populations reported in some colder regions. This limited the degree of confidence that we could attribute to estimates of over-wintering ability as determined by low temperatures. When estimating geographical variation in the potential numbers of generations, we assumed that B. tabaci MEAM1 had either over-wintered successfully in very cold locations, perhaps in a greenhouse, or was re-introduced in the spring each year. We made no correction for the modifying effects of crop canopies on microclimates because the target area included both crop-and non-crop host plants, with widely varying canopy densities.

Model fitting, verification and validation
The CLIMEX growth parameter values were adjusted to describe seasonal growth and over-wintering of B. tabaci at Brawley, California, where data on B. tabaci show high population levels under irrigation (Coudriet et al., 1985), at Fresno, northern California, where B. tabaci MEAM1 has been reported in smaller numbers (Gruenhagen et al., 1993) and at Katherine, Australia where it is only marginally suitable (De Barro, unpublished data). The model was fitted using distribution data for the Middle East, Asia Minor and the USA, and verified using data from Australia and Eastern Asia. Distribution data in Africa, Central America, the Caribbean, Europe and South America were used to validate the model. The model-fitting procedure involves an implicit cross-validation between distribution data, parameters derived from observed phenologies (e.g., Zalom et al., 1985;Zalom and Natwick, 1987) and parameters estimated from laboratory experiments (e.g., Butler et al., 1983).

Soil moisture index
The lower soil moisture limit (SM0) was set at 0.1, which approximates the permanent wilting point of many crop plants (Kriticos et al., 2003); reasoning that if the soil moisture is insufficient for host plants to grow, then B. tabaci MEAM1 is also unlikely to be able to grow. The optimal range for growth was set from 0.5 to 1.5 to cover the range of soil moisture levels at which most crops are able to grow optimally. The upper soil moisture limit, SM3, was set to 2, accounting for the ability of B. tabaci MEAM1 and its hosts to grow to a limited capacity under fairly wet conditions such as those found in Karnataka in south-western India.

Temperature index
The base temperature for development (DV0) was set to 12°C, which compares favourably with empirical findings from Zalom and Natwick (1987) fitted to population peaks of a base temperature of 10°C. The lower temperature limit for optimal growth (DV1) was set to 28°C and the upper optimum (DV2) was set to 32°C in accordance with results from growth rate studies (Zalom et al., 1985;Zalom and Natwick, 1987;Van Giessen et al., 1995;Tsai and Wang, 1996;Liu and Stansly, 1998;Muniz and Nombela, 2001;Musa and Ren, 2005). Butler et al. (1983) noted that eggs failed to hatch at a constant 36.0°C, and that aestivation occurred with fluctuating temperatures between 27 and 43°C. Accordingly, the upper limit for growth (DV3) was set to 42°C.
The development rates for B. tabaci MEAM1 and the two Australia accessions showed no significant differences in the slopes of the fitted temperature development response functions, and differed only marginally in the development rates for a given temperature ( fig. 1; De Barro and Hart, 2000). These results suggest that the optimal temperature for development is likely to be above 30°C ( fig. 1), lending further support to the choice of parameters for the optimum temperature range.

Dry stress
The dry stress threshold (SMDS) was set at 0.1, reasoning that if the host plants were unable to grow due to drought-stress, then B. tabaci MEAM1 populations would likely decline also. The dry stress accumulation rate was fitted to barely enable persistence in the locations in the Middle East where it has been reported (Israel, Jordan and western Syria). B. tabaci MEAM1 has been reported in extremely xeric locations such as Karameh in southern Iran and the cropping areas in El Centro, Southern California and Maricopa, Arizona in the USA. These locations depend on irrigation for agriculture, and hence the fitted dry stress under the natural rainfall scenario indicates that these locations are unsuitable for year-round persistence.

Wet stress
The wet stress threshold SMWS (2.0) and stress accumulation rate HWS (0.002 week −1 ) were fitted to preclude B. tabaci MEAM1 from Goa in western India, while still allowing it to persist in Karnataka where it has been reported frequently.

Hot-wet stress
An equipment failure provided an opportunity to observe that under laboratory conditions, B. tabaci MEAM1 populations thrived under hot dry conditions, but declined at similar temperatures with high humidity (De Barro, pers. obs.). The hot-wet stress parameters were adjusted to make Katherine in the Northern Territory of Australia barely suitable.

Generation time
In the literature, the criteria used to define the generation time have not been standardized, making comparison of different datasets difficult (Zalom et al., 1985;Zalom and Natwick, 1987;Van Giessen et al., 1995;Tsai and Wang, 1996;Liu and Stansly, 1998;Muniz and Nombela, 2001;Musa and Ren, 2005). Further, as generations quickly overlap each season, and there is significant variation in timing of emergence of the first and last individuals from each generation, estimates of the number of generations is necessarily an approximation. To account for this variation, we investigated the sensitivity of CLIMEX projections to changes in the thermal sum (number of degree-days,°C days) required to complete a generation. We examined the developmental threshold temperature (10, 12 and 14°C) and the thermal sum (330, 365 and 400°C days) required to complete one generation. The generation time was taken as the sum of the pre-oviposition period and the development times of each lifecycle stage, to produce 50% emergence of adults. The fitted value for PDD of 365°C days above DV0 (12°C) compares favourably with the value of 369.5°C days above 10°C noted by Zalom and Natwick (1987).

Simulations
The model was run with the CliMond database of long-term average climate data on a 30 ′ grid (Kriticos et al., 2012). It consisted of a 30-year average of monthly minimum and maximum temperatures, relative humidity at 09:00 and 15:00 h and total precipitation. This climatology (CM30_1995H_V2) was centred on 1995, spanning the period 1980-2010. CLIMEX interpolates the monthly average temperatures to weekly values in order to describe the daily temperature cycle for each week. The initial modelling was undertaken without irrigation, and locations records that could not be fitted with reasonable dry stress and MI values were scrutinized using Google Earth to identify telltale signs of irrigation (e.g., green circles or rectangles amongst a background of brown or yellow exposed soil and the nearby presence of a river or dam). Simulations were then undertaken with 2.1 and 3.6 mm day −1 top-up irrigation in winter and summer, respectively. The Global Map of Irrigated Areas Version 5 (Siebert et al., 2005;Portmann et al., 2010) was used to define areas where the irrigation and natural rainfall scenarios were included in a final composite risk map. The GMIA dataset included some dubious artefacts in desert regions, including Central Western Australia, Algeria and Libya. In these regions, the dataset indicated that irrigation was being practiced over small, proportions of the landscape. Searching these areas using Google Earth did not reveal any signs of irrigation. Consequently, we set a threshold of 5 ha per 10 ′ grid cell in GMIA V5 to indicate that irrigation was being practiced. This is equivalent to approximately 0.05%.

Known distribution and pest status
As befits the reputation of the B biotype as an invasive species, the known distribution of B. tabaci MEAM1 is extensive, being present on every continent except on Antarctica ( fig. 2) . 2). In China, both species are found far to the north of their modelled potential distribution, as far north as Beijing. We attribute this to B. tabaci MED and MEAM1 probably overwintering in greenhouses, which are abundant in this region and easily recognized using Google Earth. In the absence of cold stress, the annual growth index (GI A ) in this region is adequate for strong population development, which accords with observations that it is a major pest in this region (Ren et al., 2001). B. tabaci MEAM1's pest status in Australia seems confined to the sub-tropics, with populations further south only appearing to warrant pest control in greenhouse situations.

Bulletin of Entomological Research 5
The known Australian geographical distribution is shown in fig. 2, as are the areas where it is known not to occur. In Western Australia, B. tabaci MEAM1 occurs at Carnarvon, but has yet to establish in Kununurra. In the Northern Territory, it has been recorded as far north as Darwin, and while it is present in the Darwin rural area, it has failed to persist further inland in Katherine, despite having initially established there. In Queensland, it is found as far north as Cooktown and extends southwards in a peri-coastal band. In the south-east of the State it extends further inland, spanning the Darling Downs as far west as St George and north to Roma, Emerald, Richmond and Dimbulah. On the east coast of New South Wales, B. tabaci extends as far south as Bateman's Bay and inland to Bourke, the Sydney Basin, Namoi (Narrabri) and Gwyder (Moree) Valleys, Narromine, Tamworth and Warren. In New South Wales, numbers tend to be low except in the Tweed Valley region around Cudgen and the Sydney Basin. In the locations south of Cudgen, infestations tend to be in urban areas in close proximity to dwellings; however, in the Sydney Basin, field infestations are in close proximity to glasshouse infestations. Thus, in areas south of Cudgen, including the Sydney Basin, the protected environments afforded by urbanization and glasshouses almost certainly confound local climatic suitability. B. tabaci MEAM1 has also been recorded from glasshouses in South Australia, New South Wales, Victoria and Western Australia, usually on plant material imported from Queensland and New South Wales. These data indicate that the species is well suited to and distributed in the tropics, but is less well adapted to the cooler non-tropical regions and its presence may well be supported by the provision of overwintering sites in protected cropping and urban settings.

Phenology
In CLIMEX, the stress indices have the strongest impact on the species modelled range, whereas the growth indices have a stronger effect on the modelled suitability patterns within the range  (Siebert et al., 2005;Portmann et al., 2010) and (d) comparison between natural rainfall and irrigation scenario. The irrigation scenario is described in the main text.

Bulletin of Entomological Research 7
boundaries, and also drive directly the modelled phenological patterns (Kriticos et al., 2015). The modelled growth indices are presented in fig. 3 for three locations where the phenology is known. Naranjo et al. (2010) observes that the impacts of B. tabaci in the San Joaquin Valley (near to Fresno) are muted compared with the impacts observed in the Imperial Valley (Brawley), attributing this to the warmer climate in the southern site, especially during the early part of the season. The modelled temperature and growth indices accord with these observations ( fig. 3). At Katherine, in the Northern Territory, the apparent strong potential for growth ( fig. 3) is offset by hot-wet stress. Here, it has been observed to establish during the cooler, drier months and undergo local extinction during the warmer, wetter months (Naranjo et al., 2010).

Potential distribution and pest status
The modelled global potential distribution of B. tabaci MEAM1 reflects the belief that it is a species adapted to the dry tropics and sub-tropics. Its modelled climate suitability under natural rainfall conditions throughout much of its putative native range in the Middle East is marginal or poor ( fig. 4a), but increases substantially in those areas where irrigation is practiced (fig. 4b). The composite irrigation scenario provides a pest risk map that highlights those areas that may be at risk if suitable hosts are grown ( fig. 4c). A difference map between the natural rainfall and the composite irrigation scenarios is presented in fig. 4d. The contrast in suitability under irrigation highlights B. tabaci MEAM1's adaptation to hosts that are ruderals sensu Grime (1979), and its ability to take advantage of habitats that become suitable for only short periods, by having a short lifecycle and high per capita rate of natural increase and a durable resting stage in order to persist. The modelled potential distribution encompasses the known global distribution, save those locations where its presence is maintained through glasshouses (e.g., northern China, South Korea and northern Iraq). In these areas, the CLIMEX ecoclimatic index is zero, indicating that conditions are unsuitable for year-round persistence due to excessive cold stress during winter ( Supplementary  Information, fig. S1). However, the GI A is suitable for supporting significant population growth during the warmer months when field crops are grown ( Supplementary Information, fig.  S2). The ecoclimatic index suggests that substantial areas in the Americas, Africa and eastern Asia are likely to be at risk of further invasion by B. tabaci MEAM1.

Discussion
The known distribution of B. tabaci MEAM1 is extensive, confirming its reputation as an extremely important invasive species. Moreover, its modelled potential distribution is even greater, especially in South America, Africa and Asia ( fig. 3, bottom panel). The CLIMEX model performed very well, giving perfect model sensitivity in relation to all known qualified distribution records, and no excessive model prevalence. The areas modelled as climatically suitable, but without known distribution records all appear to be biologically plausible.
The modelled potential distribution of B. tabaci MEAM1 completely overlaps that of MED and the point distribution records for B. tabaci s.l., as well as the known country-level records globally. Further, the observed development rates of B. tabaci MEAM1 and Australia (AN) differ only slightly (De Barro and Hart, 2000). Such overlaps and similarities reflect the difficulties in distinguishing the different taxonomic entities within the B. tabaci complex (Muniz and Nombela, 2001;Qiu et al., 2009;Boykin, 2014). It is possible that the niche differentiation between the various species within the complex are less distinguished by climatic niche, and more by host preference, virus transmission patterns, etc. (Iida et al., 2009;De Barro et al., 2011;Tajebe et al., 2014). The fundamental climatic niche of the various species in the B. tabaci complex may, in fact, be extremely similar. Inter-specific competition between these morphologically indistinguishable conspecifics has been noted previously, with invaders displacing resident species (McKenzie et al., 2004;Liu et al., 2007;McKenzie et al., 2009;Legg et al., 2014). A detailed understanding of the relative competitive ability of each of the species in the complex under various environmental conditions is likely to remain elusive. Therefore, it is perhaps prudent to consider all of those areas identified as climatically suitable in the CLIMEX model as being at risk of invasion by B. tabaci MEAM1, at least until there is evidence that another member of the complex can restrict its spread through inter-specific competition.
There is growing concern that Cassava brown streak virus and Uganda cassava brown streak virus will move from East Africa to West Africa, especially to Nigeria, via whitefly introductions or the movement of infected cassava cuttings. Nigeria is the largest cassava producer in the world (FAOSTAT, 2015), and if these two viruses are introduced into the region, where cassava production is already hampered by Cassava Mosaic Disease virus, serious food security and economic implications will ensue. It is imperative to stop the spread of both whiteflies and infected planting material to the region.
The pest risk posed by B. tabaci MEAM1 has several components. The invasion risks posed under natural rainfall conditions are extended by the use of irrigation and the presence of glasshouses, which allow survival under otherwise lethally dry or cold conditions. The set of climate suitability maps, including the composite scenario, provides the basis for assessing emerging invasion risks if the use of glasshouses or irrigation were to be extended. The underlying model can also be used with confidence to project the changes that may be expected under future climate changes.
The potential threats of further invasion by B. tabaci MEAM1 suggest that there may be value in exploring the potential for deployment of Eretmocerus spp. and other agents for classical biological control into areas that are suitable for B. tabaci MEAM1, but presently remain uninvaded. Such pre-emptive regional assessments could speed up the response to new incursions, thereby reducing the production shocks in newly invaded areas. Doubtless, the new molecular diagnostic tools available for distinguishing the species within the B. tabaci complex will be critical in matching control agents to target species. It may also be economically justifiable to develop crop host genotypes that are resistant (Smith et al., 2018).