Trypanosoma cruzi reservoir—triatomine vector co-occurrence networks reveal meta-community effects by synanthropic mammals on geographic dispersal

Vector and mammal reservoir data

Georeferenced data points for triatomines and terrestrial mammal collection sites in Mexico were used to model interactions. For triatomines, we compiled a database consisting of all documented data collections since 1979 (publication of the first global monograph of Triatominae by Lent and Wygodzynsky which re-assigned many autochthonous Mexican species complexes), published in technical reports (Mexican National Health Secretary), scientific publications, or scientific collections (i.e., INSP, INDrE and IBUNAM), that includes 1600 data points for 26 Mexican triatomine species (Ramsey et al., 2015; doi: 10.5061/dryad.rq120). Mammal data points were retrieved from the Global Biodiversity Information Facility (GBIF; http://www.gbif.org), registered for Mexican collections (last accessed in February of 2011). We obtained 36,000 records for 396 mammal species, which were all used for network construction.

Interaction model based on vector-mammal co-distributions

The methodology for the vector-host interaction inference we will use is that of (Stephens et al., 2009) which we describe here for completeness. Triatomine species hosts are a subset of all terrestrial mammal species present in their distribution range I′⊆I. Let the difference of the probabilities P(B_i|I″) − P(B_i) be the probability of the presence of a particular triatomine species (B_i) occurring given the presence of a particular mammal species (I′). To calculate this probability, we used a grid size N of 0.25 ° for the cell size. The robustness of the results with respect to changes in the grid size were considered in Stephens et al. (2009). The grid covers the area of all ecoregions where the species Bi has been reported. We consider this probability relative to the null hypothesis that B_i and I″ are uncorrelated and hence P(B_i|I″) = P(B_i). Then, $P\left(B_i|\mathbf{I}″\right)-P\left(B_i\right)=\left(N_{Bi\text{and}{\mathbf{I}}^{\prime }}∕N_{\mathrm{I}}^{\prime }\right)-\left(N_{Bi}∕N\right)$, where N_{Bi and I′} is the number of spatial cells where there is a co-occurrence of the taxon B_i and the taxon I′ (the potential mammal host), and $N_I^{\prime }$ is the number of cells where that host take their stated values (a Boolean presence/absence value), within the grid. P(B_i|I_k) is the number of co-occurrences of the taxa B_i and I_k, and allows us to define the most important statistical associations between these pairs of species distributions. However, since P(B_i|I_k) is a probability, it does not consider sample size. If P(B_i|I_k) = 1, this may be the result of a coincidence of B_i and I_k in one or 1,000 spatial cells, the latter being more significant statistically. To adjust for this difference, the following binomial test for statistical significance is used: (1)${\displaystyle \epsilon \left(B_i|I_k\right)=\frac{N_{I_j}\left(P\left(B_i|I_k\right)-P\left(B_i\right)\right)}{{\left(N_{I_j}P\left(B_i\right)\left(1-P\left(B_i\right)\right)\right)}^{1∕2}}.}$ Equation (1) measures the statistical dependence of B_i on I_k relative to the null hypothesis that the distribution of B_i is independent of I_k and randomly distributed over the grid. The sampling distribution of the null hypothesis is binomial, where every cell is given a probability P(B_i) of having a point collection of B_i. The numerator of Eq. (1), is the difference between the observed number of co-occurrences of B_i and I_k, and the expected number in the case where the distribution of point collections was obtained from a binomial with sampling probability P(B_i). Since this is a stochastic sampling, the numerator must be measured in appropriate “units”. As the underlying null hypothesis is that of a binomial distribution, it is natural to measure the numerator in standard deviations of this distribution, which is used as the denominator of Eq. (1). If N_Bj ≥ 5–10, and P(B_i) and P(B_i|I_k) are not too close to zero or one, then a normal approximation for the binomial distribution is adequate, in which case ε(B_i|I_k) = 1.96 represents the standard 95% confidence interval. Naturally, as ε values increase, the hypothesis of interaction becomes stronger. Note that such a statistical association does not necessarily prove that there is a direct “causal” interaction between the two taxa, rather it allows for a statistical inference that should be validated experimentally subsequently or assessed with independent studies from the literature.

If co-occurrence patterns between triatomines and mammals are indicative of a vector-host interaction, more significant ε values should be found for those species’ pairs that are confirmed blood sources of triatomines and/or mammal reservoirs of T. cruzi. As a model performance measure, we examine the ε value for known interactions identified empirically and for those reported in the literature (see Table S1). Usually, these reports target specific triatomine species (i.e., the blood source of Triatoma pallidipennis as reported in Mota et al. (2007), etc), therefore, model performance was evaluated using those species for which evidence is available. If the vector-host interaction is an important component of T. cruzi transmission, the ε value should correlate with the probability for a mammal species to be a T. cruzi reservoir. To evaluate this, we listed all confirmed reservoir mammals and conducted the same analysis as described above.

Description of regional availability of vector hosts and interaction patterns

Vector host availability was measured as the number of mammal species in the vector’s distribution range (the grid). This number determines the potential for vector–host interactions and the opportunity for vector specialization. A simple statistical description of frequency distribution of host availability was conducted across vector species to analyze how vectors relate to mammal species richness. The same procedure was conducted to analyze the number of mammal species with which each vector had an epsilon value ≥1.96. We extracted the mean and standard error from the group of mammals with which each vector had significant interactions to describe comparatively the strength of vector-mammal interactions.

Validation of the interaction model

A literature search was conducted in Google Scholar for publications that report the presence of T. cruzi in Mexican wild mammals and in vector hosts using bloodmeal analysis. We used these records to evaluate the model performance by comparing the distribution frequency of historically confirmed mammal hosts grouped by epsilon range against the distribution frequency of all mammals with no reports of T. cruzi presence. This allowed us to describe the model’s predictability across different epsilon values.

Network construction

We used 26 Triatominae species and 396 mammal species (Table S1) for network vertices (nodes). Vertices were inter-connected representing significant triatomine-mammal geographic associations (ε ≥ 1.96). All species were classified using binary variables associated with their ecological tolerance for anthropogenic habitats (synanthropic vs. non-synanthropic) and with T. cruzi infection, the latter used for model validation (Table S1). To examine the geographic pattern of vector-mammal potential interactions, a network was projected onto a map of Mexico in which nodes were located according to their range centroid, which was determined as the geographic centroid calculated from the point collections of the node species. We used an automatic procedure to draw the network with a force-directed algorithm that assigns forces among the set of edges and nodes. Briefly, this algorithm emulates a system as if the edges were springs and the nodes were electrically charged particles. Forces are applied to the nodes, thereby pulling them closer together or pushing them further apart. This is repeated iteratively until the system comes to an equilibrium state i.e., their relative positions do not change from one iteration to the next, and the graph is then drawn. The equilibrium state represents all forces in mechanical equilibrium. In the final graph, the position of every node corresponds to their connectivity level, indicating that the most connected vertices are located in the center of the graph and vertices with fewer edges are in the periphery (Kamada & Kawai, 1989). In order to analyze geographic patterns of vector-mammal interactions relevant for Chagas disease risk, we used modular algorithms to identify network “community structure”, which use centrality indices to identify community boundaries. Instead of constructing a measure for edges that are central to communities (as in hierarchical clustering), focus was placed on the least centralized edges that are “between” communities (Newman, 2006). Rather than constructing communities by adding the strongest edges to an empty vertex set, they were constructed by progressively removing edges from the original graph. Vertex betweenness was used as the centrality measure, since it reflects the node influence in the network. The betweenness centrality of a vertex i is defined as the number of shortest paths between pairs of other vertices that run through i. To find which edges in a network are between other pairs, we generalize betweenness centrality to edges and define the edge betweenness as the number of shortest paths between pairs of vertices. If there is more than one shortest path between a pair of vertices, each path is given equal weight so that the total weight of all of the paths is unity. If a network contains communities or groups that are only loosely connected by a few inter-group edges, then all shortest paths between different communities must attach to one of these edges. Thus, the edges connecting communities will have high edge betweenness. By removing these edges, we separate groups and reveal the underlying community structure of the graph (Girvan & Newman, 2002). We further analyzed the structure of the network by disentangling the relative effect of triatomines and mammals grouped according to their disturbance tolerance (synanthropic and non-synanthropic) and parasite association (T. cruzi and non-T. cruzi), in the network structure. To do this, we analyzed the mean degree of species belonging to each group. The degree of a node is the number of connections (edges) it has to other nodes in the network. Finally, we estimated the vector species’ similarity in terms of their mammal host communities using a hierarchical cluster analysis based on the 1-Jaccard similarity coefficient (Legendre, Borcard & Peres-Neto, 2005). All analyses were conducted using the igraph package ver 1.0 for the statistical package R ver 2.15 (http://www.r-project.org). All routines were included in a script written to perform analyses (Supplemental Information 1).

Regional availability of vector hosts and interaction patterns

As would be expected, mammal species’ number per vector was moderately correlated with vector grid size, i.e., the number of spatial grid cells occupied by a given vector (Spearman rank correlation coefficient r = 0.61, P = 0.001). The number of mammals to which vectors had significant epsilon values was less correlated with vector grid size (r = 0.5, P = 0.01). Six species of the phyllosoma complex (Triatoma pallidipennis, Triatoma longipennis, Triatoma mazzottii, Triatoma picturata, Triatoma mexicana, and Triatoma gerstaeckeri), one of the protracta (Triatoma barberi), and two of the dimidiata complex (Triatoma dimidiata Pacific hg2, and T. dimidiata Gulf hg3), had their distributional ranges which coincided significantly with regions having a high number of mammal species (i.e., higher than the 95% CI [79–128]). In contrast, eleven triatomine species (Triatoma protracta, Triatoma peninsularis, Triatoma sinaloensis, Triatoma neotomae, Triatoma lecticularia, Triatoma nahuatlae, Triatoma uhleri, Triatoma zacatecensis, Triatoma woodi, Dipetatolagster maximus and Paratriatoma hirsuta) had distribution ranges in regions with a significantly lower number of mammal species (i.e., lower than the 95% CI). These latter species, conversely, had a higher proportion of significant links (>95% CI) with mammal species with highest mean epsilon values (Figs. 1A and 1B). This pattern was inverted in vectors with high availability of mammals, since they had a lower proportion of accessible mammal hosts with significant links (Fig. 1B).

Rodentia and Chiroptera were the mammal orders with greatest number of species, however orders with the highest proportion of synanthropic species were Didelphimorphia and Rodentia, followed by Carnivora (Table 1). Didelphimorphia had the highest proportion of registered T. cruzi infected species, followed by Carnivora. Although Primates, followed by Rodentia, had greatest link strength (mean epsilon), rodents, opossums (Didelphimorphia), and bats (Chiroptera) had highest epsilon values. Opossums and Lagomorphs had the greatest number of links with vector species, although no Lagomorphs have been reported infected with the parasite.

Vector-host co-occurrence and Trypanosoma cruzi transmission

All mammal species known to be T. cruzi hosts had significant epsilon values with at least one vector. The distribution frequency of host mammals was skewed positively with the magnitude of epsilon, in contrast to those mammals which have not been reported infected with T. cruzi. There was a significant difference between the proportion of T. cruzi hosts with high epsilon values, and the proportion of non-hosts with high epsilon (Fig. 2).

Network topology and the importance of synanthropic hosts

The interactions wiring patterns differed depending on the spatial position of the geographic centroid of vector and host ranges (Fig. 3). Edge density is higher in the Neotropical region, in central and southeast Mexico, as compared to the Nearctic region in the north. Synanthropic species (both triatomine and mammal) shape the Kamada–Kawai force-directed network structure, as indicated by their position in the center of graphs, while most non-synanthropic species are located on the periphery (Fig. 4). The structure of the complete array of mammals and vectors that can be considered relevant for the T. cruzi dispersal provided their geographic associations were not random, includes all significant links (ε ≥ 1.96) (Fig. 4). Mammals which are both synanthropic and T. cruzi-infected are the most central in the network. These synanthropic mammals remained in greater proportion (compared to non-synanthropic mammals) in networks with higher epsilon values (ε ≫ 2). In general, this last group of mammals (synanthropic and reported reservoirs), had proportionally more significant links than the group of registered reservoir species that are non synanthropic (OR =1.17, 95% CI [1.01–1.4]). At highest epsilon values (ε ≥ 8), the network connectivity was reduced to isolated clusters (Fig. 4), although these clusters had several nodes of synanthropic/T. cruzi reservoir species. For instance, Triatoma dimidiata Yucatan hg1 (node 365) had 21 links in the network with ε ≥ 10, and 9 (43%) of those links were synanthropic and T. cruzi reservoirs.

Network connectivity, as measured by node degree, decreased exponentially with an epsilon increase, both in vectors (Fig. 5A) and mammals (Fig. 5B). Synanthropic vectors were more connected than non-synanthropic ones (Fig. 5A). Similarly, synanthropic reservoirs of T. cruzi were more connected than other groups, particularly at low yet significant epsilon values (Fig. 5B). All groups had the same pattern of connectivity decrease with increasing epsilon, except for non-synanthropic reservoirs.

Two principal clusters emerged from the complete geographic model network. These clusters differed according to number of species and geographic region. The first cluster was composed principally of phyllosoma (seven species) and dimidiata (two haplogroups) complex species and 171 mammals, while the second cluster contained principally protracta and rubida complexes (17 vector species and subspecies) and had 183 mammals. The more dominant species of this latter group are T. rubida, T. uhleri and T. woodii. Several mammal and vector species were intermediate between these two former groups and could not be assigned to either: T. dimidiata hg1, T. barberi, T. longipennis, T. mexicana, and T. pallidipennis (Fig. 6). Hierarchical similarity of vectors constructed from mammal hosts produced two main clusters. Most synanthropic vectors were grouped into one cluster, while most non-synanthropic species were separated into a second cluster (Fig. 7).