Gradients in predation risk in a tropical river system.

The importance of predation risk as a key driver of evolutionary change is exemplified by the Northern Range in Trinidad, where research on guppies living in multiple parallel streams has provided invaluable insights into the process of evolution by natural selection. Although Trinidadian guppies are now a textbook example of evolution in action, studies have generally categorized predation as a dichotomous variable, representing high or low risk. Yet, ecologists appreciate that community structure and the attendant predation risk vary substantially over space and time. Here, we use data from a longitudinal study of fish assemblages at 16 different sites in the Northern Range to quantify temporal and spatial variation in predation risk. Specifically we ask: 1) Is there evidence for a gradient in predation risk? 2) Does the ranking of sites (by risk) change with the definition of the predator community (in terms of species composition and abundance currency), and 3) Are site rankings consistent over time? We find compelling evidence that sites lie along a continuum of risk. However, site rankings along this gradient depend on how predation is quantified in terms of the species considered to be predators and the abundance currency is used. Nonetheless, for a given categorization and currency, rankings are relatively consistent over time. Our study suggests that consideration of predation gradients will lead to a more nuanced understanding of the role of predation risk in behavioral and evolutionary ecology. It also emphasizes the need to justify and report the definition of predation risk being used.

Rivers are defined by gradients, with abiotic conditions changing from their upper to lower reaches in a predictable manner. Abiotic gradients are easily quantified and widely used to help predict and explain the ecology of freshwater communities (Vannote et al. 1980). However, biological elements, such as species richness and predation pressure, are also nonrandom features of the system that can be usefully employed to explore other ecological and evolutionary patterns (Beecher et al. 1988; Gilliam et al. 1993; Reimchen 1994; Jackson et al. 2001). In parallel with the physical conditions of the continuum, these biological gradients run from upstream stretches, which typically have lower species richness and lower predation pressure, to downstream stretches where species richness and predation pressure are almost always higher (Matthews 1998). Explanations for these gradients include dispersal limitation due to waterfall barriers, preventing some species from colonizing upstream regions, as well as habitat characteristics such as water depth, which may be too shallow in higher reaches to support larger species (Harvey and Stewart 1991).

Predation pressure can have a profound effect on the structure and composition of aquatic communities (Matthews 1998; Jackson et al. 2001). The Northern Range of Trinidad has played an important role in helping us understand the role of predation in natural ecosystems. Identified as a “natural laboratory” in the 1940s (Haskins et al. 1961), over the past 60 years, studies on guppy Poecilia reticulata populations living in Trinidad’s multiple parallel streams have provided invaluable insights into evolutionary ecology (e.g., Liley and Seghers 1975; Reznick et al. 1990; Magurran et al. 1992). Indeed, some of the best evidence for evolution by natural selection has emerged from seminatural experiments using this system. Most notably, this work includes some of the first experimental demonstrations of evolution on the timescale of years (rather than millennia) and in a vertebrate (rather than bacteria in a petri dish) (e.g., Endler 1980; Reznick et al. 1997; O'Steen et al. 2002).

Gradient versus Dichotomy

One of the reasons the Northern Range has been so fruitful in terms of evolutionary research is that some of its parallel streams are interrupted by waterfalls, which can be tens of meters high and thus prevent the upstream dispersal of certain fish species. Guppies have successfully colonized many of the above-barrier sites, and these are considered as “low risk” from the guppy’s perspective, as typically only the killifish (Anablepsoides hartii, formerly called Rivulus hartii) also exists here. In contrast, those below waterfalls are regarded as “high risk” because guppies living there find themselves cohabiting with larger predatory species, including the pike cichlid Crenicichla frenata (formerly called C. alta) and the wolf-fish or “guabine” Hoplias malabaricus. The advantages of this “barrier waterfall” phenomenon to evolutionary study was first noted in the 1950s (Haskins and Haskins 1951) and since then the “high” versus “low” dichotomy has helped scientists link differences in predation risk to differences in a huge range of traits, including size, behavior (Liley and Seghers 1975), coloration (Endler 1980) and numerous life history traits (Reznick and Endler 1982).

Traditionally, guppy research has emphasized the contrast between the extremes in this system. However, as most researchers are aware, there are finer scale differences in guppy traits along gradients of predation pressure in these rivers (Gilliam et al. 1993; Torres Dowdall et al. 2012). To understand how subtle differences in local conditions influence evolutionary outcomes we first need to be able to properly quantify the expected predation gradient. We also need to be aware that using different categories of predation pressure potentially influences the ranking of sites.

Predation Categorization

Most studies into guppy adaptation to predation do not quantify the predation pressure as such, but instead characterize different sites according to which species are found there (presence/absence), either historically, or at the time of sampling. Endler (1978) classified 6 “types” of predator assemblage. However, the 2 species that most researchers pay attention to are the wolf fish H. malabaricus and the pike cichlid C. frenata (Farr 1975).

A limitation of this approach is that we know that many other species also predate on guppies, even if to a lesser degree (Seghers 1973). Even the cohabitant most closely associated with “low”-predation sites, A. hartii, is known to consume juvenile guppies (Mattingly and Butler 1994). A key issue, therefore, is deciding which species to consider as potential predators. Despite many decades of research on this system, guppy predators have yet to be systematically ranked in terms of their actual risk to a guppy. Fortunately, we do have some information on trophic level and feeding mode which allows us to assign the fish occurring in an assemblage to categories based on the degree of risk they likely pose (Kenny 1995; Phillip et al. 2013).

Even once justification has been made for which species “count” as predators, the next question is should these be quantified in terms of their numerical abundance or their biomass? Different abundance currencies provide different insights in community ecology (Magurran 2004). In general, biomass is more commonly used in aquatic systems, and numerical abundance in terrestrial studies. Numerical abundance is widely used to explore community dynamics, whereas biomass may be a better indicator of resource allocation (White et al. 2007). Although some studies treat these abundance currencies as interchangeable (Taper and Marquet 1996), this is not necessarily the case in taxa, such as fish, where body size varies considerably amongst species and through ontogeny (Mittelbach et al. 1988). It is important to ask, therefore, whether our perception of where sites lie on a gradient of risk depends on the abundance currency being used.

Temporal Dynamics

A final consideration is that aquatic assemblages are not static, and any survey to assess predator numerical abundance or biomass is just a snapshot of a dynamic system. Indeed, temporal turnover is a feature of all ecological communities (Magurran and Henderson 2010). Thus far, relatively little attention has been paid to the temporal dynamics of the Northern Range system (Magurran 2005). For example, we know that both guppy density (Reznick 1989; Grether et al. 2001; Magurran 2005) and sex ratio (Pettersson et al. 2004) vary considerably over time, yet we understand very little about how these, or other patterns, may be related to temporal changes in the dynamics of predatory species.

Aims

Using data from 16 sites across 8 different rivers along the southern slopes of the Northern Range collected over a 6-year period, we seek to evaluate the potential of the biological gradient of predation pressure in the Northern Range as an alternative to the more commonly applied high–low predation dichotomy. We categorize predation pressure into 3 levels—narrow, intermediate, and broad, using published assessments of the diets of Trinidadian fish species (See Table 1). We use this categorization to, first, examine the evidence for a gradient in predation risk in this system. We do this by plotting site trends in predation pressure for each category and for 2 abundance currencies (numerical abundance and biomass). Next, we ask if site rankings shift across predator categories within currency, and across currency within predation categories. Finally, we assess the degree of temporal consistency in site rankings when predator category and currency are held constant.

Table 1.

Species of predator included in each category

Scientific name (and authority)	Common name	Narrow	Intermediate	Broad
C. frenata (Gill 1858)	Pike cichlid	X	X	X
H. malabaricus (Bloch 1794)	Wolf fish/Guabine	X	X	X
A. hartii (Boulenger 1890)	Hart’s killifish		X	X
Gymnotus carapo (Linneaus 1758)	Cutlass fish/Knifefish		X	X
Agonostomus monticola (Bancroft 1834)	Mountain mullet			X
A. pulcher (Gill 1858)	Blue acara/Blue coscorob			X
Cichlasoma taenia (Bennett, 1831)	Brown acara/Brown coscorob			X
Corynopoma riisei (Gill 1858)	Swordtail sardine			X
Hemigrammus unilineatus (Gill 1858)	Featherfin tetra			X
Rhamdia quelen (Quoy and Gaimard 1824)	River catfish/Silver catfish			X

Note that C. frenata was formerly called C. alta; A. hartii was formerly called R. hartii, and Andinoacara pulcher was formerly called Aequidens pulcher. “Narrow” includes only the 2 most strongly piscivorous species, “intermediate” includes an additional 2 moderately piscivorous species, and the broad category includes all piscivorous fish in the system.

Materials and Methods

Data collection

A total of 16 sites in Trinidad’s Northern Range were sampled once a year (July/August) for 6 years (2011–2016). These sites are spread between 8 different rivers along the southern slopes of the Northern Range, and each consists of a 50 meter stretch of stream (Figure 1) .

Figure 1.

Map showing the location of the 16 sites across Trinidad’s Northern Range. Site abbreviations are as follows: A1: Acono site 1; A2: Acono site 2; C1: Caura site 1; C2: Caura site 2; L1: Lopinot site 1; L2: Lopinot site 2; LA1: Lower Aripo site 1; LA2: Lower Aripo site 2; M1: Maracas site 1; M2: Maracas site 2; Q1: Quare site 1; Q2: Quare site 2; T1: Turure site 1; T2: Turure site 2; UA1: Upper Aripo site 1; UA2: Upper Aripo site 2.

Sampling involved blocking the stretch at either end using seines, before fishing with a 2-person seine net, followed by electrofishing. Dip nets were also used in the shallows. All fish were identified to species, and weighed individually to the nearest gram, with the exception of guppies, which were too small to be weighed in the field, and were instead noted as mature female, mature male, or juvenile (typical weights were then used to estimate biomass for this species). On occasions where fish were seen but not caught, length was estimated and an appropriate length–weight curve consulted to give an estimate of biomass (see Deacon et al. 2017 for more details).

Analyses

We used 3 alternative categories of predation, based on trophic classifications using available information on diet (Kenny 1995; Phillip et al. 2013). The first takes a “narrow” view of what a guppy predator is, including only 2 “strongly piscivorous” species widely acknowledged to be the main predators of guppies: H.malabaricus and C.frenata. The second “intermediate” view includes an additional 2 species that are considered “moderately piscivorous.” Finally, the “broad” view encompasses even those fish classified as “weakly piscivorous” (Table 1).

Assuming that predation risk depends both on the presence and abundance of predators, the predation pressure of each site was calculated for 6 different scenarios: narrow, intermediate, and broad predation categories, and for the 2 abundance currencies—numerical abundance and biomass. Numerical abundance was defined as the total number of individuals of the species in a given category, summed over the year’s sampling. Biomass was the total wet weight (g) of all individuals of the species in a given category. All analyses are based on these annual abundance data, and took place in R 3.3.2 (R Core Team 2016).

We first constructed box plots for each of the 6 scenarios based on mean annual predation values for each site. We used log2 numerical abundance and log2 biomass, respectively, to depict predator abundance, and positioned sites along the x-axis in order of predator abundance. The goal of these plots was to assess if the sites lie along a gradient or represent a dichotomy; the former would be represented by a continuum of sites, ranging from high to low, and the latter by distinct groupings of sites. The plots also allowed us to ask if some categories provide stronger evidence for a gradient than others.

Next we constructed tanglegrams [using the plot() function in R] to visualize the consistency in site ranking by predator category and abundance currency, and quantified this consistency using a Spearman’s rank correlation test to compare the rankings. Sites were ordered in terms of decreasing mean predator abundance over the entire 6-year sampling period, with the sites with the highest values at the top of the tanglegram.

Finally, using linear models we plotted the trend (predator abundance against year) for the 16 sites using the 6 approaches (3 categories using numerical abundance and 3 categories using biomass). We evaluated the temporal consistency of site rankings within each predator category and currency using a Kendall concordance test using [the kendall.global function in the R package vegan (Oksanen et al. 2007)].

Results

Gradient versus dichotomy

Plots of log numerical abundance and biomass illustrate how the 16 sites capture the gradient of predation risk over space (Figure 2). The narrow category of predation separates out the 2 lowest predation sites (UA1 and UA2; see Figure 1) whereas all other sites fall on a continuous gradient.

Figure 2.

Mean values, standard errors (boxes), and ranges (whiskers) for 6 time points for each site, for log predator numerical abundance (upper row) and log predator biomass (lower row) for the 3 categories of predators (narrow, intermediate, and broad). Sites are plotted in rank order along the X-axis. See Figure 1 for key to the site codes.

Predation categorization

When quantifying predation using numerical abundance of individuals, site rankings shift across the predation categories. Specifically, the ranking of sites differs greatly depending on whether the narrow category of predators is used as opposed to the intermediate. There is no significant correlation between the rankings of narrow and intermediate numerical abundances (r = 0.06; P = 0.824) or between narrow and broad (r = 0.12; P = 0.660). Note the switch from lowest to highest ranked for the Upper Aripo sites between narrow and intermediate rankings. In contrast, the intermediate and broad categories rank the sites similarly and are highly correlated (r = 0.96; P = 0.001) (Figure 3A).

Figure 3.

Tanglegrams comparing rankings of sites by mean predation risk, depending on whether a narrow, intermediate, or broad category of predator assemblages is applied and whether (A) numerical abundance and (B) biomass is used as currency. See Figure 1 for key to the site codes.

The ranking of sites is similar across predation classifications if biomass is used to quantify predator abundance, with significant correlations in all cases (narrow + intermediate: r = 0.52, P = 0.039; intermediate + broad: r = 0.56, P =  0.025; narrow + broad: r = 0.70, P =  0.002) (Figure 3B).

Sites are ranked consistently irrespective of abundance currency for the narrow category of predation (r = 0.65; P = 0.007) (Figure 4A). In contrast, when the intermediate or broad categories of predation are used, there is no significant correlation of site rankings between numerical abundance and biomass (Intermediate: r = 0.21; P = 0.443; Broad: r = −0.11; P = 0.680) (Figure 4B,C).

Figure 4.

Tanglegram comparing how sites are ranked depending on whether mean biomass or numerical abundance are used to quantify predation, using the (A) narrow, (B) intermediate, and (C) broad categories of predation. See Figure 1 for key to the site codes.

Temporal dynamics

Considerable variation in predator abundance exists over time within sites for all categories and currencies of predation. Different sites show different trends, but in general the site rankings remain consistent over time regardless of the approach to quantifying predation pressure (Figure 5). This is reflected in significant Kendall concordance scores for rankings over time for all approaches, with the exception of when biomass and the intermediate category of predation are used to rank the sites. Numerical abundance does a marginally better job at ranking the sites consistently over time than biomass, with greater concordance coefficients for all categories.

Figure 5.

Trends over time for sites based on ranks per year for log2 numerical abundance and log2 biomass, as calculated according to the 3 categories of predation. Each trend line represents a site. Kendall concordance coefficients (W) and P values are given for each plot.

Discussion

We find compelling evidence that, in terms of risk to guppy populations, river sites across the Northern Range of Trinidad lie along a continuum of predation risk. We further show that site rank along this continuum is broadly maintained over time. However, we also find that where a particular site lies on this gradient depends on how the predator assemblage is defined.

Despite the expectation that biotic gradients (including predation intensity) parallel abiotic ones in river systems (Vannote et al. 1980), few studies have utilized the gradient concept when asking questions about evolutionary ecology in streams (although see Jourdan et al. 2016). Dichotomous contrasts, pitching high risk against low risk, have been widely employed. Many of the most famous examples are from Trinidad’s Northern Range streams, a disproportionate number of which represent huge strides in our understanding of evolutionary ecology (e.g. Endler, 1980; Reznick and Endler, 1982; Magurran et al. 1995, to name just a few). This approach has the advantage of highlighting the stark differences in behavior, physiology, and ecology associated with habitats that offer sharp contrasts in risk.

After 5 productive decades focusing on the “high–low” dichotomy, a study in 2012 confirmed that the Guanapo River represents a useful gradient of predation pressure for investigations of guppy life history evolution (Torres Dowdall et al. 2012). The nature of predator communities for the Guanapo was defined in the early 1990s, when Gilliam et al. (1993) mapped the presence of different fish assemblages along its length. Torres Dowdall et al. (2012) linked this biological gradient to a fine scale pattern of variation in terms of guppy life history traits, moving beyond the traditional dichotomy to add nuance to our understanding of the effects of predation risk on guppy evolution. Our results show that this gradient is apparent between as well as within stream systems. As such it emphasizes the need for careful and consistent quantification of predation pressure.

Our data show that sites do not form discrete clusters, but rather fall along a continuum of predation risk. This is true irrespective of the abundance currency or predator category used. Temporal variation means that sites exhibit overlap in predation pressure over the time frame of this study. The one exception occurs for the narrow predation category (regardless of currency): due to the barrier waterfalls some sections of stream (in this case 2 Upper Aripo sites) never support either of the main guppy predators, therefore these extreme low-predation sites are separated from other sites. Numerical abundance appears a better discriminator of sites along this gradient than biomass (see Figure 2).

Although our sites were not selected a priori to represent the dichotomy within each river, the Upper and Lower Aripo sites (UA1, UA2, LA1, LA2) do comprise a typical and well-used dichotomous pair. The Upper Aripo is a frequently utilized “low predation” river, whereas the Lower Aripo is a commonly sampled “high predation” river (e.g. Magurran and Seghers, 1990; Botham et al. 2006). Using the narrow classification of predation (for either abundance currency), the Upper Aripo sites are consistently separated from the Lower Aripo sites at extremes of the overall gradient (Figure 2). However, when the intermediate or broad categories are employed this pattern disappears. Reasons for this difference are discussed in the following section.

Perhaps the most important finding in this study is that the manner in which predation is quantified determines where a site falls on this predation risk gradient. In other words, the ranking of a given site can shift markedly, depending on which species are considered predators, and the abundance currency being used. When the predator assemblage is expanded to include species that potentially include guppies in their diet, but are less likely to be specialized predators, sites can move from being perceived as very low risk, to having relatively high risk.

In our case, broadening the definition of predation to include 2 more species, turns the Upper Aripo sites (previously “no predation” sites) into the 2 sites ranked as highest in predation risk. This is primarily due to the inclusion of the killifish A. hartii in the “intermediate” category. Although adult A. hartii will consume juvenile guppies (Mattingly and Butler 1994; Fraser and Lamphere 2013) alongside invertebrates and allochtonous material, they are themselves prey to many of the other predators. Thus, in the absence of larger predators A. hartii are typically found at much higher densities (Gilliam et al. 1993). Moreover, the size-related nature of the predation risk posed by A. hartii differs from that posed by H. malabaricus and C. frenata both of which are more specialized predators and large enough to capture mature guppies.

Because biological gradients including species richness and predation risk tend to increase from upstream to downstream (Matthews 1998) we might have expected a correlation between richness and predation risk for one or more of our 6 categories. However, in no case did we detect a significant correlation between overall species richness and either biomass or numerical abundance across the 16 sites (see Appendix). This shows that richness cannot necessarily be used as a surrogate for risk, in this system at least, given its spatial scale and the gamma diversity involved (Deacon et al. 2017).

The 3 predation categories used here allow us to gauge the robustness of the gradient to differences in composition. We find some evidence of stability across categories, but note that the rankings of a few sites are dramatically affected. For example, using numerical abundance, the rankings for the broad category are generally consistent with those when the intermediate category is applied (a significant correlation coefficient of 0.96, P = 0.001), despite 6 additional species being included. Using biomass, all comparisons of rankings within categories show significant correlations (P = 0.05).

Although the categories we used in this study were informed by the most up to date information on the ecology of fish in Trinidad’s Northern Range, we recognize that other categorizations of predator community may be equally valid. In all, 21 species of fish occur in these sites (Deacon et al. 2017). We included the ten species recognized as being, at least occasionally, piscivorous. It is likely that guppies are also predated on by some of the omnivorous species in the assemblages (Seghers 1973), but a lack of information on the diet of most species in the system makes justifying which to include difficult. This gap is surprising given how well-studied other aspects of the system are, and further stomach content analysis of all species, supported by stable isotope work, would be an extremely valuable addition to the guppy literature. Expanding our understanding of predator diet and behavior may allow future studies to “weight” the contribution of different predators according to their relative risk from a guppy’s point of view. We also acknowledge that at some sites predation by decapods (Millar et al. 2006), birds (Seghers 1974), and even spiders (Deacon et al. 2015) contributes to the overall strength of predation pressure, and this should also be taken into account by researchers. A further issue, that we have not addressed here, is that the size distribution of predators, within and amongst species, will also affect risk.

Ecologists have long debated whether numerical abundance and biomass provide equivalent insights into the way in which species subdivide resources (Sugihara 1989; Pagel et al. 1991; Taper and Marquet 1996; White et al. 2004). Morlon et al. (2009) showed that the shape of the species abundance distribution often changes with abundance currency. In addition, studies at a local scale have found little evidence of a consistent relationship between numerical abundance and biomass (Stead et al. 2005; Saint-Germain et al. 2007; White et al. 2007; Ehnes et al. 2014). Our results support these findings and further indicate that site rankings, in terms of predation risk, are not conserved across currency. As such, they add more weight to the conclusion that numerical abundance and biomass cannot be assumed to be surrogates for one another in ecological studies. Investigators thus need to make informed decisions about definitions of predation risk in the context of local knowledge.

Within most sites, there is considerable temporal variation in both predator numerical abundance and predator biomass, which is what one would expect given natural fluctuations in population size around an average value (Magurran and Henderson 2010). However, the concordance coefficients indicate that the ranking of sites remained similar over the 6 years of our study, particularly when numerical abundance was the currency. This suggests that even sampling a site once may be sufficient to position this site appropriately on the predation gradient, meaning that short-term studies can still produce meaningful results. It also suggests that other factors, such as climate, may be affecting all sites across the Northern Range similarly, keeping rankings fairly consistent over time.

Biomass produced less temporal consistency in rankings over time than numerical abundance, regardless of the category of predation used (see Figure 5). It is possible that this is because movements of just a few large individuals in or out of a site will have a considerable impact on the total biomass, but only a minor effect on the numerical abundance (Knouft 2002). Again, the goals of a study will be key in deciding which currency to use. Biomass is thought to be a better indicator of resource allocation (Tokeshi 1993; Thibault et al. 2004; McGill et al. 2007; Török et al. 2016) whereas numerical abundance reflects population dynamics (Brook and Bradshaw 2006).

Our data represent annual predator abundance for a period of 6 years. Shorter-term variation in predation risk, such as seasonal differences, or longer-term changes over evolutionary time, are outside of the scope of this study. Nonetheless, both are worthy of further exploration in relation to the Northern Range system (Brown et al. 2006; Carroll et al. 2007; Kotrschal et al. 2017; Barbosa et al. 2018), especially for those interested in exploring the evolutionary reasons behind the guppy’s remarkable flexibility and colonization success (Deacon and Magurran 2016), and it would be interesting to further investigate how these different timescales affect the predation risk gradient.

Conclusions and Recommendations

Biological gradients tend to be underutilized in studies of evolutionary ecology, not least in the famous “natural laboratory” of Trinidad’s Northern Range. Here, we have demonstrated that the classic high-risk/low-risk dichotomy of fish assemblages in this system represents a spatial gradient, underlain by temporal variation.

Of the approaches compared in this study, using numerical abundance reveals a marginally clearer gradient and more consistent ranking of sites over time, whereas rankings using a narrow category of predation appear to be more robust to changes in abundance currency. However, before making specific recommendations for the use of one currency or category over another, it will be important to validate these patterns with actual levels of predation risk in these habitats. This is a realistic prospect in this well-studied system and would be a valuable next step. Further, it is important to keep in mind the specific goal of a study when selecting an approach, as different approaches may reveal different features of the system.

Our study suggests that a consideration of predation gradients will lead to a more nuanced understanding of the role of predation risk in behavioral and evolutionary ecology. For example, in the context of Northern Range guppies, better appreciation of temporal change in community structure within sites, as well as of differences in community composition between sites, may help explain the evolution of individual variation in behavior and decision-making. This additional knowledge could also shed new light on the evolution and maintenance of individual variation in male color patterns. Similar opportunities will exist in other systems. Understanding that predation risk is not just spatially and temporally variable, but also depends on how the predator community is defined in terms of both composition and abundance, will be key here. It also emphasizes that the definition of predation risk used in investigations should be justified and clearly reported in all cases.

Table A1.

Pearson’s correlation coefficients and P-values when species richness for each of the 16 sites is plotted against their predator abundance, using both currencies for each of the 3 categories

Currency	Category	Pearson’s correlation coefficient	P-value
Numerical abundance	Narrow	0.20	0.460
	Intermediate	−0.21	0.442
	Broad	−0.22	0.422
Biomass	Narrow	0.26	0.337
	Intermediate	0.05	0.863
	Broad	0.21	0.434