Introduction

Tsetse flies (Diptera: Glossinidae) are considered the sole cyclical vectors of African trypanosomes and are reported to occur in about 10 million km² in sub-Saharan Africa [1]. The trypanosome parasites cause Human African Trypanosomosis (HAT) or sleeping sickness in humans and African Animal Trypanosomosis (AAT) or nagana in livestock. Both diseases have a substantial negative effect on agricultural development and economic growth in sub-Saharan Africa [2,3]. HAT is a fatal disease if left untreated and, although absent in southern Africa, occurs regularly in some regions of sub-Saharan Africa with 70 million people at risk of becoming infected in 36 countries [4].

Tsetse flies are restricted to sub-Saharan Africa [5] and they have been sampled as far south as the north-eastern parts of KwaZulu-Natal Province (KZN) of South Africa (Latitude S28°31’13.44”). This southern population extends into the neighbouring Maputo Province (MP) of Mozambique [6–11]. Of the 31 described tsetse fly species and subspecies, only two species are found as far south as South Africa, i.e., Glossina brevipalpis Newstead belonging to the Fusca (forest) species group and Glossina austeni Newstead belonging to the Morsitans (savannah) species group [1,7,8].

In addition to several species of wildlife, both species feed on cattle [5] and are involved in the transmission of two pathogenic protozoa, Trypanosoma congolense and Trypanosoma vivax [12–14] that cause the debilitating disease AAT in livestock. T. congolense is the most abundant species in South Africa [13] and limited vector competence studies indicated that G. austeni was the more competent vector for T. congolense in this area [14,15].

Tsetse infested areas in KZN are mainly used for communal farming and inhabited by 426 000 humans, 130 000 small ruminants and 360 000 cattle [16]. Livestock production and agricultural development is severely hindered by these flies as vectors of trypanosomes which causes considerable stress to the farmers, not only in KZN but in the entire distribution area of the two species [14]. The G. brevipalpis belt extends from Ethiopia in north-eastern Africa southwards to Somalia, Uganda, Kenya, Rwanda, Burundi and Tanzania [10]. In southern Africa the presence of G. brevipalpis extends from Malawi, Zambia, Zimbabwe, Mozambique to the north-eastern parts of KZN [5]. G. austeni is found in East Africa from Somalia in the north, extending southwards into Kenya, Tanzania, Zimbabwe, Eswatini, Mozambique and the north-eastern parts of KZN [5,10]. G. austeni was also present on Unguja Island of Zanzibar, Tanzania, but the population was eradicated in 1997, after implementing an area-wide integrated pest management (AW-IPM) campaign that included the sterile insect technique (SIT) [17].

The most southern distribution of these two species is shared by Mozambique (Matutuíne District), Eswatini (Mlawula Nature Reserve) and South Africa (north-eastern KZN). Localised surveys and available distribution prediction models for South Africa showed that, on a micro-ecological scale, G. brevipalpis and G. austeni are confined to pockets of dense vegetation in north-eastern KZN [18–21]. This apparent patchy distribution, in addition to possible low migration potential [7,18], may have resulted in the development of localised genetically isolated population pockets [22]. Preliminary studies, using molecular and phenetic (geometric morphometrics) markers, however, suggested an absence of significant barriers to gene flow within the pockets of this southern population [23]. These data provided evidence that G. brevipalpis and G. austeni populations in southern Mozambique, Eswatini and South Africa can be considered homogenous and that an eradication strategy that is limited to one country will not be sustainable because of potential reinvasion from uncontrolled neighbouring areas [23]. However, this southern tsetse population might be geographically isolated from the main tsetse belt that starts approximately 500 km north of Matutuíne District, south of the Save River in central Mozambique [6–9,11,24]. If confirmed, this will offer an opportunity to sustainably create a tsetse-free zone in southern Africa.

The current strategy in South Africa to manage nagana relies on long-term vector suppression using live-bait technologies and the ad hoc treatment of cattle with trypanocidal drugs. For the live bait technology, the chemical amitraz (C₁₉H₂₃N₃) that has normally been used in the existing network of dipping tanks for the management of tick-borne diseases, was replaced with a wide-spectrum insecticide, i.e., the pyrethroid cyhalothrin (C₂₃H₁₉ClF₃NO₃) that is more effective against dipteran flies [25]. This cattle-dipping regime for suppressing the tsetse populations has been adopted since 2015 (personal communication Dr. L. Ntantiso, Department of Agriculture, KwaZulu-Natal, South Africa). Although pyrethroids as a dipping agent are effective for the management of both ticks and tsetse flies [26], ticks can develop pyrethroid resistance with a reduction in its efficiency as an acaricide as previously reported in the area [26,27]. This livestock dipping strategy can also only be effective if a large proportion of the tsetse population feeds on domestic rather than wild animals [28]. Managing AAT with pyrethroids in dipping tanks can therefore only be a temporary solution at best. The establishment of a tsetse-free zone in southern Africa would be a more sustainable solution to the nagana problem. A potential sustainable solution to eradicate G. brevipalpis and G. austeni from north-eastern KZN was proposed in 2007 [7] and was based on an AW-IPM strategy that includes an SIT component. In the absence of significant barriers to gene flow between the populations of MP, Eswatini, and KZN [23], the proposed strategy [7] will have to be adapted to target sequentially the populations in all three countries.

Determining the precise and potential geographic distribution and abundance of the targeted tsetse population will be vital for the success of an AW-IPM strategy. The size of the area that needs to be treated will directly affect the outcome, sustainability and cost of any control campaign. Because of the logistical and funding constraints involved in implementing surveys over such a large rural area that extends over three countries, models that predict the distribution, presence and even abundance of the tsetse populations would facilitate the development and implementation of these surveys and make them much more focussed and cost-effective [29].

The first probability of presence model developed for South Africa was based on tsetse fly sampling data collected with sticky XT traps between 1993 and 1999 [19,30]. This model incorporated climate and environmental variables and predicted a more extensive geographical distribution for both G. brevipalpis and G. austeni than was indicated by the sticky trap data [19], suggesting that the model may have overestimated the distribution [30]. Although this model was updated and refined by including data from tsetse sampled with the more effective H-trap between 2005 and 2007, there were still areas where the predicted probability of presence could be improved. Incorporating data on cattle density, human population, agricultural intensity and detailed vegetation biomes were suggested to improve model fit [20]. Furthermore, the earlier models did not include southern Mozambique and Eswatini. The objective of the study was to develop a MaxEnt species distribution model and habitat suitability maps for the southern tsetse belt of South Africa, southern Mozambique and Eswatini using available entomological survey data of G. austeni and G. brevipalpis. The updated models may be able to predict to what extent changes in land use and agricultural practises may potentially influence tsetse abundance and the occurrence of AAT.

Materials and methods

Study area

In South Africa the envisaged tsetse infested area (16 000 km²) stretches from around 10 km south of the Mfolozi River in the south, for approximately 190 km, to the border of Mozambique in the north, and from the Indian Ocean coast in the east for 80 km to the Hluhluwe-iMfolozi Park in the west (Fig 1) [7,8]. The infested area extends into the Matutuine District (8 500 km²) (MP) of Mozambique, the northern limit being the Boane and Namaacha Districts of Mozambique [6,11]. In the east, it borders with the Indian Ocean and in the west with Eswatini. In Eswatini, the Mlawula Nature Reserve was surveyed. The reserve is located west of the Lebombo Mountains (elevation 776 m), an 800 km-long narrow range of mountains that stretch from Hluhluwe in KZN in the south to Punda Maria in the Limpopo Province in South Africa in the north parallel with the Mozambique border.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Apparent density of Glossina austeni and Glossina brevipalpis collected between 2009 and 2019 from Maputo Province (MP), Mozambique, Eswatini and north-eastern KwaZulu-Natal Province, South Africa (https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/PA7U7L).

https://doi.org/10.1371/journal.pntd.0009989.g001

Farming systems are predominantly subsistence farming with numerous communal farms interspersed with several protected areas consisting of provincial and private game parks and reserves. This includes the iSimangaliso Wetland Park which was listed as South Africa’s first World Heritage Site in December 1999. These areas contain a wide variety of game animals consisting of large numbers of bigger mammals, small primates, rodents as well as birds that are potential hosts for tsetse. The South African target area contains several state forests, mostly pine and eucalyptus plantations, and sugarcane farms.

The area has a large range of land cover relevant to the presence or absence of G. brevipalpis and G. austeni such as coastal forests, savannah and agricultural areas. The climate is subtropical, except the mountainous “highveld” area in Eswatini that has a temperate climate.

Data analysis

Tsetse occurrence and density.

The data collected during these independent entomological surveys, were correlated with vegetation classes developed for the area and incorporated in models to predict suitable habitat for potential tsetse distribution. The statistical software R version 3.6.2. using RStudio Desktop version 1.2.5033 [32] was used for data analysis. A co-occurrence analysis, consisting of a binary presence-absence matrix, was carried out to assess if the distribution of G. austeni and G. brevipalpis throughout the study area was segregated, aggregated or random. A Checkerboard score (C-score) [28] that is measuring aggregation or segregation intensity (checkerboardness) was calculated. To evaluate the tsetse relative abundance, a one-way analysis of variance (ANOVA) was used to compare the mean tsetse fly AD between sites. The data were not normally distributed and nonparametric Kruskal-Wallis tests were used as a post test. Additionally, Dunn’s multiple comparison tests were used where P value < 0.05.

Spatial filtering and generation of pseudo-absence.

All datasets were spatially filtered to reduce the spatial correlation in the modelling process. One record was kept per pixel (500 m) and any absence data within a buffer of 2 km around a presence data was removed from the analysis. A probabilistic model was applied to remove non-significant absence from the analysis (5% level) [33,34]. Absence data were kept for validation and some background (pseudo-absence) was generated randomly in the grid according to the kernel based on sampling effort. While higher sampling effort increases the pseudo-absence generated, areas less sampled will have less. The kernel was built using spatial point pattern analysis incorporating all available trap information from the study area. Subsequently, one pseudo-absence per raster cell (approx. 500 m x 500 m) was selected. Finally, all pseudo-absences into a 2 km radius buffer around presence data were removed, taking into consideration the average flight range of these species. To reduce the sampling bias these pseudo-absence data were generated by taking into account the environment and the range of the efficiency of the H-traps to sample both species. A model based on a multidimensional nonparametric kernel [35] was developed in areas of high sampling effort to correct for sampling bias. The models were fitted with G. austeni and G. brevipalpis presence data collected with H-traps from the three countries. Additionally, the models were validated using presence and absence data for both species.

Remote sensing data.

Time series of high spatial resolution (100 m to 1000 m) remote sensing data were downloaded, cleaned and summarized to build relevant covariates. Ten years of Moderate-resolution Imaging Spectroradiometer (MODIS) products (January 2010 to December 2019) were used.

A detailed vegetation classification map of the tsetse distribution area in South Africa, Eswatini and southern Mozambique was derived from four Landsat 8 images (resolution 30*30m, June to July 2013) and field surveys conducted in July 2013. The field records comprised 1092 photographs taken at locations within the known tsetse distribution area of South Africa, Eswatini and southern Mozambique. The photographs were used to identify 11 main types of vegetation based on the international nomenclature that were considered relevant for the presence or absence of G. brevipalpis and G. austeni (Fig 2), i.e., savannah woodland, herbaceous savannah, shrub savannah, dense dry forest, gallery forest, tree plantations, crops (agricultural areas), urban areas, swamps, water bodies and bare ground were the mainland-cover classes. MODIS products and indices were selected to capture the complexity of tsetse habitat preference, in particular, thermal and vegetation indices. MODIS products were acquired from the NASA Earth Observing System data server, the human population density layer from the Worldpop project [34], cattle density from the FAO Database [36], FAO livestock gridded of the world project [37], and layers of protected areas from the World Database on Protected Areas [38].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Vegetation classes map with Glossina austeni and Glossina brevipalpis presence and absence collected between 2009 and 2019 (https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/PA7U7L).

https://doi.org/10.1371/journal.pntd.0009989.g002

For each variable, data were pre-processed and cleaned by re-interpolating the data at a spatial resolution of 500 m using the nearest neighbour method. A MODIS QA (pixel quality) mask was applied to remove poor quality pixels. Thermal data were filtered using the boxplot algorithm to reduce the effect of outliers [39]. Summary statistics such as mean, maximum, minimum and range were computed for the variables Normalized Difference Vegetation Index (NDVI), Middle Infra-Red (MIR), Enhanced Vegetation Index (EVI), Tree Cover (treecov), Day Land Surface Temperatures (DLST) and Night Land Surface Temperatures (NLST) (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Variables derived from remote sensing data used in the model.

https://doi.org/10.1371/journal.pntd.0009989.t001

Model selection and validation.

Ecological Niche Factor Analysis (ENFA), a variant of factor analysis, was used to explore and model species ecological niche [40]. The environmental space used by the species was compared with the available environmental space using two indicators i.e., marginality and specialization. Marginality was used to measure niche central position. It captures the dimension in the ecological space in which the average conditions where the species accrues differs from the global conditions. A large marginality value implies that the conditions where the species was found were “far” from the global environmental conditions. In contrast, specialization measures the spread and usage of the ecological space along dimensions of niche use. The higher this value, the narrower the space used by the species. Consequently, the species niche can be summarized by an index for marginality and specialization and represented on a factor map within the biplot framework.

The Maximum Entropy ecological niche model (MaxEnt), one of most common species distribution models [41], was used to predict the potential distribution of the two species. It uses a machine learning method based on the information theory concept of maximum entropy [42]. MaxEnt fits a species distribution by contrasting the environmental condition where the species is present, and the environment characterized by some pseudo-absence data also called background. The logistic output from this method is a suitability index that ranges between 0 (low suitable habitat) and 1 (high suitable habitat). This output was used to create suitability maps for both species.

Predictors used in these models were chosen according to the ecology of the two species involved. Twenty-four bioclimatic variables based on remote sensing data were built to model the distribution of these two species. For each of the remote sensing derived index, we computed some measures of position (mean, min, max) and some measures of spread (range, coefficient of variability). Multiple models were fitted varying the background of the model to assess the effect of pseudo-absence and use model averaging [43]. The effect of pseudo-absences was assessed by repeating the process to generate them several times and monitor variability and stability in model quality metrics. Furthermore, the predictive quality of the different MaxEnt models was assessed with a subset of the entomological survey data that was excluded from the analysis. Uncertainty was measured using the coefficient of variation from the predictions. The main metric used to check predictive power of the different models was the Area Under the ROC Curve (AUC). A likelihood-based metrics was also analysed to assess the importance and contribution of each variable.

Based on the predictive quality of the developed MaxEnt models, and keeping the degree of uncertainty in mind, the potential distribution of the two species was determined by extrapolation of the model to adjoining areas outside the collection limits.

Results

Maximum entropy ecological niche model (MaxEnt)

Results of the ENFA showed that G. austeni occurrence was positively associated with vegetation indices (EVI range (5.1%) and tree cover mean (2.2%)) as well as presence of protected areas (Fig 3). The same positive association with vegetation indices and presence of protected areas was observed for G. brevipalpis (Fig 3). Human population and cattle densities and most of the temperature indices showed a negative association with occurrence of both species, except the minimum night land surface temperature which was positively correlated with occurrence of G. brevipalpis. Slope and aspect were positively related to G. brevipalpis presence but negatively with G. austeni. Mean night land surface temperature influenced the habitat for the two species. Minimum MIR for G. austeni and mean MIR for G. brevipalpis accounted for most of the variance and fell outside the cloud of average conditions available in the study area (Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The ecological niche factor analysis plan.

Light grey polygon shows the overall environmental conditions available in the study area, dark grey polygon shows environmental conditions where Glossina austeni and Glossina brevipalpis were observed, and the white circles corresponds to the barycentre of its distribution.

https://doi.org/10.1371/journal.pntd.0009989.g003

For G. austeni, the variable cattle density contributed 24.8% towards model fitting, followed by the range of day land surface temperatures (16.7%) and presence of protected areas (13.9%). For G. brevipalpis, the variable cattle density contributed 42.6% towards model fitting, followed by mean night land surface temperature (12.9%) and presence of protected areas (8.9%) (Fig 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Contribution of variables to the suitability index by decreasing importance for Glossina austeni and Glossina brevipalpis.

The 95% confidence interval is indicated in red and individual values in grey.

https://doi.org/10.1371/journal.pntd.0009989.g004

For the Maxent models the mean AUC as assessed with the independent data was 0.88 (range 0.87 to 0.89) and 0.84 (range 0.81 to 0.85) for G. austeni and G. brevipalpis, respectively (Fig 5). The uncertainty grid for the habitat suitability index model (Mass analyses) indicated that the uncertainty in the predictions was low, except in areas west of the known tsetse infested zone (S1 Fig). When the model was extended northwards to the Gaza and Inhambane Provinces of Mozambique the uncertainty increased. However, the prediction values along the coastline in the Inhambane Province retained relative high certainty. The distribution of the two species was similar with high suitability areas in and along the protected areas in all three countries and the coastal areas of north-eastern KZN and MP.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

Area under the curve for the average MaxEnt model (in red) and the 10 sub models (in grey) for Glossina austeni and Glossina brevipalpis.

https://doi.org/10.1371/journal.pntd.0009989.g005

Suitable habitat for G. austeni was predicted mainly along the coastline and inside protected areas of north-eastern KZN and MP (Fig 6). The predicted suitable habitat extended to a small degree into the communal farming areas close to the protected areas and coastline. In Eswatini the suitable predicated habitat was also linked to protected areas along the eastern border of Mozambique. There was also a small suitable area in the central part of Eswatini.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. The mean habitat suitability index predicted by a MaxEnt model for Glossina austeni and Glossina brevipalpis for Maputo Province (MP), Mozambique, Eswatini and north-eastern KwaZulu-Natal Province, South Africa (https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/PA7U7L).

https://doi.org/10.1371/journal.pntd.0009989.g006

Overall G. brevipalpis displayed a wider distribution than G. austeni (Fig 1). The suitable habitat for G. brevipalpis was highly linked to protected areas but extended to a greater degree into adjoining communal cattle farming areas as compared with G. austeni. This extension was more pronounced in the MP than in KZN. Suitable habitat for G. brevipalpis was also predicted over a narrow strip at the border between Eswatini and Mozambique, and, like G. austeni, in the central part of Eswatini, but to a larger extent. In north-eastern KZN, a suitable band was predicted in the northern part in the communal farming area along the coast and border of the iSimangaliso Wetland Park. The largest area of predicted suitable habitat for G. brevipalpis in the communal farming area was west and south of Hluhluwe-imfolozi Park. The area of high probability for G. brevipalpis extended south along the coast far beyond its known historical distribution limits.

The model was extended for both species (Fig 7) to include the Gaza and Inhambane Provinces of Mozambique. The model indicated that large parts of this area were unsuitable for both species, but the habitat was more suitable for G. brevipalpis than for G. austeni. Areas of high suitability were mostly found along the coast and in a small area of the Zinave National Park close to the Save River.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Extended mean habitat suitability index predicted by a MaxEnt model for Glossina austeni and Glossina brevipalpis in Mozambique, Eswatini and South Africa (https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/PA7U7L).

https://doi.org/10.1371/journal.pntd.0009989.g007

Discussion

A sound understanding of the potential distribution of targeted tsetse fly species will be essential for the successful implementation of any AW-IPM programme that include an SIT component, e.g., it will establish the area over which sterilized males need to be released in programmes. The development of habitat suitability maps, and an assessment of the factors that regulate the presence or absence of the targeted species, will therefore be advantageous, if not essential, to ensure that these programmes will be efficient and cost-effective. Reliable prediction maps will enable more efficient planning of entomological surveys and selection of appropriate trapping sites. The present study used existing entomological survey data to develop habitat suitability maps for G. austeni and G. brevipalpis using statistical methods whilst incorporating what is known of the ecology of both species, e.g., host and habitat preferences. The methodology followed to develop these prediction maps was similar to the one used in the AW eradication programme of Glossina palpalis gambiensis Vanderplank in the Niayes of Senegal [29,41]. The MaxEnt approach was used as previous studies in West Africa [41] indicated that it better predicted suitable landscapes and tsetse presence as compared with a presence-absence regularized logistic regression model. The MaxEnt model was develop based on ADs as obtained with odour baited H-traps for both species. Although the H-trap was specifically developed for the collection of these species [44] the differential efficacy of this trap remains unknown for the two species. In order to develop better informed and robust models, systematic sampling within study area is still the most efficient way to improve this model. In agreement with previous entomological surveys in the affected area in north-eastern KZN, climate, vegetation and presence of protected areas and cattle were the key regulators that determined the presence and abundance of G. brevipalpis and G. austeni [8,16,18,21].

Previous surveys and distribution prediction models showed that the G. austeni and G. brevipalpis populations were, on a micro-ecological scale, mainly confined to pockets of dense vegetation in north-eastern KZN [18,19,21,30]. A strong relationship between tsetse relative abundance and vegetation type was observed in the current study. High relative abundance of G. austeni and G. brevipalpis was strongly associated with savannah woodland and dense dry forests, respectively. Low relative abundance of G. austeni and G. brevipalpis was associated with swamp and crop areas, respectively. This association was evidenced by the positive association between tree cover and predicted suitability for both species. These observations highlight the regulatory role that the presence of suitable vegetation may play in predicting the presence or absence of these two tsetse fly species in an area and the success of bush clearing as an earlier method for tsetse control [45].

The model predicted suitable habitat for both species along the coast of north-eastern KZN and MP (Fig 6). Previous studies conducted in north-eastern KZN suggested a relationship between the relative abundance of the tsetse populations and the range of temperature variation [8]. In accordance with model predications the relative abundance of both species was higher at coastal sites compared to sites in the interior, and this can be related to the variation in the average temperature and relative humidity, which are less pronounced at the coast as compared with the interior. The range in day land surface temperatures contributed as much as 16.7% towards model fitting for G. austeni. For G. brevipalpis the contribution of the range in day land surface temperature (1.6%) was of less importance than the mean night land surface temperatures (12.9%). This relationship was also reflected in the prediction model, as there was a negative association between probability of presence and an increase in the range of day surface temperatures.

The model indicated that the predicted presence of both species is associated with the presence of protected game areas in all three countries. Protected areas contributed 13.9% and 8.9% towards model fitting for G. austeni and G. brevipalpis, respectively. These protected areas harbour the preferred hosts of G. brevipalpis such as hippopotamus (Hippopotamus amphibious), African elephant (Loxodonta africana) and African buffalo (Syncerus caffer) as well as small game animals such as bush pig (Potamochoerus larvatus) and duikers (Sylvicapra species) that are the preferred hosts of G. austeni [5,46,47]. In general, these protected areas are characterized by denser vegetation and lusher tree cover as compared to the communal farming areas.

Notwithstanding the apparent host preferences of these two species, both will feed on cattle [5] and cattle will probably be able to sustain tsetse populations in the absence of game in areas with suitable vegetation and climatic conditions. In the current study, cattle density contributed the highest percentage towards model fitting for both G. austeni and G. brevipalpis. For both species, cattle densities, irrespective of suitable habitat, showed a negative association with fly occurrence. This apparent contradictive observation may be related to cattle being used as live bait for tsetse control in north-eastern KZN. In the past the existing extensive dipping network, that was established mainly for tick control, was modified after an outbreak of nagana to include tsetse fly control by replacing the acaricide dipping agent with an insecticide [26,27]. Once the outbreak was under control the dipping agent was changed back to the acaricide to prevent ticks from developing insecticide resistance. The current live-bait control program in north-eastern KZN started in 2015 (personal communication Dr. L. Ntantiso, KZN Department of Agriculture and Rural Development) and the impact of this control method on tsetse abundance is reflected by the predicted low suitability at diptanks in the communal farming areas where control is currently implemented.

The current control strategy in north-eastern KZN may have facilitated the creation of tsetse refuges in the protected areas. The perceived positive relationship between protected areas and tsetse occurrence and relative abundance were confirmed by the MaxEnt model. And this emphasizes a strong reinvasion potential from the uncontrolled protected areas into the communal farming areas. Over the last decade game farms, private nature reserves and other forms of wildlife-oriented land use, have increasingly become prominent features in KZN [48]. This expansion of protected areas combined with the absence or at least a general reluctance to suppress tsetse populations in these protected areas [49], increases the probability of migration into the farming areas with the associated constant threat of nagana transmission to livestock. This threat is exacerbated by the presence of buffalo and other wildlife hosts that are considered a reservoir host of nagana in some of these protected areas [12].

Previous probability of presence models of G. austeni and G. brevipalpis [20] looked at the estimated area covered by each species in north-eastern KZN at three probability threshold cut-offs, P > 0.5, P > 0.25 and P > 0.125. The total potentially infested area increased two-fold when the threshold was decreased from P > 0.5 to P > 0.125, from 5600 km² to 11 750 km² [20]. The current developed model included cattle density, human population, agricultural intensity and detailed vegetation biomes and the prediction area was extended to include Eswatini and southern Mozambique. The estimated area covered by both species at a probability threshold of P > 0.5 was 1700 km², 57.4 km² and 3901 km² for north-eastern KZN, Eswatini and MP, respectively. Therefore, the current model predicted a smaller tsetse fly-infested area than previous models. The reason for this is not clear but it can be speculated that this may be due to model refinement, changes in the habitat (e.g., bush clearing for agricultural purposes) as well as the apparent success of the current tsetse control actions in the area. Considering the sensitivity of tsetse to environmental factors this can even be indictive of the effects of climate change in the area. A more intensive study will be needed to pinpoint the exact reason for the predicted smaller tsetse fly-infested area.

The generated prediction model has been expanded to include the Gaza and Inhambane Provinces in Mozambique, an area suspected to be free of tsetse following a reduction in cattle and wildlife numbers due to human settlement and population expansion. This assumed tsetse-free area extends northwards for approximately 500 km from the northern part of the Matutuine tsetse belt up to the great Central and North tsetse belts that start at the Save River [9,11]. In agreement with historical data [10] and in view of a reduction in cattle and wild host populations in the area, the present model predicted very few areas with suitable habitat in the Gaza and Inhambane provinces of Mozambique (Fig 7). Suitable areas, with a low habitat suitability index, for G. brevipalpis (<0.6) and G. austeni (<0.2) were restricted to the coast (Fig 7). Factors that limit the presence of tsetse flies in the area may be the lack of suitable vegetation combined with a greater variation in temperature range in the interior. Historical collection data from 1984 for G. austeni, indicated that they were present 60 km west of Maputo and found to be present in the Androstachys forests [50,51].

The exact distribution of G. austeni and G. brevipalpis in Mozambique is not known, but it is suspected to be discontinuous due to natural barriers and land use changes. No tsetse control programmes have been implemented since the late 1960s and only limited surveys have been implemented in the last 20 years [6,9,11,52]. This emphasises the importance of the developed MaxEnt model as it provides indications of the most suitable areas that can be prioritised for the deployment of traps during surveys. This will significantly reduce time and cost of future surveys in this area.

Assessment of the MaxEnt model quality using an independent data set attained good performance with a high predictive power (AUC > 0.80 for both species). The available trap data however only partly validated this prediction model, as the present model seems to indicate a potential wider distribution of the two tsetse species as compared to the survey data. This underpins the importance of these models as tools not only for the planning of the surveys and the monitoring activities, but also for the suppression and later sterile male release activities.

The factors that contributed or shaped the distribution of G. austeni and G. brevipalpis as revealed by the present prediction model, indicated that changes in the climate, agricultural practises and land use can have a significant and rapid impact on tsetse presence/absence and abundance. The current trypanosomosis control strategy in South Africa (dipping with insecticides) can only be effective if a large proportion of the tsetse populations feed on cattle [28]. The observed trend in increased wildlife areas in KZN will lead to a decrease in the proportion of tsetse flies taking a blood meal from treated cattle. To manage trypanosomosis it will be crucial to supersede the current South African control strategy of using cattle as live baits with a more sustainable AW-IPM strategy. In 2007, an AW-IPM strategy that include an SIT component was proposed to establish a tsetse fly free South Africa [7]. The proposed AW-IPM strategy suggested the division of the infested area into four zones from south to north with the successive implementation of four phases (pre-suppression, suppression, SIT and post-eradication) in each zone following the rolling carpet principle [7]. In respecting the principles of an AW-IPM approach, the control effort should be directed against the entire insect population, and therefore, the proposed AW-IPM strategy of 2007 should be modified to include the tsetse populations from southern Mozambique and Eswatini.

Supporting information

Acknowledgments

We thank our co-workers in the Epidemiology, Parasites & Vectors Programme as well as the staff at the KwaZulu-Natal Tsetse Research Station for their valuable assistance.