Modelling the spatial distribution of aquatic insects (Order Hemiptera) potentially involved in the transmission of Mycobacterium ulcerans in Africa

Distribution of aquatic hemipteran families

In this study, we focused on the families Naucoridae, Belostomatidae, Notonectidae, Nepidae, Corixidae and Gerridae. The first four are considered biting insects, and thereby more likely to be involved in the transmission of M. ulcerans [8, 20]. Species of the Corixidae and Gerridae, although they are non-biting insects, have been found carrying M. ulcerans, thereby being attributed some role in the transmission of this bacteria [7, 19]. We retrieved georeferenced presence data for these hemipteran families from the Global Biodiversity Information Facility Database (GBIF) [25]. Occurrence data was cross-checked for spatial consistency (i.e. to verify that occurrence sites were correctly mapped according to country and location name) and remapped when necessary using OpenCage geocoder [26] and Google Engine Map. Geographical correction was done using the ggmap and opencage packages in R v.3.3.2.

Environmental datasets used in ecological modelling

Also from CGIAR-CSI, we obtained a raster dataset of elevation at 1 km². This elevation layer resulted from processing and resampling the gridded digital elevation models (DEM) derived from the original 30-arcsecond DEM produced by the Shuttle Radar Topography Mission (SRTM). We used the elevation layer to generate two topography-related datasets: slope of terrain and flow accumulation. Slope was obtained in degrees.

To produce the flow accumulation layer, we initially created a flow direction layer, in which the direction of flow was determined by the direction of the steepest descent, or maximum drop, from each cell in the elevation dataset. This was calculated as follows: change in elevation value / distance * 100. Flow accumulation was then calculated as the accumulated weight of all cells flowing into each downslope cell in the flow direction layer.

where a is the upslope contributing area per unit contour length (or Specific Catchment Area, SCA), which can be approached by using the flow accumulation, and β is the local slope gradient for reflecting the local drainage potential [30].

We also produced continuous surfaces of straight line distance (Euclidean distance) in km to the nearest water body and permanent rivers based on the Global Database of Lakes, Reservoirs and Wetlands [31] and Digital Global Chart [32] respectively.

Land cover data were downloaded from the GlobCover project at the European Space Agency [33, 34]. This global land cover map is derived by an automatic and regionally-tuned classification of a 300-m MERIS FR time series (19 months) and comprises 22 land cover classes according to the UN Land Cover Classification System (LCCS) [35]. We simplified this raster dataset by grouping the 22 land cover classes into 4 major groups; agricultural land (crops), grassland-shrubland, forest areas and woodlands, and others (i.e. bare soil, urbanized areas and snow/rock).

Finally, raster datasets of averaged Enhanced Vegetation Index (EVI) and land surface temperature (LST) for the period 2000–2015 were obtained from the African Soil Information System (AfSIS) project [36]. This project generates time series average products for several environmental indicators such as vegetation indices and LST using MODIS satellite image data collected by the National Aeronautics and Space Administration (NASA). The MOD13Q1 product, which is updated every 16 days at 250m spatial resolution, includes vegetation indices such as Normalized Difference Vegetation Index (NDVI) and EVI [37]. Compared to the NDVI, the EVI minimizes canopy background variations (which can lead to errors in vegetation classification), and maintains greater sensitivity over dense vegetation conditions. Day and night LST data are generated from MOD11A2 products, and have a spatial and temporal resolution of 1 km and 8 days, respectively [38].

Input grids were resampled to a common spatial resolution of 5 × 5 km using nearest neighbour approach, clipped to match the geographic extent of a map of mainland Africa, and eventually aligned to it. Raster manipulation and processing was undertaken using raster package in R v.3.3.2. Final map layouts were created with ArcGIS 10.3 software (ESRI Inc., Redlands CA, USA). Only non-linearly related covariates were considered for subsequent analysis. This variable selection was done by implementing a stepwise procedure to identify a set of non-correlated variables that will have low variance inflation factor (VIF). VIF is calculated as follows: \( {VIF}_i=1/1-{R}_i^2 \) , where \( {R}_i^2 \) is the coefficient of determination of the regression for each variable as a function of all remaining predictors. We retained a variable combination which resulted in a VIF below 10.

Ecological Niche Factor Analysis

As a first step to build the ecological niche models, an Ecological Niche Factor Analysis (ENFA) [39] was performed with Biomapper 4 [40]. ENFA analysis is widely used to explore habitat suitability and compute species distribution maps without absence data [39, 41]. This analysis was used to identify predictors most likely to limit the insects’ distribution, and for reducing dimensionality, thereby avoiding the inclusion of several correlated variables [42]. The output of this analysis includes non-correlated factors, similar to that obtained by a Principal Components Analysis (PCA), with the particularity that factors obtained always have ecological significance [39].

Employing the environmental variables described above as predictors a set of factors for each of the Hemiptera families were obtained. Factors were selected by the broken stick criterion: they were retained if they explained larger variance than the broken stick null model, which assumes that the variance has been divided at random among the factors [43]. This method is only appropriate for continuous data, and thus land cover classification was excluded for this stage.

The marginality factor (MF), the specialization factors (SF), and global marginality and tolerance values of each species were considered to support the interpretation of the niche models. MF is the first component obtained by ENFA and describes how far the family optimum niche is from the mean habitat in the study area, for each predictor. SF are extracted successively from the n-1 residual dimensions. High absolute factor loadings of a variable on SF indicate a more restricted range of the family on the corresponding variable. The global marginality takes account of the MF scores of all the predictors and gives a summary of how much the family habitat differs from the available average conditions. A low value (close to 0) indicates that the family tends to live in average conditions throughout the study area. A high value (close to 1) indicates a tendency to live in the most extreme conditions of the study area. The global tolerance value considers the eigenvalues of all the factors obtained by ENFA and indicates how specialized the family niche is in relation to the overall model area. It can have values between 0 and 1, and the higher is it, the wider is the range of environmental conditions that the species bears. Thus, a value close to 0 indicates a specialist taxon and a value close to 1 indicates a generalist one [39].

For each hemipteran family, variables with a factor loading higher than 0.3 were considered to best define the insects’ ecological niche [41]. Variables with a factor loading above this cut-off were retained to create sets of environmental factors that were used for the distribution modelling of each family.

Developing ecological models

An ensemble of distribution models was generated for each hemipteran family, based on the presence data collected and the environmental factors identified through ENFA analysis. We used seven algorithms available within the BIOMOD framework [44] to obtain those ensembles of predicted distributions: generalized linear models (GLM), generalized additive models (GAM), generalized boosted regression models (GBM), artificial neural networks (ANN), multiple adaptive regression splines (MARS), maximum entropy (MaxEnt) and random forest (RF). These models were run using the parameters set by default in the biomod2 R package [44], except for the GBM models. As per guidelines from Elith et al. [45], the learning rate (lr) and tree complexity (tc), key parameters in GBM models, were set according to the number of observations and testing different values on a subset of samples (75%), using deviance reduction as the measure of success. Overall, lr of 0.005 and tc of 5 were identified as optimal parameters, thereby enabling the model to account for up to 5 potential interactions and slowing it down enough to get the model converged without over-fitting the data. This tuning was undertaken using the gbm package in R v.3.3.2.

All these models are intended to discriminate the suitability of the environment for the presence of a particular organism, and for this they need to be trained with presence and absence records. When there are no available absences records, an alternative is to generate background points or pseudo-absences [46]. We generated sets of the same number of background points as presence data compiled for every hemipteran family. Background points were randomly generated accounting for the underlying geographical bias on the occurrence data, as previously recommended [47]. For this, we created a sampling bias surface by counting the number of occurrence records within each grid cell (5 × 5 km resolution) and then extrapolated these data across Africa using kernel density estimation using the R packages kernlab, ks and sm. Lastly, we generated the background points from random locations weighted by the sampling bias surface [48, 49].

Models were calibrated using an 80% random sample of the initial data and evaluated against the remaining 20% data using the area under the curve (AUC) of the receiver operation characteristic (ROC) and the true skill statistic (TSS) [50]. Projections were performed 80 times, each time selecting a different 80% random sample while verifying model accuracy against the remaining 20%. Verification or internal evaluation does not allow for assessment of the predictive performance of the models - independent evaluation data would be required for this purpose - but it provides a measure of internal consistency of the models. Unfortunately, the scarcity of occurrence records collected did not allow for independent cross-validation of the final models. The evaluation statistics (AUC and TSS) were used to select the models to be assembled on the basis of matching between predictions and observations. Here, models with AUC < 0.8 or TSS values < 0.7 were disregarded when constructing the final assemble model.

The final assemble model was obtained by estimating the median of probabilities across the selected models per Hemiptera family and per grid cell. This approach was used alternatively to the mean because is less sensitive to outliers [51]. The range of uncertainties obtained with the seven modelling techniques was also calculated by estimating the confidence intervals across the ensemble for each hemipteran family and per grid cell.

Data sources on Buruli ulcer occurrence

Two difference datasets were used to explore the potential geographical overlap between BU occurrence and the predicted distribution of aquatic Hemiptera insects. A data table with a list of 91 georeferenced villages from southwestern Ghana that reported BU cases during 2007–2010 was obtained from a recent publication [52]. A second dataset was provided by the National Committee for Yaws, Leishmaniosis, Leprosy and Buruli ulcer Control in Cameroon (Cameroon Ministry of Health) and included a list of 414 communities that reported BU cases during 2003–2015. We decided to retain for this analysis those communities providing accurate geographical coordinates and reporting PCR-confirmed BU cases in early stages of disease (to exclude potentially imported cases). A map displaying the georeferenced locations from Cameroon and Ghana has been included in a supplementary file (Additional file 1: Figure S1).

General description of datasets

A total of 2217 records of occurrence were compiled for the six hemipteran families targeted in this work: 439 for the Naucoridae (Additional file 2: Table S1); 233 for the Belostomatidae (Additional file 2: Table S2); 469 for the Notonectidae (Additional file 2: Table S3); 234 for the Nepidae (Additional file 2: Table S4); 234 for the Gerridae (Additional file 2: Table S5); and 608 for the Corixidae (Additional file 2: Table S6). Based on the source of data, 52.3% (1589/2217) were extracted from the GBIF and 47.7% (1058/2217) through systematic literature searches (Table 1). The proportion of occurrences by source varies across the families, as shown in Table 1. Only 18 recorded occurrences were associated with BU studies and in 10 of them M. ulcerans was found (by PCR test) within studied specimens.

Most of the collected occurrences are distributed across southern, western and northern Africa, although geographical distribution varies between families (Fig. 1). Naucoridae and Belostomatidae occurrences are geographically more widespread, extending to eastern Africa. For others, such as the Gerridae and Corixidae, recorded occurrences are restricted to northern and southern Africa. Figure 1 also shows the distribution of background points (“pseudo-absences”) which, as detailed in method section, were generated to account for the geographical bias in the distribution of occurrence records.

ENFA analysis

The Belostomatidae scored the lowest value for the global marginality index (Table 3), indicating that this family would be most likely to occupy average or normal habitats, according to the analyzed environmental predictors, throughout Africa (Table 3).

Flow accumulation was predicted to condition the distribution of the families Corixidae, Gerridae, Nepidae and Notonectidae. Considering the global tolerance index (Table 4), the Gerridae and Nepidae seem to have more restricted ecological niches in Africa than the rest of the hemipteran families.

Environmental suitability

A suite of 560 ensembles obtained by different modelling approaches (GLM, GAM, GBM, RF, ANN, MARS and MaxEnt) - 80 per each - were run with a random sample of 80% of the total occurrences and background points for each hemipteran family. Only those ensembles with a high predictive performance (TSS > 0.7 or AUC > 0.8) contributed to the final model. Additional file 1: Table S1 shows the median value of statistic indicators and 95% confidence interval (95% CI) of their distribution. Boosted regression trees (GBM) and random forest (RF) consistently outperformed the other modelling approaches, thereby primarily contributing to the final ensemble model. However, a number of other models also informed the final predictive models, as can be inferred from their TSS and AUC values (Additional file 1: Table S1).

Figures 2 and 3 display the environmental suitability for the presence of each hemipteran family across Africa and the uncertainty of prediction, this being calculated as the range of the 95% confidence interval in predicted probability of occurrence for each pixel and rescaling to a 0–1 scale. Coastal areas of West Africa, and the south-east of South Africa were predicted to be highly or moderately suitable for all families, and parts of Central Africa and East Africa were suitable for all biting families. Predicted suitability for the families Nepidae, Notonectidae and Gerridae in endemic areas of West and Central Africa is associated with high uncertainty, most likely associated with a lower number of presences recorded in these regions, but suitability for the Naucoridae, Belostomatidae and Corixidae is accompanied by moderate levels of uncertainty. The Belostomatidae was predicted to have the widest suitability range. The Sahara was predicted to be generally unsuitable for all families, but with pockets of moderate suitability for the Naucoridae and Belostomatidae. These predictions were associated with high certainty. Additional file 3: Figures S1-S24 provide more details on findings from ecological modelling for each hemipteran family: (i) maps of predicted suitability (median probability of occurrence) and 95% CI lower and upper bound maps; (ii) binary maps displaying the areas where their occurrence is more likely based on optimal cut-off; and (iii) a detailed description of marginal effect plots for GBM and RF models.

In summary, the families Naucoridae and Belostomatidae show a ubiquitous geographical distribution across Africa, except throughout drier and hotter areas such as the Kalahari and Sahara deserts. The other two families of biting Hemiptera, the Notonectidae and Nepidae, would have a more restricted distribution, being more predominant in western and southern Africa and narrowly present in middle Africa. All these four families are widely distributed across coastal areas of West Africa, although the Belostomatidae distribution would extend towards more inland ecosystems of western African countries. Of the hemipteran families considered non-biting, the Corixidae would be spread across Africa, even in desert areas from Sahel, whilst the Gerridae, a well-known cosmopolitan hemipteran family, would have a more limited distribution, according to our environmental models.

Additional file 1: Figures S2-S3 show the relative contribution (average) of selected environmental variables to the final predictive models. Although there are differences between families, land surface temperature appeared to be the main factor driving the distribution of all hemipteran families, followed by other temperature-related parameters such as mean temperature in the warmest and coldest quarters. According to the marginal effect plots (Additional file 3: Figures S3, S4, S7, S8, S11, S12, S15, S16, S19, S20, S23, S24), all these aquatic Hemiptera seem to thrive under mean temperature between 15–22 °C and near to water courses. However, the suitability for their presence declines sharply when LST exceeds 30 °C. The family Naucoridae, however, seems to tolerate even higher temperature according to their marginal effect plots (Additional file 3: Figures S3 and S4). In general, optimal ecological niches for these hemipteran families are characterized by moderate annual average precipitation (i.e. 350–1000 mm/year). The four biting Hemiptera families show preference for transformed ecosystems such as agricultural landscape and human settlements over forest or grasslands.

Figure 4 shows the potential for concurrent geographical distribution of all these hemipteran families. Our models predict the geographical overlapping of the biting Hemiptera families in broad areas of western Africa, southern and northern Africa, highlands of Ethiopia and Uganda, the northernmost and south-west coasts of the continent, the east coast of Madagascar and in pockets across eastern Africa.

Geographical overlap with Buruli ulcer

A simple exploratory analysis overlaying geolocated communities that have reported BU cases and our predictive models shows that M. ulcerans and species of these aquatic hemipteran families might coexist in the same ecosystems (Fig. 5, Additional file 1: Figure S1). Although all BU-endemic areas in both countries were suitable for at least one hemipteran family, our analysis suggests a marked heterogeneity in this relationship: BU-endemic areas in Cameroon were only suitable for two of the six hemipteran families, whereas BU-endemic areas in Ghana were predicted to be suitable for all six families (Fig. 5). There also seems to be spatial heterogeneity in the diversity of hemipteran families present in BU foci from West and Central Africa. Over 75% of BU-endemic communities from Ghana here studied would have the environmental conditions for five out of the six hemipteran families to prosper, while in Cameroon, less than 10% communities of the BU-endemic communities would have suitable conditions for 4 of the hemipteran families. This difference is even more significant when we restrict the analysis to biting insects (Fig. 6).