The study area is located on the hilly uplands (>1,600 m a.s.l.) of the Edward-George branch of the East African Rift valley in western Uganda (Figure 1). It is bordered by the Rwenzori Mountain range in the North-West, Lake Albert in the North and Queen Elizabeth National Park in the South. It contains over 90 small crater lakes, most of which formed between 4,000 and 10,000 years ago as a result of faulting and volcanic eruptions (Vinogradov, 1980; Schumann et al., 2015). Lake Ntambi, which formed c. 50,000 years ago is an exception, (Dirk Verschuren, pers. comm.).

The crater lakes of western Uganda have been repeatedly promoted as an ideal model system for studying large environmental gradients in limnological characteristics in a setting that allows meaningful comparisons due to shared geological history (e.g. Melack, 1978; Mills and Ryves, 2012; Saulnier-Talbot and Lavoie, 2018). The crater lakes region is one of the most densely populated rural areas in sub-Saharan Africa (Hartter et al., 2012), with a population growth rate of 3.3% annually (Gelorini et al., 2012). The population growth is coupled with increasing human impact through settlement, fishing, agriculture on the fertile volcano soils, swimming and water extraction for irrigation and domestic use. As a result of the high population density, the prevalence of (human) schistosomiasis in the region has dramatically increased (Kabatereine et al., 2004; Kabatereine et al., 2006; John et al., 2008; Stanton et al., 2017). In 2014, the crater lakes region covered parts of the 73 out of 112 districts of Uganda with prevalence of schistosomiasis (Loewenberg, 2014). Nationwide, four million people are estimated to be infected and almost 20 million are at risk of infection (Loewenberg, 2014). The only preventative measures in the area are massive drug administration while environmental transmission interruption is rarely emphasized (Loewenberg, 2014).

Melack (1978) classified the region into four crater lake fields (CLFs, Figure 1). The Ndali-Kasenda field is locatedin the central part of the region, ∼20 km from Fort Portal field in the North and ∼65 km from Bunyaruguru field in the South. The fourth field of Katwe-Kikorongo contains lakes with saline waters (Rumes et al., 2011), unsuitable for mollusks (Tumwebaze et al., 2019).

Our dataset combines information from a total of 56 crater lakes in the three fields, including 32 lakes of the Ndali-Kasenda field (CLF 1), six lakes of the Fort Portal field (CLF 2) and 18 lakes of the Bunyaruguru field (CLF three; see Supplementary Table S1).

We applied a machine learning approach, i.e., random forest (RF; Breiman, 2001), to assess the combined impact of the chosen set of predictors on the distribution of IH snails. Machine learning approaches, and particularly RFs, have gained prominence in classification and regression analyses across various fields of science in recent years (e.g., Huang and Boutros, 2016; Pang et al., 2017; Schonlau and Zou, 2020; Collin et al., 2021; Georganos et al., 2021; Ruiz-Aĺvarez et al., 2021). In classification problems, RFs have been demonstrated to give more accurate predictions than other approaches, such as logistic regression (Boulesteix et al., 2012; Bunyamin and Tunys, 2016; Couronné et al., 2018; Xia et al., 2019; Zhang et al., 2020). Since it is a non-parametric technique, the RF algorithm is not affected by multicollinearity among the predictor variables (Boonprong et al., 2018), which is a common problem in ecology. Also, many RF software packages come with convenient solutions to deal with missing values (Brieuc et al., 2018).

We conducted separate RF analyses to variably predict the presence of Bulinus, Biomphalaria and both genera combined. To assess variation of the potential predictors across different geographical scales, at a larger geographical extent, we ran analyses for the overall and complete dataset combining data from all the three crater lake fields. On the scale of individual crater lakes regions, we ran analyses for two subsets of CLF 1 and CLF 3. The analysis for Fort Portal field (CLF 2 subset) was not performed because it contains only six lakes. All analyses were done in the R statistical environment v. 4.0.3 (R Core Team, 2020), using the packages randomForest v. 4.6-14 (Liaw and Wiener, 2002), rfUtilities v. 2.1-5 (Evans and Murphy, 2019) and rfPermute v. 2.1.81 (Archer, 2020).

We performed imputation to fill missing data prior to further data analyses, using the function “rfImpute” in the package randomForest, which uses the RF algorithm to obtain weighted averages of the available observations. This was done for the predictors; calcium, magnesium and water depth, for which no data were available for 18, 16 and 16 lakes, respectively. Overall, missing data added up 4.7%. The “rfPermute” algorithm was used to assess variable importance in each RF model via permutation. The algorithm creates decision trees from the original dataset by random sampling of rows (i.e., lakes) without replacement. At each node, two-thirds of the rows are taken as training data to create the model, the remaining one-third is taken as so-called out-of-bag (OOB) sample and is used to make predictions and test for the performance of the model (Breiman, 1996; Breiman, 2001). Variables were permuted 100 times over 1,000 decision trees.

Model performance was additionally assessed via cross validation. This approach was chosen over the standard train–test data procedure because of the comparably low number of lakes in the dataset. Cross validation is a commonly used resampling method to assess the generalization potential of a model and to avoid overfitting (Berrar, 2019). We report here the kappa coefficient K, which determines the model’s predictive accuracy, i.e., it gives the percentage of the data that is in agreement with the model (and is thus the opposite of the validation error). We adopted the suggestion of Viera and Garrett (2005), who regarded kappa values of 0.81-1 as excellent and 0.61–0.80 as substantial agreement, to evaluate our model performances.

The relevance of individual parameters to the overall RF models was assessed using two importance metrics, i.e., the mean decrease in accuracy (MDA) and the mean decrease Gini (MDG) (Calle and Urrea, 2011; Huang and Boutros, 2016). Due to discrepancies in ranking results between MDA and MDG, the results of a single metric are not completely exhaustive (Strobl et al., 2007; Liu et al., 2011). Therefore, we included both metrics but discussed in detail only those parameters that are found significant by both MDA and MDG. We used partial dependence plots to visualize the relationships and marginal effects of individual predictor variables (Friedman, 2001; Evans et al., 2011).

Geographical variables were found to be among the most important controls for IH snail distribution, but the individual parameters depend on the scale of observation and the taxon in question. Latitude, longitude and distance to the nearest large lake are relevant for the occurrence of Biomphalaria as well as both genera combined. In contrast, Bulinus is much less influenced by geographical variables. Although latitude and altitude were identified significant, no positive relationship was found for the first (Figure 3), whereas the latter was of comparatively low importance (Figure 2). A similar trend is found in the Ndali-Kasenda crater lake field (CLF 1), where the distribution of Biomphalaria is strongly shaped by the distance to the nearest crater lake, whereas the presence of Bulinus is mostly related to the diversity of the associated mollusk fauna (Supplementary Figure S2).

The strong positive relationship between the presence of Biomphalaria and latitude as well as longitude reflects the uneven occurrence of the genus across the study area (Figure 1). Particularly the rarity of the genus in the Bunyaruguru crater lake field (CLF 3) indicates that dispersal limitations and/or environmental filtering (e.g., Hauffe et al., 2016b) might constrain its distribution.

The generally positive association between the occurrence of Biomphalaria and distance to the nearest large lake is surprising at first, considering that distance is related to biogeographical isolation and colonization potential (Covich, 2010). Generally, these snails have a high passive dispersal capacity, a high reproduction rate and short generation times (Brown 1994). This combination is pronounced particularly in pulmonate snails such as Biomphalaria and makes them prime colonizers (Kappes and Haase, 2012; Kappes et al., 2014). The crater lakes are hydrologically disconnected, the possible longitudinal dispersal means are natural vectors such as birds (Kappes and Haase, 2012) and humans, through the attachment to fishing nets and/or boats.

A possible explanation for the positive association with distance concerns the taxonomic resolution. Here, we investigated the presence of genera, which may overprint patterns of dispersal and colonization of individual species of Bulinus and Biomphalaria. Also, Euclidean shore-to-shore distances might not reflect real dispersal means or pathways, because hydrological connectivity varies considerably between the large lakes and the crater lake fields. For example, the Fort Portal field (CLF 2) is connected to the Lake George–Lake Edward system via the Mpanga River but not to the geographically closer Lake Albert. The pattern may also results from choosing only four major lakes but disregarding smaller ones and river systems. Finally, in colonization processes, stochastic components might play an important role. For instance, regular episodes of shifting mixing regimes in some crater lakes (De Crop and Verschuren, 2021) could lead to low dissolved oxygen and consequently to the demise of most aquatic life. Based on genetic relationships of Bulinus, a previous study suggested that the large lakes acted as potential sources for (re)-colonization of the crater lakes (Tumwebaze et al., 2019), a scenario that cannot be tested here given the lack of genetic data. Fast re-appearance or aestivation in fluctuating environments is an intrinsic ecological feature of many Bulinus species, including the ones found in the study area (Watson, 1958; Brown, 1994). As such, the colonization and re-colonization patterns might be more complex and faster than previously anticipated.

Altitude was found to have limited influence on the distribution of IH snails. Several previous studies have indicated a shift in recent years of intestinal human schistosomiasis transmission towards higher altitudes (Kabatereine et al., 2004; John et al., 2008; Stanton et al., 2017). This is also supported by our findings showing that IH snails occur up to approx. 1,600 m a.s.l., and thus considerably higher than the previously presumed threshold of 1,400 m a.s.l. (Kabatereine et al., 2004). Future studies will have to establish an upper limit for both forms of schistosomiasis, because both IH snail genera have been found at altitudes above 2,000 m a.s.l. in Uganda (Stanton et al., 2017, unpubl. data).

A strong control for the distribution of both Bulinus and Biomphalaria is the diversity of the associated mollusk fauna. The strong positive relationships in all three models (Figure 3) indicate that both genera are more likely encountered in generally diverse systems. One may expect the opposite, that a higher number of species (especially in a small crater lake) results in higher competition (e.g. Svanbäck and Bolnick, 2007; Hauffe et al., 2016b) and thus a lower chance of successful establishment. Instead, the positive association with diversity indicates that the IH snail genera are present in environments that provide favorable conditions for mollusks generally.

The associated mollusk fauna consists predominantly of pulmonates (Radix natalensis and several species of planorbids other than Biomphalaria), which are good colonizers in general and are characterized by a high productivity and shorter generation times (Kappes et al., 2014). Malacological surveys showed the presence of Biomphalaria sudanica, Biomphalaria pfeifferi, Bulinus forskalii, Bulinus globosus, and Bulinus tropicus (Tumwebaze et al., 2019; Tumwebaze et al., unpubl. data). All Biomphalaria species in the region are regarded susceptible to S. mansoni. Only B. forskalii and B. globosus were identified as host species for human schistosomes and the majority of Bulinus in the crater lakes region were B. tropicus, which is an important IH snail for livestock schistosomes. The non-pulmonate species are common and widespread regional species with no particular ecological requirements (Brown, 1994). More specialized taxa such as Bellamya or Cleopatra or unionid bivalves were not found in the crater lakes. The mollusk association in the crater lakes could still be seen as depauperate. The most obvious variation in habitat conditions and thus potentially determining mollusk associations is in the limnological characteristics of the crater lakes. Although we gathered a variety of different parameters concerning water chemistry and climatic conditions, the associations with individual parameters do not allow drawing a multifactorial picture. Moreover, the parameters measured for each crater lake do not cover seasonal variability (except for climate and lake area variation).

Only selected hydrochemical variables had an influence on the distribution of IH snails in our study region. Water conductivity was a significant driver across the entire study region and in CLF one for Bulinus. For Biomphalaria, oxygen was found a significant driver by one metric (MDA), but its relative importance was low (Figure 2). Similarly, MDA found the presence of both genera linked to magnesium and calcium, but again with low importance values.

Water conductivity determines the ionic strength of the concentration of dissolved solids including calcium and magnesium (Cormier et al., 2013). An expanding distribution with increasing conductivity was also suggested by Camara et al. (2012). This might reflect the presence of dissolved ions (e.g., calcium), which stimulate shell development for snail species (Dillon, 2000; see also below). The non-linear relationship with a sharp increase around ∼400 μS/cm followed by a nearly steady but weak decline indicates that a certain threshold must be met to allow the establishment of Bulinus in a lake. Higher electric conductivity is, in turn, less favorable. Conductivity as a complex factor integrating several chemical components has been identified as a determinant for Bulinus mortality and thus poses an important constraint on its occurrence (Brown, 1994; Marie et al., 2015).

The predicted probabilities for the presence of both genera in relation to magnesium and calcium concentrations, which were, however, only found significant by one importance metric (MDA), demonstrate that a certain threshold must be surpassed in both factors to promote the establishment of IH snails. Especially the relationship with calcium is not surprising, considering that a low concentration would constrain snail growth, fecundity, survival rate and reproduction, which, in turn, limit snail distribution (Dillon, 2000; Brodersen and Madsen, 2003). The concentration of calcium and magnesium and their ratio in the water both affect the presence and life cycle performances in southern African streams with thresholds at the lower and higher ends of concentrations (Brown, 1994).

Climatic factors seem to have surprisingly little effect on the distribution of IH snails in Ugandan crater lakes. Precipitation was found a significant but rather weak predictor in the model with Biomphalaria, whereas air and water temperature do not appear influential. Generally, climatic conditions and climate change are known to be important predictors for mollusk species distribution on larger geographical scales (Marcogliese, 2008; Stensgaard et al., 2019) and were also found to be relevant for IH snails. Temperature influences the survival and reproduction rates of the snails (Paull and Johnson, 2011; McCreesh et al., 2014; Kalinda et al., 2017). Higher rainfall causes more runoff into freshwater ecosystems increasing the supply of organic matter serving as food for the snails, which, in turn, promotes growth and fecundity (Madsen et al., 1987; Camara et al., 2012; Nyström Sandman et al., 2013; but see also discussion in David et al., 2018 for opposite associations). This may also explain the relationship with precipitation in our case. However, the lack of a (strong) association in our models is likely owed to the constrained geographical scale of our observations. The crater lakes are located in the same climatic zone in western Uganda, which is why precipitation and air temperature vary only little across the study region. Whereas air temperature seems to be irrelevant on that scale, the relationship with precipitation indicates that even small variations might have a significant impact on the IH snail distribution.

A series of parameters was found to be non-significant in any of our models. Here, we briefly discuss the potential reasons for these findings, especially in the light of conflicting results reported in the literature. In contrast, previous studies found relationships between the distribution of Biomphalaria and water pH, whereas both oxygen and water pH have been suggested to determine the occurrence of Bulinus (Yirenya-Tawiah et al., 2011; Stensgaard et al., 2013; Marie et al., 2015; Mahmoud et al., 2019; Alhassan et al., 2020). However, no such influences were noticed despite of the influential factors of high oxygen concentration and water pH of the lakes studied (Supplementary Table S1). Possibly, the high carbon dioxide emission from decomposing submerged vegetation and organic matter, together with the presence of other dissolved ions could have indirectly affected such relationships (Tchakonté et al., 2014). Lake area variation, being a measure of ecosystem stability, seems also of little importance. Apparently, most of the lakes show minimal fluctuations through time, and the few that occasionally dry up could not influence the distribution or recolonization, which is facilitated by the short distances between the crater lakes. The mollusk associations of mostly opportunistic species found is also supporting this interpretation (see above).

Human impact, quantified by both land use and population density around the crater lakes, has apparently no impact on snail distribution. This result is rather unexpected, because humans are often involved in introduction of snails into new environments as passive dispersal vectors (Kappes and Haase, 2012). In the Ugandan crater lake fields, extensive anthropogenic activities (e.g., multipurpose water fetching, watering for livestock, littering and pollution) are limited to lower, relatively flat and thus more easily accessible shores of most of the crater lakes. As such, one may have expected an impact on snail occurrences. Moreover, only about 26% of the crater lakes are located in national parks (Queen Elizabeth NP and Kibale NP), thus three-quarter of the lakes are accessible and utilizable by humans. Despite these constraints, some reserves and parks like Kibale NP increasingly face migration of settlers and extensive cultivation (Hartter et al., 2015). Perhaps the prevailing population pressure affects the prevalence of infected snails, but it does not influence their presence in general. Follow-up studies should focus specifically on this aspect for mollusks communities and the prevalence of schistosomiasis. Interestingly, in other non-gastropod taxa such as cladocerans (Rumes et al., 2011) and fungi communities (Gelorini et al., 2012), negative effects of land use change and anthropogenic pressure have already been demonstrated for the crater lakes region. With increasing human activities in the Ugandan crater lakes region (including tourism), we expect further changes in biological communities of the lakes as well as increasing cases of human schistosomiasis among both local communities and visiting travelers (Lachish et al., 2013) in the near future.

Despite the numerous advantages machine learning approaches like RFs offer, one disadvantage concerns the non-linearity of the approach, resulting in a limited prediction power outside the data range. Therefore, we cannot extrapolate our findings to other datasets or regions. Furthermore, small datasets (e.g. CLF 2, Biomphalaria subset with CLF3), decisively limit RF classification due to insufficient data. Optimally, future studies need to include data from a larger geographical scale, involving a greater variation in the ranges of predictor variables, to provide more general predictions for the controls of IH snail distributions. Using machine-learning algorithms like RFs on a comprehensive dataset will eventually facilitate more general conclusions about the importance of individual predictors (or sets of predictors) for the presence of IH snails. Future studies may focus on applying this approach to map infection risk areas, especially in comparison to areas with actual prevalence records and those where preventative measures are in place.