Species distribution modeling for disease ecology: A multi-scale case study for schistosomiasis host snails in Brazil

Abstract

Species distribution models (SDMs) are increasingly popular tools for profiling disease risk in ecology, particularly for infectious diseases of public health importance that include an obligate non-human host in their transmission cycle. SDMs can create high-resolution maps of host distribution across geographical scales, reflecting baseline risk of disease. However, as SDM computational methods have rapidly expanded, there are many outstanding methodological questions. Here we address key questions about SDM application, using schistosomiasis risk in Brazil as a case study. Schistosomiasis is transmitted to humans through contact with the free-living infectious stage of Schistosoma spp. parasites released from freshwater snails, the parasite's obligate intermediate hosts. In this study, we compared snail SDM performance across machine learning (ML) approaches (MaxEnt, Random Forest, and Boosted Regression Trees), geographic extents (national, regional, and state), types of presence data (expert-collected and publicly-available), and snail species (Biomphalaria glabrata, B. straminea, and B. tenagophila). We used high-resolution (1km) climate, hydrology, land-use/land-cover (LULC), and soil property data to describe the snails' ecological niche and evaluated models on multiple criteria. Although all ML approaches produced comparable spatially cross-validated performance metrics, their suitability maps showed major qualitative differences that required validation based on local expert knowledge. Additionally, our findings revealed varying importance of LULC and bioclimatic variables for different snail species at different spatial scales. Finally, we found that models using publicly-available data predicted snail distribution with comparable AUC values to models using expert-collected data. This work serves as an instructional guide to SDM methods that can be applied to a range of vector-borne and zoonotic diseases. In addition, it advances our understanding of the relevant environment and bioclimatic determinants of schistosomiasis risk in Brazil.

Fig 1. Biomphalaria presence points by species (color) and source (shape), thinned to 1 km.

A) National, B) Minas Gerais, C) São Paulo. Maps were built in R (version 4.2.2) using shapefiles from the geobr package [80].

Fig 2. Large variation in snail suitability probabilities at a national scale.

National prediction maps of B. glabrata (A–C), B. straminea (D–F), and B. tenagophila (G–I) suitability probabilities by model type (MaxEnt: A, D, G; Random Forest: B, E, H; Boosted Regression Tree: C, F, I). Each pixel shows the mean value across 10 bootstrapping iterations in which models were provided 80% of the available species presence records. Maps were built in R (version 4.2.2) using shapefiles from the geobr package [80].

Fig 3. Scale and species drive SDM performance metrics more than model type.

Plots of ten-fold spatially cross-validated, out-of-sample AUC values across species (A, B, C), scales (panels), and model types (colors). Plots display mean (point) and +/- standard error (error bars).

Fig 4. State and national models produce substantially different state-level prediction maps.

Minas Gerais prediction maps of B. glabrata suitability probabilities by model type (rows) and model geographic extent (columns) Each pixel shows the mean value across 10 bootstrapping iterations in which models were provided 80% of the species presence records available at a given scale. Parallel prediction maps of B. straminea in Minas Gerais and B. tenagophila in São Paulo can be found in Figs C and D in S1 Text. Compared to national models (Fig 2), at smaller geographic scales it becomes more obvious that suitability probabilities can be highly localized, producing points of high suitability surrounded by areas with low suitability. Maps were built in R (version 4.2.2) using shapefiles from the geobr package [80].

Fig 5

Examples of marginal effects of covariates on suitability probabilities that vary across model type (A), geographic scale (B), and species (C). Partial dependence plots for three covariates (columns) across model types (color), species (top two rows vs. bottom two rows), and scale (first row vs. second and third row vs. fourth).

Fig 6. Variable importance of land use/land cover (LULC) variables can increase at smaller scales.

Proportion of total variable importance averaged across all training folds attributable to distance to high population density and proportion of temporary crop cover. Displayed for all species (A,B, C) and model types (color).

Fig 7. Expert collected and public GBIF data produce visually different suitability maps for B. tenagophila in São Paulo across model type.

Predicted suitability maps with varying input data (columns) supplied to all model types (rows). Compared to national models (Fig 2), at smaller geographic scales it becomes more obvious that suitability probabilities can be highly localized, producing points of high suitability surrounded by areas with low suitability. Maps were built in R (version 4.2.2) using shapefiles from the geobr package [80].

Fig 8