Climate change favours connectivity between virus-bearing pest and rice cultivations in sub-Saharan Africa, depressing local economies

Target species and study area

Chaetocnema Stephens, 1831 is a flea beetle genus widespread in all zoogeographical regions. It consists of over 300 described species. About the Afrotropical region, Chaetocnema genus includes over 100 known species (Biondi, 2002; Biondi & D’Alessandro, 2006a, 2006b, 2018). The species of this genus are often associated with moist environments and plants in many botanical families, particularly the Chenopodiaceae, Polygonaceae, Cyperaceae, Gramineae (Biondi, Urbani & D’Alessandro, 2015). A small group of very similar species (Biondi & D’Alessandro, 2008), including Chaetocnema pulla Chapuis, 1879, C. nkolentangana Bechyné, 1955, C. subquadrata Jacoby, 1897, and C. vanschuybroecki Biondi & D’Alessandro, 2008, widespread in the Afrotropical region, are here considered. The occurrence localities of pulla group were mainly obtained from Biondi & D’Alessandro (2008) and integrated with new data (M. Biondi, 2021, personal data), obtaining the most complete dataset for this taxon to date.

About the target plant species, two main species of rice are cultivated in Africa. The first is Oryza sativa L., native to Asia and was introduced in Africa about 450 years ago. Another lesser-known species, O. glaberrima Steud, is native to Africa and was domesticated in the Niger River Delta (in northern Mali) over 3,000 years ago. As a result of their evolution and domestication, both species have distinct and (Carpenter, 1978) complementary advantages and disadvantages for African farming systems. Asian rice (O. sativa) is characterized by good yields, no entrapment and crushing of grain, and high fertilizer yields, unlike its African counterpart (O. glaberrima). However, unlike Asian rice types, native varieties of O. glaberrima often have good competitiveness and weed resilience against major African biotic and abiotic stresses (Koffi, 1980; Jones, Mande & Aluko, 1997). The study area is sub-Saharan Africa, where the two target taxa both occur.

Ecological niche modelling

To estimate the current and future climatic suitable area of both pulla group and rice, Ecological Niche Models were built; these have recently gained popularity among many fields of natural sciences, as they permit to better address conservation policies (Iannella et al., 2018; Brunetti et al., 2019), understand biogeographic patterns (Iannella, Cerasoli & Biondi, 2017; Console et al., 2020) and predict possible distribution of invasive species and pests (Fick & Hijmans, 2017).

For this aim, 19 bioclimatic variables for current climatic conditions were downloaded from the Worldclim.org repository, ver 2.1 (Fick & Hijmans, 2017), at 2.5 arc-min resolution. To avoid possible correlation among predictors, which may lower model’s performance, both Variance Inflation Factor (VIF, threshold set = 10 following Guisan, Thuiller & Zimmermann (2017)) and Pearson’s r (|r| < 0.85, following Dormann (2007), Elith et al. (2006)) were assessed through the ‘vifstep’ and ‘vifcor’ functions of the ‘usdm’ R package (Naimi, 2015), and a subset of predictors was then used for models’ calibration.

Ecological Niche Models were built through the Maxent (Phillips, Anderson & Schapire, 2006) implementation in “SDMtoolbox” 2.4 version (Brown, Bennett & French, 2017), in ArcMap 10.0 (Esri, Redlands, CA, USA). This software permits to take advantage of the powerful machine-learning algorithm of Maxent and to integrate into the modelling process some details (not available on the “standard” version of Maxent), which increase models’ discrimination power by lowering spatial biases. In particular, the ‘Spatial jackknifing’ procedure implemented in this software permits to reduce spatial autocorrelation among data by splitting the study area (and the corresponding presence occurrences and background points) into training and test datasets, iteratively (Brown, Bennett & French, 2017). For the present study, five spatial jackknife groups were set; moreover, considering the great extent of the study area, a bias mask to account for and correct the latitudinal bias was generated through the ‘Bias File for Coordinate Data (BFCD) in MaxEnt’ algorithm (Brown, Bennett & French, 2017) and applied to the modelling process. Before being used as input data for the models, both pulla group and rice occurrence data were processed through the ‘Spatially Rarefy Occurence Data for SDMs’ algorithm (resolution = 10 km) so as to make comparable the resolution of occurrences with one of the predictors (Sillero & Barbosa, 2020). To obtain reliable estimates, 10 replicates for each model parameter class were calculated, and a Regularization Multiplier = 5 was set, for a total of 250 single models for each species. Also, spatially segregated groups were chosen, as this option permits the models’ training and evaluation in potentially non-analogous environments (Brown, Bennett & French, 2017).

Models trained on current climate were projected to future conditions using three different general circulation models (GCMs) to reduce uncertainties due to their differences (Stralberg et al., 2015). For this aim, 2030, 2050 and 2070 GCM projections of BCC-CSM2-MR (Wu et al., 2019), the IPSL-CM6A-LR (Boucher et al., 2020) and the MIROC6 (Tatebe et al., 2019) were chosen, for two Shared Socioeconomic Pathways (SSP), the SSP245 and the SSP370. These were chosen based on their link between the possible future trends of land-use (see below) and two Representative Concentration Pathways (RCPs) scenarios, with a more optimistic (RCP 4.5) and a “middle of the road” one (RCP 7.0). Moreover, the Multivariate Environmental Surface Similarity (MESS, following Elith & Leathwick, 2009), which measures the projected models’ degree of extrapolation (i.e., the divergence of environmental predictors used for models calibration with respect to the ones used for models projections), was taken into account. This information was incorporated into the MEDI, an index used to proportionally weigh future models’ projections based on their degree of extrapolation measured through the MESS (Iannella, Cerasoli & Biondi, 2017). Models’ discrimination power was assessed by both the Area Under the Curve of the Receiver Characteristic Operator obtained from the Maxent output (Phillips, Anderson & Schapire, 2006) and the continuous Boyce index (ranging from −1 = counter prediction to +1 = optimal prediction), which is the best choice for evaluating presence-only ENMs (Boyce et al., 2002; Hirzel et al., 2006). This was calculated through the ‘ecospat.boyce’ function of the ‘ecospat’ package (Di Cola et al., 2017) in R environment (R Core Team, 2016).

Land-use simulations and GIS analyses

The ENMs obtained were further analysed in a GIS environment to refine the predictions in more realistic distribution scenarios of both target taxa (Iannella et al., 2019). All analyses of this section were performed in QGIS 3.10.11 (QGIS.org, 2021).

We performed future land-use-land-cover (LULC) simulations in sub-Saharan Africa using the InVEST suite, a geospatial modelling framework tool that evaluates the impact of land-use change on ES (Sharp et al., 2018; De Simone et al., 2020).

Furthermore, we evaluated rice production under the influence of pulla group in different future scenarios (SSP2 and SSP3 narratives) within a time range of 40 years (2030–2050–2070).

The basic LULC we have chosen for all the analyses performed is the Global Land Cover 2010 (used as a current scenario), distributed by the Copernicus Climate Change Service (C3S), which classifies the earth’s surface into 22 classes with a spatial resolution of 0.002778° (~300 m).

In addition to the map described above, to identify the rice crops more precisely, we used the SPAM (Spatial Production Allocation Model) map 2010 (Yu et al., 2020), isolating the “rice” class from other crops. Using GIS techniques, it was fused with C3S’s “rainfed/irrigated cropland” class to create an additional land-use category (called Rice) representing potential rice crops of the current situation.

To obtain simulations of future LULC maps, the ‘Scenario Generator proximity-based’ (SGpb) of the InVest suite (Sharp et al., 2018) was used. This tool allows to create simulations of possible alternative futures using the data provided by the user.

Following the creation of the new Rice class, we translated management policies into spatial data using “guidelines” (“storylines”). To this end, the SSP was used. In particular, the SSP database (SSPdb) of Riahi et al. (2017) was used, which quantifies the mitigation policies. This database was created to model and compare possible future trajectories of different LULC and climatic scenarios (deriving from RCP + SSP) and quantify such changes.

Future LULC scenarios have been built to strengthen the link between the SSP narratives and the new RCP scenarios, to have more precise information about future projections.

The socio-economic scenarios chosen were SSP2 4.5 and SSP3 7.0, and the new Rice class previously created was selected for the future simulations. Given the unavailability of the SSP3 7.0 scenario in SSPdb, we used the baseline SSP3 data because Fujimori et al. (2017) indicated that the case with a forcing level of 7.0 W/m² roughly corresponds to the SSP3 baseline cases.

To extract the real values of the rice crops area for our study region, we used the FAOSTAT database (FAO, 2016). To obtain the future growth rate of rice crops area, we performed a linear regression analysis on the data obtained by FAOSTAT (1973 ÷ 2010). We then performed a ratio between the cultivated areas of the SSPdb and those of the current LULC to obtain the data for future simulations of the land-use scenarios.

We considered the Croplands variation in millions of hectares for two SSP narratives: SSP2 and SSP3. In SSP2 (midway), land-use change is incompletely regulated, i.e., tropical deforestation continues, albeit with slow decline rates; crop yield rates slowly increase over time (Fricko et al., 2017).

The second SSP narrative considered, SSP3 (regional rivalry), shows a resurgent nationalism, concerns about competitiveness and security, and pushes countries to focus more and more on national or regional issues. Land-use change is hardly regulated and crop yield rates decline sharply over time, mainly due to the minimal transfer of new agricultural technologies to developing countries (Fujimori et al., 2017).

The values for the correct setting of the SGpb tool were obtained by calculating the variation (2010 ÷ 2030; 2030 ÷ 2050; 2050 ÷ 2070) between the present and future rice areas.

To simulate the future LULC, we applied the two SSP projections to cultivated land in the LULC classes that are most likely to be modified (see below) closer to existing cultivated land. Three types of land cover were chosen: (a) “Focal ground cover” (SGpb tool converts classes from the edge of focal land cover zones) was chosen based on the most likely categories to increase agricultural areas. It is more likely to extend agriculture from existing crops than to create new ones. (b) “Convertible land cover” represents the types of land cover that can be converted, chosen based on the classes that have changed the most in the last 20 years, and obtained from the analysis of these LULC variations (ESA, 2017). (c) “Replacement land cover” was chosen because the agricultural land cover class was the most frequently converted compared to the other LULC categories (De Simone et al., 2020).

To identify the areas of co-occurrence for the correct MEDI models (pulla group and rice), we polygonized and intersected them with the simulations of future scenarios (LULC_2030_SSP2_4.5, LULC_2050_SSP2_4.5, LULC_2070_SSP2_4.5, LULC_2030_SSP3_7.0, LULC_2050_SSP3_7.0 and LULC_2070_SSP3_7.0).

From the spatial data obtained, we calculated the resulting rice crops area values to estimate the surface area most affected by the pest and to evaluate the co-occurrence between the two taxa within each sub-Saharan African country.

Landscape connectivity

To assess the potentiality of the study area in allowing the current and future colonization of rice cultivations from where pulla group currently occurs (i.e., measuring landscape connectivity for pest invasion), Circuitscape v.5 (Anantharaman et al., 2020) in Julia 1.5.3 programming language (Bezanson et al., 2017) was used.

Circuitscape applies both random walks and circuit theory to model the ecological connectivity in a target area (McRae et al., 2008; McRae, Shah & Edelman, 2016) and is particularly useful and robust, with respect to other connectivity software, for applications in which the species has no previous knowledge of the route to use for the colonization (McClure, Hansen & Inman, 2016), as in this case.

For this aim, the ENMs outputs were converted into friction maps (the inverse of suitability maps) through the ‘Raster calculator’ tool in ArcMap 10.0 and used as resistance layers (i.e., the resistance which a landscape has to the movement of target individuals/taxa). The ‘Advanced’ mode was used to enable the activation of independent sources (pulla group occurrence localities) and grounds (the destination patches, in our case, the current and future predicted rice cultivations in raster format). Considering the great geographic extent of the study area, the ‘solver = cholmod’ (a solver method applying a Cholesky decomposition) was used (Bezanson et al., 2017). The outputs obtained are raster maps depicting a current value, which were subsequently managed in ArcMap 10.0.

Pest aggression index (PAI) and provisioning ecosystem services

To go further into the understanding of the consequences of the connectivity which the landscape offers to the invasion of pulla group to rice cultivations, we developed an appropriate index, which we call “Pest Aggression Index” (PAI), calculated through the corridor edge (CE) and Mean Corridor Magnitude (MCM, whose rationale is reported in Fig. 1), the parameters of which are introduced and detailed below.

When referred to target countries (PAI_Country), to obtain a single value for a target area (in our case, a country, even though the PAI can be applied to any territory), we first process the connectivity maps inferred in the previous steps. Considering that these are continuous rasters spanning the whole study area (Fig. 1A), we remove zero values to isolate each corridor from another (Fig. 1B). Then we calculate the mean cell raster value in R (R Core Team, 2016), obtaining a single number that represents the magnitude of an ideal vector which symbolizes the corridor, a value we name Mean Corridor Magnitude (MCM) (Fig. 1C). The MCM is then used to set the threshold of the continuous corridor rasters. Concurrently, we consider the length (in km) of the corridor edge (CE), which is in contact with the destination rice patch, representing the “frontline” which the pests will first find, moving through the landscape corridor (Fig. 1C).

To standardize the results and compare the PAI for each country, we refer to the multiplication of these first two terms to the total rice cultivation (RP) and agricultural area (AR) of the target country through the formula:

${\mathrm{P}\mathrm{A}\mathrm{I}}_{Country}={\displaystyle \frac{{\sum }_{i=1}^{n}C{E}_i\ast \mathrm{M}\mathrm{C}{\mathrm{M}}_i}{TotalR{P}_{Country}\ast TotalA{R}_{Country}}}$

The PAI_Country is divided by the maximum value of each year_SSP scenario and further normalized by subtracting the obtained value from one. Then, the value obtained is multiplied by the rice production of each country, obtaining a “weighted” one, which mirrors the potentiality that the pulla group populations have to reach and affect cultivations.

This index ranges from zero (no landscape connectivity among the pest’s populations and target cultivations) to one (the highest PAI, the highest expected losses to cultivations in each country).

Considering the data achieved from the previous analyses, we calculated the future potential rice production (FPRP) in the different scenarios (SSP2 4.5 2030–2050–2070 and SSP3 7.0 2030–2050–2070) for each sub-Saharan country. We calculated the FPRP by multiplying the areas of the rice patches (in co-occurrence with the pulla group) with the average rice yield (ton/ha). As previously done to obtain the growth rate of rice crops area, to calculate the rice yield’s future growth rate, we performed a linear regression analysis on the data obtained by FAOSTAT (1973 ÷ 2010).

To obtain the final values of loss in rice production (provisioning ecosystem service), we calculated for each country the difference between the total FPRP (without pulla) and the FPRP values obtained under the influence of pulla. Furthermore, to standardize these values, we multiplied the final value by the PAI values (by country and LULC scenario). We have condensed the loss in rice production (%) obtained by the countries in a heatmap achieved through the Seaborn library in the Python environment (Waskom, 2021).

To achieve a more accurate view of the economic losses by country, we also calculated the ratio between the economic loss and the gross domestic product (GDP) for each nation. We used the obtained economic loss data and the projected GDP values to the future taken from the SSPdb, Faostat and WFP VAM (World Food Programme (WFP), 2020) database. A selection of countries was made both for the results obtained and for the availability of the various data sources (rice price, GDP).

Models calibration and evaluation

After the variables’ selection process, nine bioclimatic variables were selected for the models’ calibration, namely BIO1 (mean annual temperature), BIO2 (mean diurnal range), BIO3 (isothermality), BIO8 (mean temperature of wettest quarter), BIO13 (precipitation of the wettest month), BIO14 (precipitation of driest month), BIO15 (precipitation seasonality), BIO18 (precipitation of warmest quarter) and BIO19 (precipitation of coldest quarter). The models built on this set of variables resulted in high discrimination power, with an AUC = 0.817 for Oryza and AUC = 0.807 for pulla group (Fig. S1A) and the Continuous Boyce Index obtained for Oryza = 0.983 and for pulla group = 0.957 (Fig. S1B). Moreover, low spatial standard deviation resulted for both the models calibrated on the current climatic conditions (Figs. S1C, S1D).

Current distribution and future co-occurrences

The occurrence localities of pulla group are clustered in specific areas, mainly Central-Eastern and Western Africa, as well as Madagascar, while rice is widely cultivated (obtained from GBIF (2020); accessed on November 18, 2020), even in territories where pulla group does not occur (Fig. 2).

These distributions are partially mirrored in the areas obtained through the ENMs. The combinations of the low, medium, and high suitability classes report broad areas in which cultivations in both Western and some territories of Eastern Africa may not currently suffer the presence of pulla group (Fig. 3).

The general trend of vast areas hosting the conditions mentioned above disappears in future projections (Figs. 4A–4F). Instead, an increase in the Mid-Mid, Mid-High and High-Mid suitability classes (intended as pulla-rice) is observed. For instance, the SSP_2.45 shows a +68% of the Mid-High class in the 2030 scenario (Fig. 4A), with some peaks observed in the 2070 SSP_3.70, where the Mid-Mid class reaches an increase of +84% (Fig. 4F).

Connectivity of pests to rice cultivations and provisioning ecosystem services

The landscape connectivity inferred from future scenarios show corridors allowing pulla group to reach rice cultivations in specific African territories, such as Nigeria, Cameroon, Sudan, South Sudan, Democratic Republic of Congo, Zimbabwe, Republic of South Africa, and Madagascar (Fig. A.2).

The Pest Aggression Index (PAI) values ranged from 0 to 0.992; these values were used to obtain the rice production loss (in percent) for each SSA country, as summarized in the heatmaps reported in Figs. 5A, 5B. Different results are observed according to the SSP scenario (2.45 and 3.70). The production loss for the SSP 2.45 and SSP 3.70 narratives (Figs. 5A, 5B) shows that the countries that will suffer the largest production losses (>−90%) in all the scenarios analyzed (2030, 2050, 2070) are Liberia and Sierra Leone. In the SSP 2.45 scenario, other countries such as Madagascar, Guinea, Côte d’Ivoire, Ghana and Guinea–Bissau (only in 2070) show losses between 65% to 83%. In the SSP 3.70 scenario (Fig. 5B), the same countries recorded a similar loss of rice production to the SSP 2.45 scenario, with the addition of the Central African Republic, Togo and Guinea–Bissau. Some countries, such as the Democratic Republic of Congo and Cameroon, show average loss values (46–54%) only in the SSP 3.70_2070 scenario.

Concerning the ratio between production loss and GDP (Figs. 5A, 5B; stacked-bars), the condition is still different, as emerging from the heatmap. In the SSP 2.45 narrative (Fig. 5A), the countries showing the most considerable losses (~−0.95%) are the Central African Republic, Guinea–Bissau, Liberia, Madagascar, Gabon, the Democratic Republic of Congo and Mozambique (respectively, in descending order). In the SSP 3.70 scenarios (Fig. 5B), harsh economic losses are predicted to reach a net value of −4.8%, with the Central African Republic (highest absolute values), Guinea–Bissau, the Democratic Republic of the Congo, Chad, Liberia, Mozambique, Madagascar and Gabon being the more affected by the loss. The Central African Republic is the only country showing a notable increase in 2070 (+0.82%) (Fig. 5B).