Small-scale phenotypic differentiation along complex stream gradients in a non-native amphipod

Selective landscapes along complex stream gradients

Environmental gradients provide a unique opportunity to evaluate the relative importance of various forms of both natural [1, 2] and sexual selection [3, 4] for creating intraspecific phenotypic variation [1–7]. Changes of environmental conditions can be accompanied by divergence of phenotypic traits when it allows local populations to reach a new fitness peak within a spatially heterogeneous selective landscape [6, 8]. Studies focusing on environmental gradients acknowledge the multitude of different ecological gradients and described, for instance, gradual variation of abiotic selection factors across latitudinal [9–11], or altitudinal [12] gradients, as well as gradients formed by environmental stressors like temperature [13], salinity [14, 15], or acidification [16, 17]. Other studies focused on biotic selection factors such as predation [18], as well as intra- [19] or inter-specific competition [20]. Beside spatial variation, selection pressures can also vary substantially across time [21, 22].

Streams and the associated limnic ecosystems are shaped by a widespread form of an environmental gradient in which an array of abiotic and biotic selection factors vary systematically in space and time, e.g., from source regions over smaller tributaries to slow-flowing lowland river sections [23–25] or between seasons [21, 23]. Predictable longitudinal changes include, amongst others, changes in primary carbon sources, with allochthonous organic material typically representing the main food source for headwater communities and autotrophic primary production becoming more important at downstream sites [23, 25]. Accumulation of various forms of nutrients usually results in a higher nutrient availability in downstream sections of rivers [23], while in anthropogenically transformed landscapes also pollutants tend to accumulate in downstream sections (e.g., pharmaceuticals [26], pesticides [27], or plastics [28]). Temporal changes in different river sections include, amongst others, flow regimes and water temperatures: source habitats usually show little seasonal fluctuation with respect to abiotic conditions—such as water temperatures [29]—and closely mirror the conditions seen in the groundwater bodies feeding them [23]. The stability of abiotic conditions rapidly decreases with increasing distance from the source(s), especially when the stream is no longer covered by forest canopy, such that water temperatures will attain greater variance due to solar input [25, 30]. These river sections often experience pronounced fluctuations in abiotic conditions, as they undergo recurrent catastrophic flooding (e.g., after snow melt [31]), or temporal desiccation either during hot/dry seasons [32, 33] or as a result of climate change [34]. Abiotic conditions reattain a more stable state in downstream river portions, where multiple tributaries interconnect to form an extensive wetland system [23, 25].

Despite the generality of some of these characteristics of river gradients at a larger geographical scale [23–25], the situation in reality tends to be more complex, not only because numerous landscape features alter the hydrogeological setting at a local scale [35, 36], but also because several forms of anthropogenic stressors [37, 38] contribute to the complexity of selective landscapes along river gradients [23–25]. Obviously, it is close to impossible to empirically capture all abiotic and biotic components characterizing a given stream gradient. Moreover, several factors—including those that were empirically assessed and those not assessed—can be tightly interrelated (e.g., oxygen content and temperature, or the occurrence of certain pollutants [39]). Empirical studies reporting intraspecific trait divergence along stream gradients and ascribing this divergence to one particular component of the gradient are, therefore, prone to overlooking that divergence may have been driven instead by another, interrelated selection factor, or by multiple factors [6, 40].

Using amphipods to study phenotypic differentiation along stream gradients

Acknowledging the complexity of selective regimes in stream gradients [23–25, 37, 38], where multiple selective forces may have complex synergistic, non-additive, or antagonistic effects [80] (potentially overriding certain evolutionary responses, e.g. [81]), we refrain from providing unidirectional predictions for the different trait suites studied herein; an overview of potential predictions can be found in Table 1. We considered physiological traits that are likely to be chiefly under natural selection from abiotic agents of selection (e.g., gill surface areas should vary in dependence of local temperature and oxygen regimes [46]), reproductive life history traits known to be under natural selection from both abiotic and biotic selection factors (brood size/fecundity and offspring size [18, 82, 83]), and traits that are primarily under sexual selection (e.g., size of the male gnathopods, which are involved in mate securing and defense [84]; Table 1).

Intraspecific phenotypic differentiation in invasive amphipods

In our present study we examined a non-native species that reached Central Europe approximately 150 years ago [41] (see Additional file 1: Material S2 for additional information). While the examination of a species that colonized Central Europe relatively recently already drew the focus of our study towards questions related to contemporary evolution [85, 86] or adaptive phenotypic plasticity [87], we included yet another level of questions related to rapid evolutionary change/adaptive plasticity and compared streams with and without anthropogenic thermal pollution (starting approximately 50 years ago [88, 89]). Studying biotic responses to thermal pollution has become increasingly popular to forecast the potential effects of global warming [90]. Thermal pollution should primarily affect naturally selected traits (e.g., via lower oxygen availability), but could (mediated, e.g., by longer reproductive seasons) also affect the strength and direction of sexual selection (Table 1). Given that studies on inter- and intra-specific body size variation have been at the forefront of research on climate-related evolution for centuries (Bergman’s rule, originally formulated for endotherms [91] but partly also applied to ectotherms [47]), we assessed body size as another dependent variable (Table 1).

Our study was motivated by the idea that freshwater amphipods—including invasive species or those that expanded their distribution ranges relatively recently—could be established as a future model for several questions in evolutionary ecology, especially questions related to local adaptation [6, 8], contemporary evolution [85] or adaptive phenotypic plasticity [87], life-history evolution [18, 82], and climate change biology [92]. Our present study is centred on two major questions:

(1) Do we find small-scale phenotypic differentiation in non-native G. roeselii approximately 150 years upon their arrival in Central Europe [41] along repeated river gradients? To answer this question, we sampled amphipods from 16 sites along two streams in Germany (Fig. 1) and assessed an array of phenotypic traits (Table 1).

(2) If phenotypic differentiation is uncovered, can we ascribe at least parts of this variation to known drivers of both natural and sexual selection processes (or triggers of adaptive plasticity), including selection factors arising from anthropogenic activities (Table 1)?

Study species and sample collection

We focused on a ‘naturalised’ [93] non-native species, Gammarus roeselii, which has spread throughout Central Europe and established stable populations since the nineteenth century [41]. Gammarus roeselii is usually considered as a downstream-adapted species [94, 95]. However, several populations in downstream sections of Central European streams were recently replaced by the more recent invaders Dikerogammarus villosus and Echinogammarus ischnus [96–98]. On the other hand, G. roeselii is nowadays found in at least some Central European first order streams at altitudes above 400 m (this study).

We collected specimens in two subcatchments within the Rhine catchment (Fig. 1). The Erft is located in a heavily exploited, urbanized area and in downstream sections receives cooling water from several nearby coal power plants [88, 89], which significantly increases average water temperatures (mean ± SD during winter: 6.1 ± 1.1 °C at thermally non-polluted sites and 12.4 ± 1.2 °C at thermally polluted sites; during summer: 17.5 ± 3.0 °C at thermally non-polluted sites and 23.7 ± 0.3 °C at thermally polluted sites), while our sampling sites in the Kinzig are located within a long-term ecological research area (LTER [99]), and its tributaries do not receive thermally polluted water (Additional file 1: Table S1). We caught amphipods using a multi-habitat kick-sampling method [100]. To account for potential seasonal variation of the traits considered here, we sampled all sites twice, once in February (winter sampling) and once in May of 2017 (summer sampling). All captured specimens were immediately preserved in 96% ethanol until they were processed in the laboratory. As G. roeselii has been reported to occasionally co-occur with the congeners G. fossarum and G. pulex [94, 101], we verified species identity of all specimens under a stereomicroscope (Olympus SZX12 with an Olympus DF PLFL 0.5X PF objective) using standard identification keys [102, 103]. Altogether, we dissected n = 1648 G. roeselii from 16 sampling sites (n = 373 from thermally polluted and n = 1275 from thermally non-polluted sites; Additional file 1: Table S3).

Environmental parameters

Life history and morphometric data

We collected information on male and female life-histories (see Additional file 1: Material S3 for additional information) and morphological traits from 50 to 128 individuals per population (Additional file 1: Table S3). All measurements of distances or areas were conducted under a stereomicroscope (OLYMPUS SZX12) with an OLYMPUS DF PLFL 0.5X PF objective and an attached OLYMPUS SC30 camera connected to a computer. We used the software Cell^A (Olympus) for all linear and area measurements.

Sexing, adult body size and weight

Specimens were sexed according to external sexual characteristics: males were identified by the presence of genital papillae (Fig. 2a) which lie ventrally in the middle of the bases of the 5th pereopods [103]. Female gammaridean amphipods usually have four pairs of oostegites between pereopods 2–6 [105]. In some cases, intersexual individuals [106] with both genital papillae and oostegites were observed (Fig. 2b; see Additional file 1: Material S4 for additional information and a posteriori analyses on data from our present study).

We determined each specimen’s body length [mm] from the tip of the rostrum to the telson tip (Fig. 2c). To do so, alcohol-stored specimens were placed in lateral position in a wax-filled Petri dish and were carefully straightened using preparation needles. We ensured that all specimens were similarly stretched before measuring their body length. We also measured head capsule length as an alternative proxy for body size [107]. However, both traits were strongly collinear (Pearson’s r = 0.92; Additional file 1: Figure S1), and so we used body length throughout. Before we assessed somatic dry weight [mg], we removed all eggs from the female marsupium, after which we dried specimens for 24 h at 60 °C in a drying oven (Heraeus Type B5042). We stored dried individuals in a desiccator containing dry silica gel to prevent water-uptake while cooling down to room temperature. Afterwards, we weighted specimens to the nearest 0.1 mg on a micro-scale (SI-234 Denver Instrument).

Offspring-related phenotypic traits

For each female, we determined fecundity by carefully removing all eggs from the brood pouch (the marsupium [66, 105]; Fig. 2g) with a fine brush and counting them. Data from females harbouring stage 7 embryos (see below) were excluded, as some of the juveniles may already have left the marsupium. To estimate egg volumes [mm³], we measured the longest and shortest axis of each egg in 2D view and used an ellipsoid formula to approximate egg volume (see [67]). We identified embryonic developmental stages according to earlier descriptions for Gammarus spp. [67, 108, 114]. A detailed description can be found in Additional file 1: Material S5, and exemplary photos are shown in Fig. 3.

Statistical analyses

We asked whether and how environmental parameters used to characterize stream gradients and condensed into principal components (Table 2) drive divergence of multiple phenotypic traits. To this end, we applied generalized least squares models (GLS) using the gls-function, implemented in the nlme package in R [115] fitted by maximum-likelihood estimation and assuming a Gaussian error distribution. To account for spatial autocorrelation, we fit our GLS models with an autoregressive (AR-1) correlation structure [116]. Each of the following variables was modelled separately: ‘body length’, ‘dry weight’, ‘length of 1^st antennae’, and ‘length of 2^nd antennae’. For females, we additionally analysed ‘fecundity’ and ‘egg volume’ and for males ‘gill area’ and ‘gnathopod size’. We included the explanatory variables ‘environmental PCs 1 – 4’, ‘sex’ (in all cases where data were available from both sexes) and ‘body length’ to account for ontogenetic changes in trait values. For the model on egg volume, we additionally considered ‘developmental stage’ as an explanatory variable to account for ontogenetic changes in egg/embryo size during development. The statistical significance of explanatory variables in the linear models was evaluated by Type III ANOVAs (using the Anova-function implemented in the car package [117]). Assumptions of normality of residuals was assessed visually by inspecting QQ-plots and plots of residuals vs. fitted values (Additional file 1: Figure S3). In the case of significant effects, we calculated and visualized predicted marginal effects of environmental variables on dependent variables using the ggeffects package [118]. All statistical analyses were performed in R 3.5.1 [119].

Body size and condition-related adult traits

Adult body size decreased significantly from upstream to downstream sites (i.e., along environmental PC 1; χ² = 23.92, p < 0.001; Table 3a; Fig. 4a) and with increasing conspecific densities (PC 2; χ² = 4.0, p = 0.046; Fig. 4b). Likewise, adult body weight decreased significantly along environmental PC 2 (χ² = 9.89, p = 0.002; Table 3b; Fig. 4c) and showed seasonal fluctuation, with greater body weight being observed in winter (PC 4; χ² = 6.33, p = 0.012; Fig. 4d).

Physiological traits

Gill surface areas decreased from upstream to downstream conditions (environmental PC 1; χ² = 15.20, p < 0.001; Table 3c; Fig. 5a). Moreover, gills surface areas were greater during winter than summer (PC 4; χ² = 10.37, p = 0.001; Fig. 5b).

Traits involved in intersexual communication and mate guarding

We found a pronounced sexual dimorphism in the length of the 1^st (χ² = 152.65, p < 0.001) and 2^nd antennae (χ² = 1381.67, p < 0.001; Table 4), with males possessing longer antennae than females. Moreover, antennae length showed gradual variation along all four environmental PCs (see significant effects of PC 1 – PC 4; Table 4). Specifically, both the 1^st (χ² = 63.08, p < 0.001; Fig. 6a) and 2^nd antennae (χ² = 48.98, p < 0.001; Fig. 6b) became longer along environmental PC 1. The significant interaction term of ‘sex × environmental PC 2’ reflects that the length of the 1^st antennae of males increased to a greater extent along PC 2 (i.e., with increasing population densities and increasingly female-biased sex ratios) than those of females (χ² = 5.66, p = 0.017; Fig. 6c). Both antennae increased along environmental PC 3 (1^st antennae: χ² = 34.23, p < 0.001; Fig. 6d; 2^nd antennae: χ² = 11.39, p = 0.001; Fig. 6e) and were longer during winter compared to our summer sampling (PC 4; 1^st antennae: χ² = 39.96, p < 0.001; Fig. 6f; 2^nd antennae: χ² = 13.62, p = 0.001; Fig. 6g). Notably, we did not find significant divergence in male gnathopod length along environmental gradients (Table 4c).

Reproductive characteristics

Female fecundity (numbers of eggs per brood)—size-corrected by inclusion of the covariate ‘female body length’ (Table 5a)—decreased significantly from up- to downstream sites (i.e., along environmental PC 1; χ² = 20.08, p < 0.001; Fig. 7a) and with higher conspecific densities (PC 2; χ² = 78.20, p < 0.001; Fig. 7b). Egg size did not show gradual variation along any of the components of the stream gradient (Table 5b).

Ecological significance of phenotypic differentiation

What are the potential ecological implications of our findings? Recent years have seen a steep increase in numbers of publications that either analyse empirical data (e.g., [120]) and/or use mathematical modelling (e.g., [121]) to understand ecological phenomena over large geographic scales (macroecology [122]). These approaches, however, usually treat species as “evolutionarily inert” entities with fixed (species-specific) ecological traits affecting ecosystem processes like nutrient and energy flow [123–125]. Nevertheless, both theoretical considerations and empirical evidence suggest that intraspecific phenotypic variation can affect ecosystem processes [126–128]. For example, divergence in foraging traits of predatory fish species can alter the structure and dynamics of ecological communities by decreasing the mean body size, total biomass, and species richness of crustacean zooplankton [129]. The consequences of phenotypic divergence for ecological processes should be particularly strong when the species under consideration acts as a keystone species or ecosystem engineer [130] or when it is simply very abundant [127]. Amphipods are considered to be key components of limnic food webs, as their shredding activity accelerates leaf fragmentation, produces faecal pellets, and transfers nutrients into secondary production, all of which are vital for maintaining diverse aquatic food webs [131–135]. Hence, rapid evolutionary divergence of amphipods along stream gradients is likely to result in multiple feedbacks between evolutionary and ecological processes [136–138]. For example, G. roeselii shows the strongest degree of body size-dependency in leaf litter decomposition rates amongst several species of native and non-native amphipods in the Rhine drainage (see Additional file 1: Figure S2 in ref. [139]). Hence, the strong decline in body size across the stream gradient reported here will have an effect on decomposition rates and food web dynamics. Ecological studies have merely just begun to acknowledge the role of rapid phenotypic divergence in the provisioning of ecosystem services [128, 136, 140], and our present study identifies invasive amphipods as prime candidates for future studies in this direction. Specifically, future studies could quantify decomposition rates [139] of phenotypically divergent populations to explicitly integrate differences in functional responses as a consequence of the observed trait variation and to understand how they translate into differences in ecosystem functioning [140–142].

Drivers of phenotypic divergence along the examined stream gradients

The native distribution range of G. roeselii is characterized by its karst topography, with spatio-temporally diverse ecological conditions along stream gradients [143, 144]. Therefore, G. roeselii might have a high potential to inhabit variable environments, either by phenotypic plasticity (see below) or through adaptive evolutionary divergence (local adaptation; see also Additional file 1: Material S1). In support of this idea, a substantial amount of the observed phenotypic trait divergence we found in our dataset could be explained by selection factors that we condensed into principal components. We will base our discussion largely on a priori predictions of how single components of the river gradient should affect trait divergence (Table 1), but acknowledge that additional factors that were not assessed here (e.g., predation pressure [145, 146], competition with congeneric species [147]) may have driven parts of the observed divergence.

Outlook and future research

Phylogeographic studies on amphipods are becoming common practice to study how past geological events affected macroevolutionary processes [166–168]. By contrast, phenotypic divergence within and between taxonomic groups—as reported here on an intraspecific level for G. roeselii—received little attention so far. Uncertainty remains regarding the underlying mechanisms that caused the observed phenotypic variation, and future studies are warranted that address the potential heritability of the traits studied herein and the role of adaptive phenotypic plasticity [87, 169]. Additionally, parts of the observed phenotypic divergence could also arise from non-adaptive plasticity in response to stressful environments [87] (Additional file 1: Material S1). Future studies could provide broad-sense heritability estimates by rearing individuals under common-garden conditions [170], or make an attempt to rear individual broods from known parents to provide narrow-sense heritability estimates (h_n², via mid-parent – offspring regressions [171]). To experimentally address some of the newly generated hypotheses, future studies could use the sister genus Hyalella, which is easy to breed under laboratory conditions. Hyalella spp. are routinely used in North America as a model organism for ecotoxicological assays [172] and contemporary evolution to anthropogenic sources of selection (pesticide exposure) has been described [173]. On the other hand, the relative role of plastic responses could be addressed experimentally by assessing the variance of traits, for example, in reciprocal transplant experiments [174, 175]. Finally, genome- and transcriptome-screens could provide insights into the underlying mechanisms of phenotypic divergence, and the widespread occurrence of G. roeselii offers several independent populations to test for parallel adaptive divergence. Altogether then, our present study identifies amphipods are promising models to study phenotypic diversification along ecological gradients made up by both natural and human-induced selection factors. Furthermore, their key-role in freshwater ecosystems renders amphipods an excellent system in which to investigate feedbacks between evolutionary and ecological processes.