Study area

The study area was defined as the Chagas endemic region, which extends from northern Argentina and Chile to the southern United States of America. The study area for each individual species was defined as the area encompassing all reports of that species since the year 2000 plus a buffer zone of 5 degrees (approximately 300km).

Species occurrence and background points

The primary source of vector species data was a database of vector occurrence locations, which was supplemented with additional species presence points derived from a database of infections in vector species. Data on vector occurrence was extracted from DataTri, a publicly available database that reports the presence of a given triatomine species, the date of collection (if available) and geographical coordinates for each collection [32]. Additional vector occurrence data was added using a database of T. cruzi infections in triatomines that also provided the vector species found, the date of collection (if available) and geographical coordinates for each collection [33]. Any data points from DataTri that were duplicated in the second data set were removed before vector occurrence data from the infections database was added to the DataTri data set. Data points before the year 2000 were removed because the aim was to investigate vector distributions in the current era.

The available vector occurrence data is usually referred to as presence only data. Techniques for modelling such data often involve augmenting the presence data with pseudo-absence or background points, which requires a source of appropriate background data [34–36].

Here we use a target-group background (TGB) approach by choosing background data that exhibits similar sampling bias as the occurrence data [37]. This approach can reduce the bias introduced by preferential sampling of the presence locations. It was successfully used to map geographical distributions of malaria hosts and vectors [38] and predict infection risk zones of yellow fever [39]. In simulation studies this method also performed well when compared to approaches using presence-absence data [37, 40]. As with all models of presence-only data, the maps produced using the TGB approach represent relative rather than absolute probabilities of species occurrence.

We constructed one TGB dataset for each vector species as outlined below and illustrated in Fig 1 for Panstrongylus megistus:

1. The presence locations of vector species k = 1, …, K (target-group) were extracted from the database and a convex hull containing all presence locations was constructed (panel (A) in Fig 1).
2. This hull was extended by a constant width of 5 degrees in all directions to allow for uncertainty with respect to the range of the species being modelled (extended hull; panel (B) in Fig 1).
3. The presence locations of all other species within the extended hull were defined as background points (blue dots in panel (C) of Fig 1).
4. Duplicate observations at the same site and the same year were removed.
5. At the modelling stage, the background points were weighted such that their total weight is equal to the number of presence observations (cf. [37, 38]).

Download:

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

Fig 1. Construction of background points.

Illustration of the construction of background points using the TGB approach for species Panstrongylus megistus. Panel (A): A convex hull is constructed around the presence locations of the species. Panel (B): The hull is extended by a fixed width of 5 degrees (extended hull). Panel (C): Background points are added using presence locations of all other species within the extended hull. Panel (D): Blocks of width w_k are allocated randomly across the extended hull. The observations in blocks numbered 1-4 are used as training data, and the observations in blocks numbered 5 are assigned to the test data. One fold consists of all blocks sharing the same number. Figure created by the authors using R package tmap [41].

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

The weighting of background points means that if presence and background points were randomly distributed the predicted relative probability would be about 0.5. Thus probabilities > 0.5 indicate that it is likelier to observe presence than background, rather than an absolute probability of occurrence.

For some species, only few observations were available in the dataset. It was suggested that approximately five [42] or ten [43–45] events (presences) per predictor are required to reliably fit a logistic regression. Given that we use up to 30 predictors, this would imply a sample size of n ≈ 150 and ≈ 300, respectively. In our data, 14 and and 9 species fulfilled this (approximate) requirement. For completeness, we fit models for all species with n > 50, but obviously results must be interpreted with care as the sample size (number of presence observations) decreases (see Results section). Species with fewer than fifty observations in the training and test data were not modelled, however, their presence locations were used as background points for the species that were modelled.

Environmental variables

Previous work has shown that vector distributions are influenced by climate, land cover types and rural/urban classifications [19–24]. Environmental variables for these three data types were obtained at a resolution of 5 × 5 kilometres. The climatic variables used were land surface temperature (annual; day, night and diurnal difference) [46], two measures of surface moisture (annual) [47], rainfall (annual) [48], elevation (static) and slope (static) [49]. The variables used for land cover were the 16 IGBP land cover classes (annual) [50] and an enhanced vegetation index [51]. Finally the variables used to distinguish rural, peri-urban and urban areas were urban footprint (static) [52], nighttimelights (static) [53], human population (annual) [54] and accessibility (static, based on road networks and distance to cities) [55]). Annual environmental variables were not always available for all time periods in which occurrence data was available. In this case, we instead used values from the closest year. A full description of each variable is given in the supplement (Table A, S1 File).

Model evaluation

In the context of spatial analysis, data available for modelling often only encompasses few locations in areas for which predictions are generated. Therefore, the model is usually evaluated on out of sample data to avoid over-fitting, to ensure transferability to new locations and to obtain realistic estimates for the goodness of fit. Standard approaches to model evaluation, however, can yield over-optimistic metrics of the model predictive ability unless the spatial nature of the data (and the model) is taken into account [56–58]. To address these concerns, the data was initially split randomly into train-test data (80%) and evaluation data (20%), stratified by species. The latter dataset is not utilised during model building but later used to evaluate the models ability to interpolate and the final prediction. Additionally, the train-test data was split into five folds. Following recommendations in [59] each fold consisted of multiple spatial blocks, where the block size w_k for species k = 1, …, K was set such that approximately 50 blocks (10 per fold) would cover the extended hull of that species and defined as , where a_k is the area of the extended hull of species k. Fig 1 (panel (D)) depicts the resulting blocks and folds for species Panstrongylus megistus. For each species, blocks one through four were assigned to the training data, while blocks numbered five (grey shade) were only used to obtain out-of-sample test errors. Allocation of blocks was spatially random to avoid systematic bias of presence and background locations in any of the folds, but stratified with respect to species presence such that the proportion of presence and background points was approximately equal in all folds. The spatial blocking for all species considered in our analyses are provided in [60]. Model performance was evaluated by the area under the receiver operator curve (AUC), which measures the models ability to discriminate between presence and background points. After model evaluation as reported in the Results section was complete, the best model was refit on all data for the final prediction.

Modelling

To estimate the triatomine species distributions we fit a logistic regression to the target-group background (TGB) data using a generalised additive model (GAM). This modelling framework is comprised of two components: observations (1) and a linear predictor (2). We consider a Bernoulli process to model the background/presence (y_k,t,i ∈ {0, 1}) of each species k in year t ∈ {2000, …, 2016} at location . Within this framework, we specify the Bernoulli model (1) where i = 1, …, n_k are the observations per species. For each species k = 1, …, K, we set a spatial domain delimited by the extended hull formed by the spatial locations of the corresponding species (see Fig 1 for further detail). The relative probability of occurrence π_k,t,i is estimated by a logistic GAM with linear predictor (2) (2) where f_k,p(x_p,t,i) is the species specific, potentially non-linear, effect of the p-th covariate estimated by a penalised thin-plate spline [61] and GP_k,i(ℓ) is a two-dimensional, species-specific Gaussian process (GP) with range parameter ℓ evaluated at location s_i (the smoothness parameter was set to 1.5). The number of covariates P_k can vary by species as some of them might not have enough unique values within the spatial extent of the species to be relevant for analysis. Here, covariates were only included if the number of unique values was at least twenty. The correlation function of the GP was defined by C(x, x′) = ρ(||x − x′||), where ρ(d) = (1 + d/ℓ) exp(−d/ℓ) is the simplified Matérn correlation function with range parameter ℓ = max_ij||x_i − x_j|| as suggested in [62] and implemented in [63].

The model was estimated by optimising the penalised restricted maximum likelihood (REML) criterion (3) (3) using a double shrinkage approach where ψ is a vector of all coefficients associated with the smooth functions f and GP, D(ψ) is the model deviance and ϕ(⋅) and ϕ*(⋅) are range space and null space penalties of the model coefficients ψ [61, 64]. The first penalty (range space) controls the smoothness of functions f_k,p and GP_k, while the second penalty (null space) enables the removal of individual terms from the model entirely. The γ parameter can be used to globally increase the penalty and thus to obtain smoother, sparser and therefore potentially more robust models. Practical estimation was performed using techniques introduced in [65–67] to increase computational speed and reduce memory requirements.

Six model specifications (Table 1) were considered for this analysis, varying by the definition of covariate effects in Eq 2 and whether the global GP term GP_k,i(ℓ) was included. For each species, the final model (out of the six candidate models in Table 1) was selected based on its performance (AUC) on the test data (fold 5). Model 1 has no tuning parameters and was fit directly to the complete training data (folds 1-4). Models 2 through 6 were first tuned with respect to the global penalty γ ∈ {1, …, 4} based on 4-fold cross-validation on folds 1 through 4. Based on the value of γ that yielded the highest average AUC, the models were refit on the complete training data (blocks 1–4). These models were used to calculate the out-of-sample extrapolation and interpolation error (see “Results” section for details). After model evaluation the best model for each species was refit on all data (blocks 1 through 5 and the random hold-out data) to create final predictions.

Download:

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

Table 1. Model specifications considered in the analysis.

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

Implementation

All calculations were performed using the R language environment [68]. Thematic maps were created using package tmap [41]. Data munging and pre-processing was performed using packages dplyr [69] and tidyr [70]. Spatial cross-validation was set up using package blockCV [71]. Package mgcv was used to fit the GAMs [63].

Vectorial capacity of the mapped species

For each triatomine species that was mapped, information related to its importance in transmitting T. cruzi to humans was collated. The prevalence of infection with the T. cruzi parasite was calculated using the data from an existing repository [33]. Collections of less than twenty individuals of a species were excluded and the mean prevalence was calculated for all species where the remaining number of collections exceeded ten. Relevant behavioural data for each vector species was extracted from the published literature.

Species distributions

A total of 30 species were mapped. Triatoma infestans is predicted to occur from northern Argentina and Chile to southern Bolivia and Peru, overlapping in part with the predictions for Triatoma guaysayana, although T. guaysayana isn’t predicted in Peru. Panstrongylatus lutzi, Psammolestes tertius, Rhodnius nasutus, Rhodnius neglectus, Triatoma brasiliensis and Triatoma psuedomaculata are all predicted to primarily occur within Brazil. Triatoma sordida is predicted to occur in Brazil, Paraguay and Bolivia while Triatoma rubrovaria is mainly predicted to occur in Uruguay. Eratyrus mucronatus, Panstrongylatus geniculatus, Panstrongylatus rufotuberculatus, Rhodnius pictipes and Rhodnius robustus are predicted to overlap to differing degrees across a broad area that encompasses northern Bolivia and Peru, northwestern Brazil, Ecuador, Colombia, Venezuela, Guyana, Suriname and French Guiana. Triatoma maculata predictions are restricted to the northern part of this area, while Rhodnius prolixus is predicted even further north in Colombia and Venezuela, and Panstrongylatus chinai is only predicted at the far west of this area within Peru and Ecuador. The predicted distributions of Panstrongylatus geniculatus, Panstrongylatus rufotuberculatus, Rhodnius pallescens and Triatoma dimidiate all extend from northern South America to Central America. Triatoma barberi, Triatoma longipennis, Triatoma mazzotti, Triatoma Mexicana and Triatoma pallidipennis were all predicted to occur in Mexico only. Triatoma gerstaeckeri, Triatoma protracta and Triatoma rubida were predicted to occur from northern Mexico to the southern United States of America and Triatoma sanguisuga was predicted to occur exclusively in the southern United States of America.

A summary for all species that were modelled is provided in Table 2, including the specification of the model selected on training data as well as the AUC of this model evaluated on test data (fold 5) and the AUC obtained on the 20% randomly selected hold-out data (denoted by AUC*). The former is an indicator of the model’s transferability and ability to predict into new areas with potentially unseen covariate values or combinations, within the area that was modelled (cf. Fig 1). This is important because this is precisely the goal of the TGB approach and other modelling strategies that account for preferential sampling. The latter value indicates how well the model interpolates.

Download:

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

Table 2. Summary table for all species considered in the analysis ordered by number of presence observations.

https://doi.org/10.1371/journal.pntd.0008411.t002

The AUC values were all well above the 0.5 random classification threshold (mean: 0.85, SD: 0.12), indicating the maps usefulness to identify areas of higher probability of presence relative to background points. The comparatively low AUC values for species T. brasiliensis and T. pseudomaculata could be partially due to an overlap with many other species, thus making it difficult to discriminate between presence and background. The AUC* values were on average higher and had a lower variance (mean: 0.88, SD: 0.08) but generally consistent with the AUC values obtained on the spatial hold-out test data (fold 5). The predicted values (including 95% CI) in raster format for all species listed in Table 2 are given in [72], (.gri file format). Respective visualisations, i.e. map images, are available from [73].

Predicted distributions of four example species.

In this section, we present the results for four diverse species that provide examples from different genera including dominant and secondary vectors with domestic, peri-domestic and sylvatic habits. In addition, we provide detailed discussion of the results for one of the most important of these vectors, T. infestans, in the following section. Together distributions of these species cover most of the endemic region from northern Argentina and Chile to the south of the United States. Predictions are presented alongside bivariate maps that display prediction and (un)certainty in one map. To do so, predictions and uncertainty (defined by the width of confidence intervals) are divided into intervals (here [0, .25), [.25, .5), [.5, .75), [.75, 1], and [0, .075), [.075, .15), [.15, .3), [3, 1], respectively. The cut off for uncertainty was chosen because a CI width of ≥ .3 means that the upper and lower CIs fall into different categories of the probability intervals. The legend in the bivariate maps indicates which colours correspond to which combination of prediction and uncertainty.

Triatoma gerstaeckeri colonises homes and kennels and is known to bite humans (S2 File). The predicted distribution of this species from the southern USA through much of Mexico is shown in Fig 2a. Uncertainty in these predictions is high at the northern boundaries of this species where sampling is particularly sparse (Fig 2b and [60]). Triatoma dimidiata colonises homes, as well as sylvatic and peri-domestic habitats (S2 File), and is a dominant vector of T. cruzi. It’s distribution predicted using data for the time period from 2000 onwards ranges from southern Mexico through Central America into northern Venezuela and Colombia (Fig 2c and 2d). Predictions are high with low uncertainty in Central America whereas predictions are lower with higher uncertainty in northern Venezuela and Colombia. It is important to note that the block width for this species is smaller than an estimation of the range of spatial auto-correlation, thus the AUC values may be optimistic. Rhodnius pictipes colonises palm trees and has been implicated in the contamination of products for human consumption (S2 File). The predicted distribution of R. pictipes in northern Brazil and Bolivia, French Guiana, Suriname, Guyana, Venezuela, Colombia, Ecuador and Peru is shown in Fig 2e. It is important to note the areas of high uncertainty within the region of higher predicted probability of presence for this sylvatic species (Fig 2f). The confidence intervals for this species are higher than those seen for Panstrongylus geniculatus in the same region, which reflects the lower volume of data available for R. pictipes. Panstrongylus geniculatus colonises trees and rodent nests, has been found in urban areas, and is known to bite humans (S2 File). Its predicted distribution, which overlaps with that of R. pictipes in South America and extends further north as far as Honduras in Central America, is shown in Fig 2g. The uncertainty in these predictions is typically low but is higher at the fringes of the area where the predicted relative probability of presence is high (Fig 2h).

Download:

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

Fig 2. Predicted relative probability of occurrence.

Predictions for 4 selected species at a resolution of 5 × 5 km within the respective extended hull of species occurrence (left panel) and bivariate map (right panel), where darker colours indicate higher predicted probabilities while the transition from white/pink to turquoise/blue indicates increased uncertainty. Row 1 (A, B): Triatoma gerstaeckeri; row 2 (C, D): Triatoma dimidiata; row 3 (E, F): Rhodnius pictipes; row 4 (G, H): Panstrongylus geniculatus. Figure created by the authors using R package tmap [41].

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

Distribution of Triatoma infestans.

Arguably, the most important T. cruzi vector species is T. infestans making it a key target for indoor residual spraying campaigns that have the potential to alter the distribution of this predominantly domestic species. Intervention coverage data was not available to our models so it is particularly important to consider the uncertainty in the predictions for T. infestans. Fig 3 shows the predicted probabilities for this vector species using data from the year 2000 onwards (left panel) alongside a bivariate map that highlights the (un)certainty of the estimation (right panel). The model predicts areas of high relative probabilities of presence with higher certainty in southern Bolivia, northern Chile and northwestern Argentina, which aligns well with the observed presence points. The uncertainty is usually high in unsampled areas, e.g., along the border of Argentina and Chile, and eastern parts of Paraguay. Low probabilities of occurrence are predicted in Brazil and northeastern Paraguay, with varying levels of certainty.

Download:

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

Fig 3.

Left panel: Final predicted map for T.infestans. Right panel: A bivariate map of the predictions that indicates areas of high vs. low probabilities together with the model uncertainty. Darker colours indicate higher predicted probability. Transitions from white/pink to turquoise/blue indicate higher uncertainty.

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

Ranking importance of the environmental variables

Caution is needed when drawing conclusions from information on which environmental variables were selected by each model because i) many of these variables are highly correlated, for example temperature, rainfall, surface wetness and elevation, and ii) some important variables may not have been made available to the models as discussed above. It is, however, interesting to note the three variables selected as most important by each species model and these rankings are given in Table D, S3 File. Forest cover was one of the three most important variables for three species and Eratyrus mucronatus, Panstrongylus rufotuberculatus and Rhodnius pictipes are all known to colonise trees (S2 File). Eight other species models selected other vegetation cover variables as important and variables that define the urban-rural gradient (urbanicity, accessibility and human population) were selected as being among the most important by seven species models. Looking across all species models (Table E and Figure A, S3 File), the most important variables for predicting T. cruzi vector species distributions were temperature, the Gaussian process (the spatial component), evergreen broadleaf forest cover, the vegetation index, rainfall and elevation, followed by variables defining the urban-rural gradient and other types of vegetation cover. Nine land cover classes were never selected by any model. Unsurprisingly these were water, needleleaf forest cover, built up areas (information that was provided to the model by other variables that were selected), snow and ice, barren areas and unclassified land (which is a rare occurrence in the land cover data). For most species there are rarely single covariates with a contribution of more than 50%, meaning that each prediction is comprised of smaller contributions from many variables (Figure A, S3 File).

Vectorial capacity of the mapped species

The current state of knowledge on factors related to the capacity of each of the mapped species to transmit T. cruzi to humans (vectorial capacity) is summarised in Table 3. Specifically, mean infection prevalence, confirmation of human blood meals in natural vector populations, and confirmation of colonisation or invasion of homes (including in urban areas) are listed. Less information is available on the feeding-defecation interval or defecation location for each species, and these values may be influenced by the different experimental conditions used, so these variables are not included in Table 3 but sources of evidence are listed in the supplement (Table B, S2 File). The 30 most commonly reported species mapped here encompass the five most important dominant vectors that frequently colonise homes (P. megistus, R. prolixus, T. brasiliensis, T. dimidiata and T. infestans) as well as species that often colonise peridomestic habitats such as chicken coops, rats nests, boundary walls, wood piles, palm trees, and livestock housing. These species encompass a range of mean T. cruzi infection prevalences from 0.8% in T. sordida to 55.6% in T. longipennis, although for 13 of the most commonly reported species there was insufficient data to generate a reliable mean infection prevalence value. In addition to differences among the mean values for species, there is also considerable variation within each species likely resulting from heterogeneities in factors such as intervention deployment and local host species. When viewing the summaries in Table 3, it is important to note that not all regions or species have been sampled or tested equally and a lack of published evidence for a specific component of vectorial capacity cannot be taken as definitive evidence of its absence. For example, no infections have been reported in Eratyrus mucronatus or Psammolestes tertius, but only 28 and 143 individuals have been tested, respectively, compared to 335,467 T. sordida individuals. In addition to the five important dominant vectors, there is evidence that many of the species mapped in this work are potential vectors of T. cruzi. Almost all of the 28 species that have been found to be infected with T. cruzi are known to invade homes, and at least 17 have been found to have fed on humans (Table 3). The sources of evidence—130 published articles in total—are given in the supplement (Tables B and C, S2 File).

Download:

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

Table 3. Infection prevalence and behaviour of selected species.

https://doi.org/10.1371/journal.pntd.0008411.t003