Current and future scenarios of suitability and expansion of cassava brown streak disease, Bemisia tabaci species complex, and cassava planting in Africa

Occurrence data collection

We obtained occurrence data in Africa from the Global Biodiversity Information Facility (GBIF), including 1,422 records for cassava (accessed 14 June 2021, http://dx.doi.org/10.15468/dl.98jqqy), and 1,169 records of Bemisia tabaci species (accessed 17 June 2021, http://dx.doi.org/10.15468/dl.755ck5). For quality assurance, we excluded records without geolocation and only kept unique records to avoid duplication. The occurrence records of CBSD were integrated from published literature (Alicai et al., 2019), GBIF, and country-wide surveys from the Cassava Diagnostic Project (CDP) in East Africa (accessed 21 November 2021 via AgShare.Today). The CDP involved seven countries, namely Tanzania, Kenya, Uganda, Rwanda, Mozambique, Malawi, and Zambia, following a published sampling protocol (Alicai et al., 2019). Through the Tanzania Agricultural Research Institute (TARI), we only accessed data from Tanzania and Uganda, which included CDP surveys from the National Crops Resources Research Institute (NaCRRI, Uganda) between 2004 and 2017. The data collection protocol are described in Alicai et al. (2019). In these surveys, cassava fields were randomly selected along motorable roads at intervals of 7–10 km and up to 20 km, depending on the density of cassava plants. A farmer was identified and asked for consent to survey the field at each location. Fields with crops between 3–6 months after planting were selected for the surveys, as CBSD foliar symptoms become apparent at this stage and before leaf shedding. Field location coordinates were collected using handheld GPS devices. In each field, 30 plants of the predominant variety were surveyed along two diagonal transects in an X-shape, with 15 representative plants each transect. Each sampled plant was scored for severity of foliar and stem symptoms on a 1–5 scale, where 1 indicates no visible symptoms and 5 corresponds to pronounced/extensive vein yellowing, chlorotic blotches on leaves or severe lesions, streaks on stems, or defoliation and die-back (Alicai et al., 2019).

To assess the impact of Bemisia tabaci species lumping, we divided whitefly occurrence records into datasets based on the geographic distribution of members of Bemisia tabaci species. The first dataset includes occurrence records from East Africa (Uganda, Tanzania, Kenya, Rwanda, Burundi, Malawi, Madagascar, Mozambique, and South Sudan). This dataset represents members of Sub-Saharan Africa 1 (SSA1) that are prevalent and widespread in this region (Mugerwa et al., 2012). The second dataset included occurrence records for Central and West Africa (Cameroon, Nigeria, Benin, Central Republic, DR Congo, Egypt Zambia and South Africa). This dataset represents members of Sub-Saharan Africa 2-5 (SSA2-5) and other Bemisia tabaci species that are widespread in this region. Then, the two datasets were compared with the full dataset from the African continent.

Environmental variables

The extent of our current study is Africa, the second-largest continent in both area and human population, stretching 30.4 million km² in land area with six distinct climate zones: the equatorial, humid tropical, tropical, semi-desert (Sahalian), Mediterranean, and desert (Beck et al., 2018). North Africa has an arid desert climate characterised by high temperatures and little precipitation. Equatorial West and Central Africa have a monsoon climate characterised by high temperatures, soaring humidity, and heavy seasonal rains (Chemura, Schauberger & Gornott, 2020). East Africa is characterised by dry and rainy seasons, while the southern part of Africa is generally more temperate (Nicholson, 2017). For the current climate (Hijmans et al., 2005), we used the 19 interpolated bioclimatic variables of 1970–2000 from the WorldClim database (https://worldclim.org/) at 10-min (∼340 km²) resolution (Table 1). We also included other biologically relevant predictors, including cassava harvested area (CHA), elevation (elev) and seven measures of soil quality (sq1-7) (Table 1). The CHA is an essential determinant of disease presence, and we used a standardised, high-quality, representative cassava map from Szyniszewska (2020), accessible via Figshare repository (10.6084/m9.figshare.9745118). Elevation can influence the occurrence and dispersal of Bemisia tabaci species by altering precipitation, temperature, vegetation, including crops, and the angle, direction, and intensity of the solar radiation (Mudereri et al., 2021; Ramos et al., 2018). Elevation data was downloaded from the digital elevation model (DEM) of the shuttle radar topographic mission (https://srtm.csi.cgiar.org/), available at approximately 90 m pixel size with a vertical error of less than 16 m. Although cassava is a resilient crop, adaptable to diverse and poor soils, the quality of soils substantially improves the crop’s productivity. We used seven key soil qualities that are important for crop production: nutrient availability, nutrient retention capacity, rooting conditions, oxygen availability to roots, excess salts, toxicities, and workability, denoted as sq1 to sq7, respectively, from the FAO soils portal (https://www.fao.org/soils-portal) at a 30 arc-sec (about 1 km²) resolution (FAO, 2000).

To assess the potential future distributions of cassava, its pests Bemisia tabaci species, and the CBSD outbreaks, we used the “BCC-CSM2-MR” climate model developed by the Beijing Climate Center of the China Meteorological Administration (Wu et al., 2019). The model was chosen based on a comprehensive assessment of the predictive capacity of eight GCMs (Figs. AC1-8, AD1-8, & AE1-8) on the three target species (cassava, whitefly, and CBSD). Our assessment revealed no marked differences in their ability to predict the future distribution of the target species. Although this model is part of 49 climate models included in the 2021 IPCC sixth assessment report (AR6), only eight models with complete coverage were available from the WorldClim 2.1 database (https://worldclim.org/) (Table 2). We used predictions from this model under two SSPs (SSP1-2.6 and SSP5-8.5) for year 2050 (mid-century average for 2041-2060) and year 2070 (near late twenty-first-century average for 2061-2080). The SSP1-2.6 scenario is part of the “sustainability” SSP1 socio-economic family, representing the best case where the best policies are implemented. Conversely, the SSP5−8.5 scenario represents the worst-case scenario with high fossil fuel consumption throughout the 21st century and without climate mitigation policies (Meinshausen et al., 2020), resulting in global warming ranging from a low of 3.1 °C to a high of 5.1 °C by 2100.

Species distribution modelling

Before building the SDM, out of 28 predictors we selected a subset according to the variance inflation factor (VIF) to avoid high level of collinearity that increases the uncertainty in model parameters and decreases the efficiency and power of model predictions (Naimi & Araújo, 2016; De Marco & Nóbrega, 2018). We calculated the VIFs of all predictors using the “usdm” R-package and excluded the one with the greatest VIF sequentially till all remaining predictors have VIFs <10 (Naimi & Araújo, 2016). As a result, 15 predictors were selected for cassava, 9 predictors for whitefly and 10 predictors for CBSD.

We developed the SDM using the “sdm” package (Naimi & Araújo, 2016) in the R platform (R version 4.0.5; R Core Team, 2021). The sdm package provides an object-oriented, reproducible, and extensible platform, capable of handling an ensemble of models, and we evaluated a total of twelve algorithms for their ability to predict the potential distribution of the three species of interest (Table 3). The models are classified based on the nature of data used i.e., presence-only vs presence-background data (Elith et al., 2006). As the species occurrence data consists of presence-only records, we randomly added 1,000 pseudo-absences throughout the study area (Thuiller, Georges & Engler, 2013). Variable importance is a metric used to determine the contribution of predictor variables in explaining the species distribution (Naimi & Araújo, 2016). In the “sdm” package, this is handled using “getVarImp” function. To assess model performance, we used the conventional metrics of the area under the receiver operator characteristic (ROC) curve (AUC), which measures the discriminatory ability of each model (Jiménez-Valverde, 2012 and Nahm, 2022). AUC is a threshold-independent performance measure that reflects the probability of a randomly chosen presence site ranking above a background site. The AUC value is usually divided into five levels (Guan et al., 2021): 0–0.6 (fail), 0.6–0.7 (poor), 0.7–0.8 (fair), 0.8–0.9 (good), and 0.9–1 (excellent). In this study, the models with AUC outputs in the range of 0.97 <AUC<1 were selected for further analysis. A fivefold cross-validation technique with five repetitions was used to assess the model performance (Thuiller, Georges & Engler, 2013). Specifically, models were calibrated on a random sample of 80% of the occurrence data and evaluated on the remaining 20% (Ramos et al., 2018; Ramos et al., 2019).

Model selection and performance

Of the 28 predictors, different sets were selected after eliminating the effect of multicollinearity. For cassava, 15 predictors (Bio02–04, 08, 13–15, 18–19; Sq1–5, 7) were selected for modelling potential distributions. Nine predictors were selected for modelling whitefly distributions (Bio02–03, 12–14, 18–19; CHA, and elev), while 10 were selected for modelling CBSD (Bio02–03; 08, 11–14, 18–19; and CHA). Of the 12 SDMs evaluated (Table 3), nine of the models had an AUC score above 90%, however, only four models (brt, maxent, gam, and mars) performed the best (AUC>97%) and were selected for further analyses. These four top models produced, on average, an AUC of 98.0% for cassava, 98.3% for Bemisia tabaci species and 99.5% for CBSD. The performance of the ensemble models on two groups of species (SSA1 and SSA2-5) is provided in Appendix F. We excluded the other models from the ensemble partly because we considered them either overfitting (AUCs closest to 100%, including random forest and Mahal.dismo) or their AUCs less than 97% (including, GLMNET and Bioclim). The average of these four top models was used for mapping the current and future suitability habitats of targeted species.

Predictors explaining distributions of cassava, whiteflies, and CBSD

For cassava, isothermality (Bio03, relative importance: 31.6%), temperature seasonality (Bio04, 20.3%), precipitation of wettest month (Bio13, 13.0%), and mean temperature of wettest quarter (Bio08, 6.5%) were the top predictors explaining the observed distribution (Table 1). For whiteflies, elevation (23.6%), precipitation of the driest month (Bio14, 13.6%), precipitation of the wettest month (Bio13, 13.5%) and precipitation of the coldest quarter (Bio19, 9.3%) strongly impacted the current distribution (Table 1). In contrast, cassava harvested area (CHA, 14.6%), precipitation of driest month (Bio14, 9.6%), annual precipitation (Bio12, 9.7%) and mean diurnal range (Bio02, 6.9%) were largely responsible for the distribution of CBSD outbreaks (Table 1).

The probability of occurrence responded to different environmental predictors in a nonlinear way (Fig. 1). Specifically, for cassava (Fig. 1A), the occurrence probability tends to increase with isothermality (Bio03), and precipitation of the wettest month (Bio13), but decline with mean temperature of wettest quarter (Bio08). For whiteflies (Fig. 1B), the occurrence probability tends to increase with precipitation of wettest month (Bio13) and precipitation of driest month (Bio14), while decline with cassava harvested area (CHA). In contrast, cassava harvested area (CHA) and precipitation of driest month (Bio14) showed a positive relationship with the occurrence probability of CBSD (Fig. 1C), while the occurrence probability showed a negative relationship with annual precipitation (Bio12) and mean diurnal range (Bio02).

Current and potential future suitable habitats for cassava in Africa

Cross-validation showed that the predicted range of suitable habitats based on 80% training records is consistent with the observed range from 20% testing records of cassava (Fig. 2A), with 54.6% (16.2 million km²) of the continent suitable for cassava production. These suitable habitats (i.e., suitability above 0.2) were predicted to locate predominantly in Sub-Saharan Africa. The regions with moderate suitability (i.e., 0.4–0.6), corresponding to 10.7 million km², were also found mainly in Sub-Saharan Africa. Although our data do not contain records from countries such as South Sudan, Sudan, Somalia, Botswana, and Zimbabwe, suitable habitats were predicted in these places. For habitats that are classified as ‘very suitable’ (suitability above 0.8), a total area of approximately 2.1 million km² was identified, predominantly in Coastal Guinea, Sierra Leone, Ivory Coast, Ghana, Togo, Benin, central to southern Nigeria, Cameroon, south-west Central African Republic, southern-western parts of South Sudan, Rwanda and the northern parts of Zambia (Figs. 2A–2B). Under the two SSP scenarios (SSP1.2-6 and SSP5.8-5) a potential increase of 59.6% and 56.6% in suitable habitats for cassava was projected for the mid-term (2050), predominantly in Sub-Saharan Africa (Figs. 2C–2D & Table 3). This trend was found to decline slightly towards 2070 (an increase of 57.8% and 55.3% compared to the current distribution). Specifically, highly suitable area (suitability above 0.8) for cassava production is expected to expand, from the East Coast of South Africa to Mozambique, and northern Madagascar (Appendix C; Fig. C1-8; for contracted areas, see Table 4).

Current and potential future regions at risk of whitefly invasion

Approximately 37.3% of Africa’s land area is predicted suitable for whiteflies, with a total of approximately 3.6 million km² highly suitable (>0.8), covering Uganda, southern Ethiopia, Rwanda, some parts of Tanzania, the Northern province of Zambia, southern DR Congo, central Cameroon, and some parts of Nigeria (Figs. 3A & 3B). For both SSP scenarios (SSP1-2.6 and SSP5-8.5), the ensemble model from the four selected SDM algorithms predicted an increase of highly suitable areas by 44.1% and 44.8%, respectively, by 2050 (Figs. 3C–3D), although the expansion was predicted to decline towards 2070 (an increase by 43.4% and 40.9% compared to the current distribution) (Appendix D; Fig. D1-8). The region with the highest future potential of whitefly infestation includes primarily East African countries, specifically South Sudan, Ethiopia, Zambia, Malawi, Angola, Central Cameroon, south-west Nigeria, Ghana, and Ivory Coast (cf. Table 5 for a summary of current and potential future suitable regions in Africa).

Current and potential future hotspots for cassava brown streak disease

The area currently susceptible to CBSD covers approximately 33.7% of Africa’s land area (10.2 million km²; Fig. 4A). The countries most vulnerable to this disease include the east coast and Lake zones of Tanzania, Uganda, and Southeast DRC (Table 5). Although no occurrence records were available in western Africa, the ensemble models highlighted suitable conditions in this region, specifically in Ivory Coast, Ghana, Nigeria (highest), and Cameroon, and East Africa emerged as a hotspot for CBSD outbreaks (Appendix E; Fig. E1-8). Under two SSP scenarios (SSP1-2.6 and 5-8.5) CBSD is expected to expand its range to 55% (16.6 million km²) and 56.6% (10.2 million km²) of Africa’s land area, respectively, by 2050 (Figs. 4B–4C), and this predicted range remains largely unchanged towards 2070 (Appendix E; Fig. E1-8), with countries most susceptible including those where the disease is yet to be reported, including Ivory Coast, Ghana, Benin, Nigeria, and Cameroon, all located in West Africa which is key for current cassava production. Other suitable areas were also found in southern Africa, specifically southwestern DRC, northern Angola, northern Zambia, and the eastern shores of Mozambique (Table 5).

Impact of Bemisia tabaci species lumping in predicting the current and future suitable habitats

Sub-Saharan Africa 1 (SSA1) Bemisia tabaci species are currently distributed and widespread in East Africa, and are likely to remain restricted to this region in future (Appendix F; Figs. F1A & F1D). On the other hand, the other Bemisia tabaci species, Sub-Saharan Africa 2-5 (SSA2-5) are currently distributed across Central and some countries in West Africa. With increasing warming in future, the suitable habitats of SSA2-5 species are likely to shrink to Ivory Coast, Ghana, Togo and Benin (Appendix F; Figs. F1B & F1E). Compared with maps produced using occurrence records from two geographical regions, the SDM model built using full Bemisia tabaci species records yielded consistent predictions for the current and future scenarios (Appendix F; Figs. F1C & F1F).

Factors of habitat suitability

We further found that different predictors are differentially driving the distribution of all three species of interest under current conditions. Only Bio 14 (precipitation of the driest month) was, to some extent, driving the distribution of all three species. This was more important for whitefly than for cassava. For cassava, we found temperature fluctuations and seasonality to have the strongest influence. Being a drought-tolerant plant, cassava can survive a wider range of environmental conditions, hence fluctuations, than insects (El-Sharkawy, 2004). The reproductive rate of whitefly is also said to decline with extreme temperatures (Aregbesola et al., 2020). As the carrier of the virus, the implication is that virus incidences will likely be reduced under unfavourable climate conditions, especially high temperatures (El-Sharkawy, 2004). Such high temperatures and precipitation events reported in Western Africa in this current term (Almazroui et al., 2020; Campo, Hyman & Bellotti, 2011) can explain the low prevalence of the vector and virus in this and other regions depicted by the current distribution.

As found in our study and consistent with the literature, elevation had the strongest influence on the distribution of whitefly but not cassava and CBSD. Jeremiah et al. (2015) found that whiteflies were more abundant in mid-altitude areas, i.e., 1,000–1,500 masl. However, areas with higher altitudes (>1,500 masl; Fig. 1B) are less suitable for whitefly, given their low temperatures (MacFadyen et al., 2018). Nevertheless, the presence of cassava fields will, to an extent, influence the occurrence of whitefly and spread of the disease (Donnelly & Gilligan, 2020). This was evident in our study by the contribution of cassava harvested area to the current distribution of both species. This contribution was at a lower magnitude for whitefly but was the highest driver of CBSD. It again amplifies the fact that under favourable environmental conditions, other aspects like infested cassava planting materials can continue to maintain the disease even in the absence of whitefly if the diseased cuttings are not adequately disposed of McQuaid, Gilligan & Van den Bosch (2017). This also explains the high vulnerability to the disease in East Africa.

Future scenarios of cassava planting

With all models and climate scenarios, Sub-Saharan Africa especially will become increasingly suitable for cultivation of cassava. This includes even small parts east of southern Africa. This is positive news in view of the importance and extensive use of cassava as a food security crop and for biofuel production (Nweke, 2004). This observation is consistent with other studies, which reported that suitable habitats for Cassava would increase across Africa (Jarvis et al., 2012). However, the above authors also reported that cassava suitability habitats would decrease in Central and West Africa. This contrasts with our findings which depict expansions across a greater proportion of Saharan Africa. A few reasons can be attributed to this discrepancy. Jarvis et al. (2012) used future climate information derived from an older phase of the CMIP project (CMIP 3, SRES-A1B emissions scenario). According to the climate modelling community, the CMIP3/CMIP5 over- or under projects climate changes in certain areas (Brown et al., 2013). Also, it was observed that while CMIP5 models did not capture the observed increase in global mean surface temperature between 1998 and 2013, the historical CMIP6 simulations are able to reproduce these observed recent temperature increases (Bock et al., 2020). More generally, CMIP6 has a higher ability to reproduce observed large-scale mean surface temperature and precipitation patterns compared to their CMIP3 and CMIP5 counterparts (Bock et al., 2020). Additionally, CMIP6 factors social-economic and policy aspects (Tebaldi et al., 2021; Zelinka et al., 2020). These authors also only used an EcoCrop model, while we used a more robust ensemble of four models.

The increased cassava suitability habitats reported in this paper also highlights Cassava’s ability to withstand climate change relative to other crops like maize, millet, sorghum, banana, and beans (Jarvis et al., 2012; Kotir, 2011). Nevertheless, this does not imply that climate mitigation actions should be ignored in favour of increasing the suitable habitats for Cassava, as our results do not speak to how final agronomic outputs will be affected. Crops are not only affected by temperature, but also generally exhibit different responses to increased CO₂ concentration. For Cassava, it is postulated that temperature increases could negatively affect biomass and yields due to poor assimilation (Amuji, 2021). Climate reports indicate that under these shared socio-economic pathway scenarios, temperature levels will already be well above 1.5 °C and 2 °C above pre-industrial periods (Osima et al., 2018). These increases are firmly cautioned against in the Paris Agreement on climate change (UNFCCC, 2016).

Increases in suitable habitats for cassava also correlate positively with slight increases in whitefly density and slighter increases in CBSD density, as seen in outputs from our future projections. This supports the likelihood of westward spread of the disease, especially along the coastal parts of Western Africa (Bundi et al., 2022). However, the tropical rainforest and monsoon conditions found in the DRC, especially central DRC, and zones along the Sahel belt may continue to have a positive impact in reducing the vulnerability of the region to both whitefly and hence CBSD. This is consistent with what was postulated by Campo, Hyman & Bellotti (2011)—that whitefly distribution is limited by extreme dryness along the Sahel belt and extreme and prolonged wetness in parts of central Africa. More so, even suitable areas for Cassava would be much reduced in central DRC and along the Sahel hence further limiting disease outbreaks. Adequate pest mitigation measures at international borders, coupled with good farm management strategies, would therefore largely favour a reduced vulnerability to whitefly and CBSD in these regions. The conclusions from this study corroborate a modelling study (using a different approach) looking at historical climate change in the East Africa region (Kriticos et al., 2020). The maps for the current distribution of Bemisia tabaci species show large areas suitable for whitefly species but without GBIF records. One of the reasons is attributed to the robustness of our model prediction. The ensemble modelling approach used could extrapolate suitable habitats for whiteflies beyond their current geographical range. This was also true for the prediction of cassava brown streak disease. Indeed, CBSD can be introduced into new areas via planting materials. Using cassava harvested area (CHA) as a covariate we found suitable areas in West Africa where the disease is currently absent (Appendix B; Fig. B4). This agrees with other modelling studies evaluating the role of mixed modes of infection in the spread of CBSD. They found that infected planting material through trade was the key long-distance pathogen dispersal mechanism (McQuaid, Gilligan & Van den Bosch, 2017).

The current distribution of the entire African species complex seems to be a composite of the current distributions of SSA1 and SSA2-5. However, there is a possibility that the southern strip in the distribution from Zambia through Angola to the coast of Namibia for the species complex has been overestimated. This overestimation might be attributed to the southern records of SSA2-5 being geographically connected with those of SSA1, resulting in an inflated southeast strip of high suitability in the distribution of the species complex.

It is worth noting that many SDMs target taxonomic units above the species level, and statistically, there is no limitation on using SDMs beyond the species level (Stas et al., 2020). While this approach is commonly used when dealing with a functional guild, it has been questioned when dealing with cryptic species. Overall, lumping species of Bemisia tabaci could lead to an overestimation of habitat suitability in areas with records of mixed lineages. This overestimation could be problematic as these mixed lineages may respond to climate change independently. Therefore, there is a need for more precise species identification in these areas.

One of the limitations of this study is the lack of species-specific data for whitefly, Bemisia tabaci, and cassava brown streak virus on GBIF. This led to lumping all Bemisia tabaci species and species causing CBSD together. However, Bemisia tabaci is now recognized as a pest species complex of over 25 species (Tay et al., 2022). In East Africa, two distinct species of Bemisia tabaci have been identified, including sub-Saharan Africa 1 (SSA1) and one from South West Indian Ocean (SWIO) islands (Mugerwa et al., 2012). The availability of these data in GBIF would improve the prediction of species distribution models. Given the generic approach used in this study (Campo, Hyman & Bellotti, 2011; Mudereri et al., 2021; Ramos et al., 2018), we believe that the resulting maps provide a good approximation of the true distribution of the species complex, but can be refined in the future for species-specific mapping.

Although Africa is one of the largest cassava producers, occurrence records in online repositories do not reflect what is known about cassava production (Szyniszewska et al., 2021). It is reported that less than 4% of the total records published on GBIF pertain to Africa. To bridge the gap of data deficit, efforts should be increased to encourage data capturing and sharing in order to promote more evidence-based studies. Future studies should also focus on how well these cultivars can perform under these climate change conditions under controlled greenhouse conditions to understand if the increased suitability habitats will be commensurate with yield increases.