This study was carried out in KZN province, South Africa. It is located in the southeast of the country, with a long shoreline on the Indian Ocean and sharing borders with three provinces (Free state, Eastern Cape, and Mpumalanga) and three countries; Mozambique, Eswatini, and Lesotho. It is the second-most populous province in South Africa with a population of about 11 million people and an area approximately 94,361 km² (Khumalo et al., 2022). KZN comprises of ten district municipalities: Amajuba, iLembe, Ugu, uMgungundlovu, uThukela, uMkhanyakude, King Cetshwayo (formerly uThugulu), uMzinyathi, Harry Gwala (formerly Sisonke), and Zululand; and one metropolitan municipality: eThekwini (Figure 1). KwaZulu-Natal’s climate differs from subtropical to temperate with four seasons: rainy (December–February), post-rainy (March–May), cold/dry (June–August), and hot/dry (September–November). KwaZulu-Natal province is characterized by major rivers such as Mfolozi, Nyalazi, Tugela, Pongolo, and uMvoti amongst others. The province consists of three distinct geographical areas, namely, the lowland Indian Ocean coastal region, the central Natal Midlands and the mountainous areas of the Drakensberg and Lebombo Mountains (Ndlovu, 2013). Primary economic activities in the province include agriculture, mining and quarrying, manufacturing, trade, and tourism.

We a conducted malacology survey in the 11 districts of KZN province for 12 months (September 2020–August 2021). Sampling sites were determined by proximity to schools (randomly selected) where parasitological surveys were carried out and where human water contact activities occurred. This was done because this study is a subset of a larger study and we needed to be able to trace infections. Snail sampling were done once in each season (that is, quarterly) at each site for 15 min by trained field staff. The trained field staffs used the combination of handpicking and the scoop nets described Appleton and Miranda (Appleton and Miranda, 2015) to collect the snails. The collected snails were identified morphologically using the key developed by Brown and Kristensen (Brown and Kristensen, 1989). Bulinus globosus and B. pfeifferi snails from each site were labelled, put in a container with water from the water body and transported to the processing site. Global positioning system was used to record the geographical coordinates of each site with the occurrence of B. globosus and B. pfeifferi. A total of 45 occurrence points of human schistosomiasis intermediate snail hosts were obtained in KwaZulu-Natal comprising of 38 occurrence points B. globosus and 17 occurrence points for B. pfeifferi (Figure 1).

Bioclimatic variables are important in determining species’ habitats. Therefore, a total of 17 bioclimatic variables were downloaded with a spatial resolution of 30 arc-second (923 m) from the AFRICLIM version 3 database (https://webfiles.york.ac.uk/KITE/AfriClim/GeoTIFF_30S/) for current and future projections of suitable habitats for B. globosus and B. pfeifferi (Table 1). These variables were downloaded because they are known to have influence on the habitat suitability of human schistosomiasis intermediate host snails (Mahmoud et al., 2022). We used the AFRICLIM ensemble model that comprised of ten general circulation models (GCMs) (GCM: CCCma-CanESM2, MPI-M-MPI-ESM-LR, CNRM-CERFACS-CNRM-CM5, ICHEC-EC-EARTH, NOAA-GFDL-GFDL-ESM2M, CSIRO-QCCCE-CSIRO-Mk3-6–0, IPSL- IPSL-CM5A-MR, MIROC-MIROC5, MOHC-HadGEM2-ES, and NCC-NorESM1-M), reduced to five regional climate models (RCM: CCCma-CanRCM4_r2, CLMcom-CCLM4-8-17_v1 (4 GCMs), DMI-HIRHAM5_v2, KNMI-RACMO22T_v1 (2 GCMs), SMHI-RCA4_v1 (10 GCMs)), the WorldClim contemporary baseline, under two representative concentration pathways of the IPCC-AR5 (RCP4.5 and RCP8.5). This ensemble model was used because it provides robust assessments of expected climate variability and reduces bias effectively (Platts et al., 2015). To predict the future distribution of habitat suitability of intermediate host snails in KwaZulu Natal, the Representative Concentration Pathways (RCPs) was used. RCPs are trajectories of four greenhouse gas emissions and concentrations adopted by the International Panel on Climate Change (IPCC). These pathways help to determine four possible future climates based on the amount of greenhouse gases emitted by 2085. RCP 2.6 is used in cases of minimum greenhouse gas emissions scenarios. RCP 4.5 and RCP 6.0 represents the intermediate scenarios of greenhouse gas emissions while RCP 8.5 represents the case of very high greenhouse gas emissions (Change, 2023). If additional efforts are not made to reduce emissions, the pathway will range between RCP 6.0 and 8.5 while RCP 2.6 represents a scenario where all efforts to reduce emission works and keeps global warming below 2°C above pre-industrial levels. Since RCP 2.6 is almost impossible to attain, this work uses RCP 4.5 and an extreme case of RCP 8.5 to predict changes in the spatial distribution of suitable habitats for schistosomiasis intermediate host snails.

The MaxEnt version 3.3.3k (https://www.cs.princeton.edu/∼schapire/maxent) was used to create the habitat suitability maps for B. globosus and B. pfeifferi. For each species, the training data comprised of 75% of the presence data that was randomly selected while the test data comprised of the remaining 25% which was used in predicting the ability of the model (Phillips, 2007). Variable contribution analysis that gives a heuristic estimate of the relative contributions of the environmental variables to the MaxEnt model was carried out. In addition, a jackknife procedure, implemented in MaxEnt, was used to quantify each environmental variable’s explanatory power (importance) (Manyangadze et al., 2016). It reveals the predictive power of a model when trained under different scenarios (without a variable, with only a specific variable, and with all variables).

To assess the performance of the fitted MaxEnt model in appropriately discriminating suitable and unsuitable habitats for B. globosus and B. pfeifferi, the area under the curve (AUC) of the receiving operating characteristics (ROC) was used. AUC values range from 0 to 1, where a model with an AUC value of 1 indicates perfect discrimination whereas where a model with AUC values ≤ 0.5 indicate no difference in the discriminative power of the model versus a random prediction. A model with an AUC value closer to 1 has better discriminant power. A model with an AUC above 0.75 indicates sufficient discrimination and is useful (Ferson et al., 2000; Phillips and Dudík, 2008).

The habitat suitability index for B. globosus and B. pfeifferi were calculated. These values range between 0 and 1, where 0 implies that the habitat is unsuitable for the species to be found there and higher habitat suitability index values indicate higher habitat suitability for species occurrence (Du et al., 2021).

Among the 13 bioclimatic variables used in the prediction of habitats suitable for the distribution of B. globosus, the four highest variables in terms of contribution were mean temperature warmest quarter (37.6%), the number of dry months (32.6%), mean diurnal range in temperature (10.8%), and isothermality (6.7%), which quantifies the degree of day-to-night temperatures oscillation relative to the summer-to-winter (annual) oscillations. On the other hand, annual moisture index (34%), mean temperature warmest quarter (21.5%), isothermality (20.5%), and the number of dry months (7%) were the four variables that contributed most in the prediction of suitable habitats for the distribution of B. pfeifferi out of the 10 bioclimatic variables included in the model. Table 1 shows all the variables that were considered in this study and presents the percentage contribution of each bioclimatic variable in predicting the habitat suitability for B. globosus and B. pfeifferi distribution. Variables that had zero percentage contribution to the model were excluded (Table 1).

The result of the jack-knife test of variable importance is shown in Figure 2 below. The four environmental variables with the highest regularized training gain when used in isolation for modeling the habitat suitability of B. globosus are mean temperature warmest quarter (BIO10), number of dry months (DM), length of longest dry season (LLDS), and mean annual temperature (BIO1). Number of dry months (DM), mean temperature warmest quarter (BIO10), annual moisture index (MI), and Rainfall wettest quarter (BIO16) are the four environmental variables with the highest regularized training gain when used in isolation for modeling the habitat suitability of B. pfeifferi. This implies that the top four important variables above have the most information that is not present in the other variables. The environmental variable that decreases the gain the most when omitted is BIO10 for B. globosus and B. pfeifferi, which therefore appears to have the most information that is not present in the other variables. In Figure 2, the light blue bars (without variables) are not longer than the red bars (with all variables), this shows that the predictive performance of the model is better when the corresponding variables are used.

The area under the curve (AUC) of the receiving operating characteristic (ROC) for B. globosus and B. pfeifferi are 0.791 and 0.896, respectively (Figure 3). The results obtained were above 0.75 indicating that the MaxEnt model had excellent predictive accuracy and satisfactorily reflects the distribution and habitat suitability for B. globosus and B. pfeifferi (Ferson et al., 2000; Phillips and Dudík, 2008).

The habitat suitability of B. globosus and B. pfeifferi in KwaZulu Natal under current climate conditions is presented in Table 2 and Figure 4. A total area of 14 096.01 k m 2 ; 15.71% and 14 973.82 k m 2 ; 16.68% in KZN is suitable habitat for B. pfeifferi and B. globosus, respectively. Current suitable areas for the distribution of B. globosus are areas along the coastlines (uMkhanyakude, King Cetshwayo iLembe, eThekwini, and Ugu districts), the northern part of the province (uMkhanyakude and Zululand districts), the western part of the province (uThukela district). The mountainous areas of the Drakensberg and Lebombo mountains (Harry Gwala and Amajuba districts) and the central Natal Midlands (uMgungundlovu and uMzinyathi districts) are currently unsuitable areas for the distribution of B. globosus. Habitats in the eastern part of the province (uThukela), districts along the coastlines excluding eThekwini, and the northern part of uMkhanyakude district are suitable for the distribution of B. pfeifferi. There are more suitable habitats for B. pfeifferi compared to B. globosus in uThukela district. The non-suitable habitats for B. pfeifferi are like that of B. globosus. We used the “maximum training sensitivity plus specificity” as the threshold to distinguish “suitable” from “unsuitable” habitats to obtain the best results as recommended by Hu and Jiang (Hu and Jiang, 2011). This was 0.542 for B. globosus and 0.480 for B. pfeifferi.

The habitat suitability predictions for B. globosus under the RCP4.5 and RCP8.5 climate scenarios in 2085 are shown in Figure 5. A significant difference was seen between the current suitable habitats and those predicted in 2085 across both RCP4.5 and RCP8.5 climate scenarios. The prediction shows a great expansion in suitable habitats for B. globosus towards the midlands (uMzinyathi, uMgungundlovu, and the entire uThukela districts). The predicted suitable habitats for the RCP8.5 climate scenario show more expansion compared to the RCP4.5 climate condition with the former projection having districts situated in the north-western and south-western regions of the province as suitable habitats. Under the RCP4.5 and RCP8.5 climate scenarios approximately 95.5% and 97.2% of the currently suitable areas will remain suitable while 32.4% and 69.3% of the currently unsuitable habitats will be suitable for B. globosus distribution in RCP4.5 and RCP8.5, respectively (Table 3).

In 2085 under the RCP4.5 and RCP8.5 climate scenarios, there will be contraction, shifts, and expansion in the habitat suitability of B. pfeifferi, as regions that were formerly suitable such as Ugu, eThekwini, iLembe, and King Cetshwayo districts may no longer be suitable. On the other hand, when there is extreme greenhouse gas emission (RCP8.5), some parts of King Cetshwayo district are likely to become suitable for the distribution of B. pfeifferi. Conversely, uThukela and uMkhanyakude districts, will always be suitable habitats regardless of the amount of greenhouse gas emissions. If the greenhouse gas emissions increase moderately (RCP4.5), districts such as uMgungundlovu, eastern part of Harry Gwala and parts of uThukela districts that are currently not suitable may become suitable habitats for intermediate host snails. In addition, under the RCP4.5 and RCP8.5 climate scenarios approximately 65.1% and 59.7% of the currently suitable areas will likely become unsuitable while 11.9% and 13.8% of the currently unsuitable habitats may become suitable for B. pfeifferi distribution in RCP4.5 and RCP8.5, respectively (Table 3).

The four key variables influencing the distribution of B. globosus are mean temperature warmest quarter, number of dry months, mean diurnal range in temperature, and isothermality (that is precipitation and temperature). Number of dry months, mean temperature warmest quarter, annual moisture index, and annual wettest quarter (summarily, precipitation and temperature) were the four key bioclimatic variables identified to influence the distribution of B. pfeifferi. Malacological surveys that investigated the influence of bioclimatic variables on B. globosus and B. pfeifferi have shown similar findings (Moodley et al., 2003; Stensgaard et al., 2013). The effect of temperature on the survival of snails have been well documented (Kalinda et al., 2017). The optimal temperature where growth and reproduction is most favourable for B. globosus is between 200°C and 30°C while that of B. pfeifferi ranges from 220°C to 28°C (Adekiya, 2018). At lower temperatures, metabolism slows down, reproduction rate reduces, and the snails become less active. At high temperatures, usually above 30°C, the snails become stressed which increases their susceptibility to diseases and eventual mortality (McCreesh et al., 2014). Thus, high temperatures in an environment will reduce the distribution of snails in that region (McCreesh and Booth, 2013). Since temperature affects vegetation growth, oxygen level, water availability, this in turn affects snail’s habitat suitability (Oso and Odaibo, 2021). Seasonal variation also affects the abundance of snail in an environment (Gurarie et al., 2017). Warm and humid summers have been found to be favourable condition for both B. globosus and B. pfeifferi. These warm temperatures also favour their lifecycle from egg laying, to hatching and to the development of larvae as the snails become more active (Nandy and Aditya, 2022). Similarly, it is expected that the water levels in dams, lakes, ponds, and reservoir will reduce due to increase in average surface temperature in the coming years (Kalinda et al., 2017). This will result in decrease in the population of B. pfeifferi. However, in regions and seasons where temperature is between 20°C and 25°C, it is expected that there will be an increase in the population of B. pfeifferi. However, at high temperatures and in areas where B. pfeifferi’s population is small, there is the possibility of them dying.

The amount of water present in an environment determines the abundance of snails. When there is no water, snails cannot survive and where there is too much water, the abundance of snail population reduces. Therefore, when there is rainfall, temporary snail habitats are created as well as new habitats due to heavy rainfall that moves snails from one environment to another (Manyangadze et al., 2021). This implies when there is heavy rainfall, the population density of snails in that environment reduces due to transportation. If the speed of flowing water exceeds 0.3 m/s, snails will be washed away to a new environment (Ofulla et al., 2013). Cold and dry seasons as well as post rainy seasons have been found to be suitable habitats for B. globosus and B. pfeifferi while hot dry seasons reduces the abundance of snails due to pools and rivers drying up. Therefore, long period of rainfall or absence of rainfall as well as the intensity of rainfall will determine the abundance of snails in a region (Rubaba et al., 2016). However, it is quite difficult to accurately determine the spatial relationship between rainfall, snail population dynamics, and infection transmission due to varying rainfall effect on species as well as geographical location (Brooker, 2007). Therefore, in areas where rainfall plays important role in snail population, population of snails reduces after each rainy season as well as their transport from one location to another, thereby creating temporary habitat (Manyangadze et al., 2021).

The observation that all the eastern part of the province currently has suitable habitats for B. globosus may be attributed to the average annual temperature in that part of the province. Temperature increases south-easterly in South Africa, hence KwaZulu natal province situated in the eastern part of the country is generally warmer with temperatures of about 25°C that create suitable habitats for B. globosus (Ebhota and Tabakov, 2021). Our findings corroborate the findings from a parasitology study carried out in KwaZulu-Natal as the districts with current suitable habitats for human schistosomiasis snail hosts were found endemic for schistosomiasis (Nemungadi et al., 2022). Western parts of uMkhanyakude district and some parts of uThukela district also shows current habitat suitability of B. globosus. As observed for B. globosus, suitable habitat for B. pfeifferi can also be found in the eastern part of the province except around iLembe district and eastern part of uMkhanyakude. Compared to B. globosus, the current habitat is more suitable for the distribution of B. pfeifferi. Other parts of the province are not suitable habitats for the distribution of both B. globosus and B. pfeifferi. However, some of these districts may have suitable habitats depending on the season. It is expected that during spring and summer months temperatures are high, the districts might become suitable for both B. globosus and B. pfeifferi. These suitable areas for both species are characterised by temperate or subtropical climates with warm to hot summers and mild winters. Hence, the warm and humid climate might be the factors influencing the distribution of both species as the habitats that are currently unsuitable for the distribution of B. globosus and B. pfeifferi such as Amajuba, Sisonke, and uMgungundlovu districts are characterized by very cold temperatures during the winter and high temperatures in the summers.

In KZN province, annual average temperature and annual average rainfall have been projected to increase for the period 2021–2050 under the RCP 8.5 scenario. Increase in the number of rainfall days may result in an increase in flood and storms. This projected increase will affect the environment, socio-economy, water, and food quality (Le Roux et al., 2017).

Temperature affects the growth, reproduction and population of snails which influences the rate of parasite and host development (Kalinda et al., 2017). Climate change predictions may play an important role in determining high-risk areas as increased temperature increases the development of hosts and parasites but limiting their survival. However, there have been conflicting reports on how increase in temperature will affect snail population. McCreesh, Nikulin (McCreesh et al., 2015) suggested that increase in temperature will result in increase in suitable habitats resulting in increased disease prevalence while Pedersen, Midzi (Pedersen et al., 2014) suggested that increased temperature will result in decreased suitable habitat leading to a reduction in the transmission of the disease. There is also the possibility of snails adapting to higher temperature which will increase habitat suitability and disease prevalence (Stensgaard et al., 2019).

These contrasting views may all be valid for certain areas. For regions where temperatures were too low and hence presenting unsuitable habitats, increased temperature due to climate change may result in the habitat becoming suitable and areas currently presenting suitable habitats may be affected by increased temperature resulting in the areas becoming unsuitable. This is supported by Kalinda, Chimbari (Kalinda et al., 2017) who reported that snails can survive at about 15.5°C during cold seasons but with increased temperatures, there will be high snail proliferation resulting in increased disease prevalence.

Figure 5 shows the predicted habitat suitability of B. globosus for both RCP 4.5 and RCP 8.5 indicating that habitats around the northern part of Amajuba district are always not suitable. However, with increased temperature, the unsuitable areas may drastically reduce as observed in Figure 5D. Also, with moderate increase (RCP 4.5) in temperature, most parts of Harry Gwala district are likely to remain unsuitable for B. globosus and this may be the situation in most parts of uMgungundlovu, iLembe, and uMzinyathi districts. Under extreme increased temperatures (RCP 8.5), most of these areas will become suitable habitats for B. globosus. According to the Koppen Geiger climate classification in South Africa, these areas are currently classified as always unsuitable habitat for B. globosus fall under the temperate, cold summer with no dry season (Cui et al., 2021). With increased temperature, it is expected that the cold summer will be warmer resulting in increased snails and increased disease prevalence. Most of uThukela district and parts of Zululand, uMzinyathi, uMgungundlovu, and uMkhanyakude that are currently not suitable will become suitable in the future. With extreme global warming, more parts in the districts will become suitable habitats for B. globosus. Only a small part of King Cetshwayo that is currently suitable will become an unsuitable habitat under RCP 4.5 while small parts of Ugu districts currently suitable will become unsuitable under extreme conditions. The eastern part of the province and the western part of uMkhanyakude as well as parts of uThukela district will always present suitable habitats for B. globosus. According to the Koppen Geiger classification, this region falls under the temperate hot summer with no dry season (Cui et al., 2021). Since there is no dry season, water is always present, and temperatures are conducive for survival of snails.

Habitat suitability of B. pfeifferi for both RCP 4.5 and RCP 8.5 show that Harry Gwala, Amajuba, Zululand, uMzinyathi, parts of King Cetshwayo, uMgungundlovu, and uMkhanyakude districts are always not suitable habitat for B. pfeifferi (Figure 6). This may be attributed to the weather condition in this region. However, with increased temperature, some of these regions like parts of Ugu, uThukela, and uMkhanyakude districts that currently not suitable will become suitable with moderate increase in temperature. However, under extreme cases, only parts of uMgungundlovu district and parts of Amajuba district that are currently unsuitable will become suitable. Much of the eastern part of KwaZulu natal province currently suitable may later become unsuitable. Since the climate classification classifies this area as temperate, hot summer with no dry season, with increased temperatures, the summer months can become very hot such that the existing snails will be desiccated leading to a reduction in disease prevalence. The modelled climate change has different effects on the habitat suitability of B. globosus and B. pfeifferi. There is likely to be expansion in the habitat suitability of B. globosus while there will be contraction and shifts in the case of B. pfeifferi. This may be because of the high tolerance of B. globosus to higher temperatures compared to B. pfeifferi (Joubert et al., 1986; Manyangadze et al., 2021). Parts of uThukela and uMkhanyakude districts were shown to always have suitable habitats for the distribution of B. globosus and B. pfeifferi. Therefore, priority should be given to these districts when prevention and control strategies are implemented to interrupt the transmission of schistosomiasis.

A limitation of our study is that we only investigated the influence of climatic variables on the habitat suitability of B. globosus and B. pfeifferi. As such, some important predictor variables such as topographic (altitude and slope), biotic (competition and predators), and behavioural variables which are important in determining the habitat suitability for B. globosus and B. pfeifferi distribution might have been omitted.