Improving climate suitability for Bemisia tabaci in East Africa is correlated with increased prevalence of whiteflies and cassava diseases

Projected climate changes are thought to promote emerging infectious diseases, though to date, evidence linking climate changes and such diseases in plants has not been available. Cassava is perhaps the most important crop in Africa for smallholder farmers. Since the late 1990’s there have been reports from East and Central Africa of pandemics of begomoviruses in cassava linked to high abundances of whitefly species within the Bemisia tabaci complex. We used CLIMEX, a process-oriented climatic niche model, to explore if this pandemic was linked to recent historical climatic changes. The climatic niche model was corroborated with independent observed field abundance of B. tabaci in Uganda over a 13-year time-series, and with the probability of occurrence of B. tabaci over 2 years across the African study area. Throughout a 39-year climate time-series spanning the period during which the pandemics emerged, the modelled climatic conditions for B. tabaci improved significantly in the areas where the pandemics had been reported and were constant or decreased elsewhere. This is the first reported case where observed historical climate changes have been attributed to the increase in abundance of an insect pest, contributing to a crop disease pandemic.

www.nature.com/scientificreports/ niche of these cryptic species may have a large degree of overlap, and their realised niches reflect differences in biotic factors such as the distribution of hosts and natural enemies, interspecific competition and interactions with plant diseases [66][67][68] . The range boundaries for a given taxa tend to reflect spatial demographic processes, and competitive advantage can be conferred by even slight differences in ecophysiological performance 69,70 .
In this paper, we apply the 'end-to-end' method of Rosenzweig et al. 18 to explore the role of historical climatic changes in the cassava disease epidemic in parts of sub-Saharan Africa. To confirm that the CLIMEX model fitted to B. tabaci MEAM1 is relevant to the SSA B. tabaci we compare the temperature development response patterns for species of B. tabaci found associated with the cassava disease outbreaks in sub-Saharan Africa with those observed for MEAM1 and MED. We run the B. tabaci MEAM1 CLIMEX model of Kriticos et al. 64 with the CLIMEX Compare Locations/Years module, using a set of gridded monthly climate time series to test the hypothesis that the modelled annual climate suitability for B. tabaci MEAM1 is correlated with a 13-year time series of abundance of B. tabaci in Uganda 71 , and that the correlation is stable through time. As a further check of model relevance, we compare the CLIMEX model outputs for long-term average climate with data on B. tabaci abundance collected from cassava fields across Uganda, Tanzania and Malawi 72 . Finally, we use the time series of climate suitability to test whether the climate suitability for B. tabaci in East and Central Africa has been increasing when and where B. tabaci SSA, CMD and CBSD have been observed to increase in prevalence 28,30 . In so doing, we explore whether we can attribute the observed epidemiological and ecological changes to changes in modelled climatic suitability for B. tabaci SSA species.

CLIMEX.
The CLIMEX Compare Locations model simulates a species population response to climate on a weekly timescale; tracking the potential for population growth during favourable seasons, and population decreases during the stressful seasons. CLIMEX calculates an annual Growth Index (GI A ) to describe the potential for population growth as a function of weekly soil moisture and temperature during favourable conditions (Eq. 1).
where GI W is the weekly Growth Index, composed of the weekly Temperature Index multiplied by the weekly Moisture Index. The Temperature and Moisture Indices are specified with a functional form that accords with the Law of Tolerance [73][74][75] , and are combined in a multiplicative manner that accords with the Law of the Minimum 76 . CLIMEX employs up to four stress indices (cold, wet, hot, dry) and their interactions (cold-wet, cold-dry, hot-wet and hot-dry); respectively CS, WS, HS, DS, CWX, CDX, HWX, HDX) to estimate the ability of the population to survive unfavourable conditions (Eqs. 2, 3).
(1)  64 ). Outlying location records in northern China are thought to represent populations that overwinter in glasshouses in areas that are otherwise unsuitably cold for population persistence of B. tabaci MEAM1. The climate data is the CM30 1995H V2 dataset 65 . Map produced using ArcMap 10.6 (ESRI, Redlands, Ca., esri.com). www.nature.com/scientificreports/ The interaction stresses are used infrequently, and often only one of these is employed in any single model. The Ecoclimatic Index (EI) summarises the balance between the opportunity for the species to grow during the favourable season(s) and the requirement to survive inclement season(s), and scales from 0 (unsuitable) to 100 (optimal) (Eq. 4).
Very large values of EI are unusual and are restricted to environments that are climatically highly stable, e.g. the wet tropics.
HWX 100  64 ). The climate data is the CM30 1995H V2 dataset 65 . Pandemic fronts digitised from Legg, et al. 28 . Map produced using ArcMap 10.6 (ESRI, Redlands, Ca., esri.com). www.nature.com/scientificreports/ For the Compare Locations model a 30-year average of monthly climate data from 1981 to 2010 was derived from the time series dataset. We calculated the 30-year average from the abovementioned minimum and maximum temperature and precipitation. We then calculated relative humidity at 9:00 and 15:00 h by (1) calculating saturated vapour pressure at these times using the Magnus equation (Eqs. 4, 5 and 6, Kriticos et al. 65 ); (2) obtaining relative humidity by dividing vapour pressure by saturated vapour pressure (Eq. 1, Kriticos et al. 65 ); and (3) averaging 30-years by month. The minimum and maximum temperature, precipitation, and relative humidity were collated into an Sqlite database for use in CLIMEX Compare Locations.
The Compare Locations/Years model was run with a 39-year monthly time-series climate data provided by the University of East Anglia's Climatic Research Unit (CRU) data 77,78 . Climate data with a spatial resolution of 0.5° × 0.5° were extracted for monthly averages of daily maximum temperature, maximum temperature, vapour pressure and monthly totals for precipitation.
Comparing development rates of Sub-Saharan African Bemisia tabaci with MEAM1. In order to assess whether the CLIMEX model parameters fitted for B. tabaci MEAM1 were relevant for modelling the climate suitability for SSA B. tabaci, we assessed the development rates of three B. tabaci SSA taxa [identified as SSA1-SG1, SSA1-SG2 and SSA2 based on their partial mtCO1 sequences in line with Legg et al. 28 subgroupings] as a function of temperature under laboratory conditions. Until recently, the known geographical range of SSA taxa was scant due to limited molecular biological assays of B. tabaci in Sub-Saharan Africa 45,46 . Furthermore, their distribution is likely to be complicated by inter-specific competition from within the B. tabaci complex, and the importance of crop hosts as a likely range-discriminating factor between these species 43,79 . A consequence of this is that we are unable to reliably infer any stress parameters for these taxa based on their geographical distribution. On the other hand, the soil moisture parameters are strongly associated with the hosts, and factors such as the minimum soil moisture for growth are closely associated with the permanent wilting point, which does not vary much between species. The SSA1-SG1 and SSA1-SG2 B. tabaci colonies used in this study were established from B. tabaci collected from Kayingo, Uganda in February 2016. The SSA2 B. tabaci colony used in this study was established from B. tabaci collected from Kiboga, Uganda in August 2013.
Six temperature treatments: 15, 20, 25, 30, 35 and 40 °C which were constant for day and night were used in the experiment. Each treatment was replicated five times across five cassava plants of cultivar Ebwanateraka, each in a Lock and Lock (LL) container. In some cases, one or two replicates out of the five failed, and hence reliable data was only obtained for three or four replicates for such treatments. The cassava cultivar Ebwanateraka was used because it was the predominant variety grown by farmers before the cassava mosaic disease pandemic in the early 1990s 80 .
Ten pairs (10 female and 10 male whiteflies) were introduced onto a cassava plant in an LL container. The LL container with the introduced whiteflies was left at room temperature (25 °C) for 24 h in the insectary at 14:10 light and darkness, at a relative humidity of 60%. After 24 h of feeding and ovipositing, all the adult whiteflies were removed from the cassava plant in the LL container and the number of eggs laid was recorded. Cassava plants in the LL containers were maintained at room temperature in the insectary for another 9 days and on the tenth day the numbers of emerged nymphs on each leaf were recorded. After recording the number of emerged nymphs, plants in the LL containers were transferred to A1000 growth chambers (Conviron Europe limited, UK) set at a given temperature treatment, 14:10 light and darkness, relative humidity of 60% for 5 days. On day 16, cassava plants in the LL containers were removed from the A1000 growth chambers and maintained at room temperature 25 °C in the insectary. Cassava plants in the LL containers were monitored every 2 days for emerged adults. The number of emerged male and female progenies were recorded, and the emerged whiteflies removed from the plant.
Insect development rate was calculated as the inverse of the time taken for the insect to complete development. The mean time taken for whitefly adult emergence (3-5 replicates) was plotted against temperature for each whitefly taxa, and the response patterns were assessed.
Climate suitability for Bemisia tabaci MEAM1 through time. The CLIMEX Compare Locations/ Years model 53 was run on a monthly time series climate database (Climate Research Unit of the University of East Anglia) for 1978-2017 to simulate and visualise the seasonal and inter-annual spatio-temporal patterns of climate suitability for B. tabaci MEAM1 64 . We used the R statistics package 81 to fit linear time trend models at each spatial grid point to estimate the average annual change in EI, GI A , TI and MI across the 39 years in the series. At each point, standardised effect sizes (t value) for this annual change were plotted spatially, as were the 95% lower and upper confidence intervals for the annual change.
Comparing modelled climate suitability for Bemisia tabaci MEAM1 with field abundance of Bemisia tabaci in Uganda, 2004-2017. A comprehensive set of data on the field abundance of B. tabaci and Cassava Brown Streak Disease (CBSD) in Uganda was published recently 71 . This dataset includes field survey data collected from 2004 to 2017 with a gap for 2016, with 7 627 field summaries across 96 districts. While the survey effort changed in intensity and spatial pattern from year to year, it remains the most useful dataset of its kind. The field data show a clear pattern of invasion of CBSD into cassava fields in Uganda through this period, and these invasion dynamics confound any relationship between prevalence of CBSD and climate suitability for B. tabaci sensu lato.
To test the relationship between observed abundance of B. tabaci (field counts) and modelled climate suitability for B. tabaci MEAM1 throughout the same period as the dataset of Alicai, et al. 71  The geocoded locations for the cassava fields were spatially intersected with CLIMEX results for East Africa. A logistic model was fitted in R (using the glm function with a binomial link) to the probability of presence of one or more whiteflies on a cassava plant as a function of CLIMEX EI. We ran the CLIMEX model using a 30-year climate average dataset centred on 1995 using the CRU data described above. The intent in comparing the modelled climate suitability with abundance data was to discern the patterns of agreement, rather than to be able to use modelled climate suitability to predict abundance with any great precision.

Results
Development rates. The observed pattern of development rates for B. tabaci SSA taxa (Fig. 3) 82 . We can therefore feel confident in using the previously developed CLIMEX model to adequately simulate climate suitability for the SSA taxa. Globally, the potential distribution for B. tabaci MEAM1 extends throughout the tropics, sub-tropics, and warm temperate climates, with marginal suitability in the Mediterranean region (Fig. 1). Location records in temperate locations highlight its potential for recurrent invasion from glasshouses, which allow overwintering. The climate suitability for B. tabaci MEAM1 in eastern Africa modelled using long-term average climate is quite high in the north west (Uganda) and South-East (coastal Tanzania) of the study area, with a moderate depression in Central Tanzania (Fig. 2) 64 . The modelled climate suitability for B. tabaci MEAM1 accords with the distribution records for all of the sub-Saharan African species (sensu 45 ) within the B. tabaci species complex reported in Legg et al. 28 .
Field Abundance of Bemisia tabaci and modelled climate suitability Uganda 2004-2017. The relationship between observed field abundance of B. tabaci in Uganda 71 and modelled climate suitability for B. tabaci MEAM1 64 was statistically consistent throughout the 13 years of surveys, with the variation in slope between years appearing random, with 2009 a standout (Fig. 4).

Probability of presence of Bemisia tabaci in East Africa 2015-2016. The field abundance patterns
(probability of presence of B. tabaci on a cassava plant) displays a sigmoid pattern with respect to modelled EI at each site (Fig. 5). There were pronounced differences in prevalence of B. tabaci SSA adults between sites at moderate EI values, perhaps reflecting non-climatic influences on population abundance of B. tabaci such as crop type, age of the cassava, crop types used across each landscape, natural enemies, competition with other herbivores, and insecticide use. Figs. S2, S3). Despite this variability, there has been a significant trend in climate suitability throughout much of Eastern Africa during the 39 years from 1979 to 2017 (Fig. 6). For each cell there were 38 degrees of freedom. Throughout this period, there has been a slight decrease in climate suitability for B. tabaci in lower-lying coastal regions (Fig. 6a), though apart from an area in Mozambique, only a few isolated areas appear to have experienced a statistically significant decrease in suitability (Fig. 6c). Conversely, there has been an increase in suitability throughout a large portion of the area, which is most pronounced in the Democratic Republic of the Congo, Uganda, Rwanda, Burundi, and western parts of Tanzania and Kenya (Fig. 6a,b). These trends are apparent in both the EI and the GI A results against a backdrop of interannual variation (Fig. 7a,b). Exploring the temporal data in more detail for selected locations reveals that most of the increase in suitability is due to increasing temperature suitability (Fig. 7c) and decreases in the EI are due to decreasing soil moisture suitability and drought stress (Fig. 7d). At Luweero and Mwanza, in the centre of the cassava disease epidemic zone, climatic conditions for both temperature and soil moisture improved through this period.

Discussion
Both the field and laboratory studies support the application of the CLIMEX model originally fitted to development data and field distribution records for B. tabaci MEAM1 to examine the historical change in climate suitability for the B. tabaci SSA species. This lack of variation in the fitted model parameters between the different B. tabaci species is perhaps not so surprising. Members of the complex are often found to have sympatric ranges. Bemisia tabaci MED and MEAM1 frequently coexist 64,83 , and in East Africa the SSA taxa frequently co-occur 43 . It is likely that the non-climatic ecological factors such as host relations (e.g., 83 ) provide bases for niche differentiation within the B. tabaci species complex.
The Ugandan time-series data shows both a statistically consistent relationship between B. tabaci abundance and modelled climate suitability (Fig. 5), and an increasing suitability trend across time. The apparent increasing climate suitability for SSA B. tabaci in parts of East and Central Africa across the 39 years of the study (including the 13 years of the Ugandan dataset) accords with the reported increased prevalence of both B. tabaci 33 and cassava diseases caused by viruses vectored by B. tabaci 26,27 . This sustained increase in climate suitability may help explain the increasing abundance of SSA B. tabaci 27 , in combination with increased use of cassava varieties that support B. tabaci species, and changes to cropping systems that provide additional host plant resources 40 . An increase in suitability across time was most obvious around the lake zone of Southern Uganda and Northern Tanzania, and west into Rwanda. This matches a more than five-fold increase in observed abundance of B. tabaci in this area between 1994 and 2009 33 . In contrast, parts of Malawi have seen little change in the climate suitability during this historical time period. Some sites such as Lira in Northern Uganda and coastal Tanzania have seen a slight decrease in suitability for B. tabaci; between 1994 and 2009, Jeremiah, et al. 33 observed a three-fold decrease in abundance at coastal Tanzanian sites, which accords with the model results for Mbawala and Bagmoyo (Fig. 7).
In order to attribute these observed changes to anthropogenic climate change, the observed changes need to satisfy the following conditions: As observed in Rosenzweig et al. 18 , attributing climatic changes to natural systems requires special treatment through 'joint attribution' . It seems clear that the recent reports of increased abundance of whiteflies and cassava www.nature.com/scientificreports/ diseases are not due simply to natural variation (internal variability). Cassava has been grown in the region for a long time prior to the recent reports of disease pandemics, and there have been a range of native whiteflies present in the region 43 . By using the CLIMEX model of B. tabaci MEAM1 64 , we have satisfied the second and third criteria for climate change attribution. The CLIMEX model was fitted using long-term average historical climate data and experimental data on MEAM1 response to climate variables (cross-checked with SSA responses), thus satisfying the plausibility criterion. When we applied the time-series climate data (natural forcing), the resulting pattern of statistically significant changes in climate suitability matched the observed patterns of changes in whitefly abundance and cassava disease. The trend was clearly sustained throughout the 39 years of the study, and was apparent despite the substantial inter-annual variability in climate suitability; further supporting the notion that it was not due to internal variability in the system, nor due to shorter-term variability in climate (e.g. inter-decadal) 15 . The null hypothesis for climate change is in-effect the baseline CLIMEX model for B. tabaci MEAM1, applied to the historical climate dataset centred on 1995. By drawing upon a process-oriented model such as CLIMEX in this manner, we are essentially utilizing a variant of the 'end-to-end' method of climate change attribution 18 .
In attributing the B. tabaci abundance and cassava disease changes to climatic changes we also need to consider other potential explanations and confounding effects. As observed by McQuaid et al. 85 , the movement of diseased cassava material can play an important role in the spread of a disease, especially over longer distances. However, the subsequent maintenance of an epidemic relies upon either a local cycling of the disease or constant re-introduction of diseased plant material. While farmers may spread diseased material locally, this is likely a small effect. Conversely, there is compelling evidence that B. tabaci plays a strong role in spreading the viruses within the landscape and maintaining epidemics 85,86 , so we can be confident that while the spread of diseased cassava cuttings into the epidemic zones was a factor in the spread of the disease 71 , the observed epidemic was most likely driven by B. tabaci dynamics and abundance.
The statistically significant trend of improving climate suitability for B. tabaci in the cassava disease pandemic area in East and Central Africa, centred mainly on Uganda appears to be mostly driven by increasing temperatures (Fig. 7c). At some sites at the centre of the recent cassava disease epidemics (e.g. Luweero, Uganda and Mwanza, Tanzania), decreasing rainfall also contributed to the increasing favourability for B. tabaci via the increase in the Moisture Index (Fig. 7d). This accords with the general perception of B. tabaci having a preference for hot, dry conditions 87 , and references therein. The preference for dry conditions is at least partly due to the susceptibility of nymphs to being dislodged from leaves by raindrops, a fact that has led to overhead irrigation being used as a means of suppressing B. tabaci 82,88 .
Where there were trends toward decreased suitability (e.g., in sub-coastal Tanzania and northern Kenya, Fig. 6a,c) the decreases also appear to be associated with decreasing rainfall. However, these predominantly rangeland areas were already poorly suited for cropping due to inadequate rainfall and became even drier and less suitable for B. tabaci during the period of this study. In East Africa, the trends in increasing climatic suitability   www.nature.com/scientificreports/ arthropod 91,92 . This result then begs the question of how to use climate change scenarios to assess the potential extent of these viral epidemics, though this is beyond the scope of the present study. www.nature.com/scientificreports/ Given the apparent increasing climate suitability trend for B. tabaci in parts of East and Southern Africa, we might expect that there could be more generations of B. tabaci and greater survival rates, and hence abundance. African Cassava Mosaic Virus (ACMV) prevalence is related to the density of whiteflies 93 , so we might reasonably expect that the increasing climatic suitability for B. tabaci could translate into greater disease prevalence. While Jeremiah et al. 33 did not find a positive correlation between B. tabaci abundance and incidence of Cassava Brown Streak Disease (CBSD), based on epidemiological theory we might also expect that increasing B. tabaci density would have a proportionally greater effect on the levels of cassava diseases with relatively poor transmission rates (e.g. CBSD) compared with those with high transmission rates (e.g. CMD or ACMV) because the latter would saturate at lower density levels of B. tabaci. Logically, because transmission rates are observed on a perinsect basis, higher insect density results in higher plant infection rates, up to the point where the availability of uninfected hosts becomes limiting.
The sustained trend in climatic suitability for B. tabaci in parts of East and Central Africa suggests that a similarly sustained effort will be required to develop and maintain cassava varieties that are resistant to current and emerging strains of cassava diseases. Otherwise, the contribution of cassava production toward African and global food security will likely be compromised. Clearly, such plant breeding initiatives need to be paralleled with efforts to control B. tabaci using biological and cultural control methods. Methods such as those employed by Pardey et al. 94 to estimate the economically appropriate amount to invest in perpetuity in developing and maintaining resistance to wheat stem rust (Puccinia graminis) might be applied to the protection of cassava production from arboviruses. Instead of framing this problem as a battle to be won-a project with a start and end date-it can be framed as a perpetual investment stream to develop cassava varieties that are resistant to the contemporary challenges from the evolving patho-system.
The newly available Compare Locations/Years tool in CLIMEX Version 4, coupled with gridded time series climatic data allowed us to explore the effects of climate variability on this important agricultural pest system. It enabled us to look not only at the resultant trends, but also to identify which components of climate were driving these trends. Linking climate dynamics to a process-oriented niche model is potentially an important mechanism for exploring trends in global change biology. The process-oriented nature of CLIMEX, combined with time series climatic data, lends it to serving a useful role in climate change attribution in many biological systems.
Our modelling identified that East and Southern Africa includes areas that are highly climatically suitable for B. tabaci SSA species, and that this suitability has been increasing for many years in certain regions, highlighting an additional potential biosecurity and food security threat for the region. At present, the native B. tabaci SSA species already pose a significant threat to cassava production, and there are other related whitefly species also present in cassava fields. The role that native African whitefly species play in terms of transmission of a large diversity of cassava viruses is sometimes unclear. Africa's notoriously porous land borders 95,96 make it difficult to control the spread of invasive species and emerging infectious diseases. West African cassava production is greater than that in East Africa. If Cassava brown streak virus disease or Ugandan cassava brown streak virus were to be introduced to West Africa, for example via the movement of infected cassava cuttings, the threat to cassava production there could be even greater than that observed in East Africa if a suitable whitefly vector is either already present, or invades that area. Clarifying the ability of West African whiteflies to transmit CBSV and UCBSV, may help define a range of strategies to manage these disease risks to cassava production in this region.