Uncovering the Grinnellian niche space of the cryptic species complex Gammarus roeselii

Sampling

We re-sampled sites studied by Grabowski et al. (2017), in order to investigate the stability of local populations of the G. roeselii species complex. In addition, to get an impression about small-scale distribution of each genetic lineage, we further sampled nearby rivers. Furthermore, we also considered additional sample sites in the more recently colonized central Germany (Jourdan et al., 2019; Weigand et al., 2020). Exact coordinates of the 42 sampling sites and their corresponding reference to Grabowski et al. (2017) is given in Table S1. Sampling took place during a field campaign in September 2021. Individuals were sampled by a multi-habitat kick-sampling using hand nets (Bioform V2A; mesh size 500 µm). We preserved all captured individuals in 96% ethanol and stored them at ~10 °C until they were processed in the laboratory for molecular analyses. All sampling was done in close cooperation with local partners; the species are not subject to any protection status. Permission was granted by the Hellenic Ministry of Environment and Energy (Permission No: YПEN/ΔΔΔ/7316/280) and the Albanian National Agency of Protected Areas (Permission No: 194, issued on 16.02.2021). Maps were created using ArcGIS Desktop (v10.8; Esri, Redlands, CA, USA).

Molecular methods

We performed DNA extraction on two individuals of each sampling site using the standard protocol for human or animal tissue and cultured cells of NuceloSpin® Tissue Kit (Machery-Nagel GmbH, Düren, Germany). From each individual two to three pereiopods of one side were carefully removed and placed in a 1.5 ml reaction tube. A total of 25 µL Proteinase K and 180 µL Buffer T1 was added to the sample and they were incubated for 4 h at 56 °C. The lysis, binding and washing of the DNA was carried out as stated in the protocol for human or animal tissue and cultured cells of NuceloSpin® Tissue Kit (Machery-Nagel GmbH, Düren, Germany) for the “High yield, high concentration”-approach.

To identify the different MOTUs cytochrome c oxidase subunit I (COI) gene fragments were amplified using polymerase chain reaction (PCR) with three different primer pairs depending on the MOTU evaluated by Grabowski et al. (2017). Primers used were LCO1490 and HCO 2198 (Folmer et al., 1994), UCOIF and UCOIR (Costa et al., 2009) and COIGrF and COIGrR2 (Grabowski et al., 2017). The full list of primers per individual used is presented in the Table S2. We performed PCR in 30 µL volume, containing 0.4 µL BSA, 3.6 µL MgCl2 (25 mM), 3.0 µL 10 × Taq-Buffer (BioLabs, Boston, MA, USA), 2.4 µL dNTPs (40 mM), 0.4 µL Taq Polymerase (Biolabs, Boston, MA, USA), each primer with a volume of 1 µL (10 mM), 14.2 µL of ddH2O and 4 µL DNA template. Amplification of a 650-bp fragment of the COI gene was carried out using the (PCR) conditions stated in Mamos et al. (2014): Initial denaturation at 94 °C (3 min), 35 cycles of denaturation at 94 °C (20 s), annealing at 50 °C (45 s) and elongation at 65 °C (1 min) and a final extension at 65 °C (2 min). PCR products were purified using the protocol of the PCR clean-up kit of NuceloSpin® (Machery & Nagel GmbH, Düren, Germany). Successful PCR amplification was verified using an aliquot of 5 µL in a GelRed® (VWR) 1.0% agarose gel. The one-sided forward sequencing of the PCR products was performed at Eurofins Genomics’ GATC service “LightRun Tube” (Eurofins Genomics Germany GmbH, Ebersberg, Germany). All sequences were checked using the BLASTn search in GenBank (Altschul et al., 1990). Site 14 did not contain G. roeselii and is therefore not considered in the results. The remaining sequences were then assembled with sequences published in Grabowski et al. (2017) for each MOTU and trimmed to a length of 522 bp in MEGA X (Kumar et al., 2018; Ver. 10.2.6). Alignment was performed with ClustalW (Larkin et al., 2007). With a neighbour-joining tree the 82 sequences were assigned to the MOTUs according to the MOTUs published in Grabowski et al. (2017) (see Fig. S1 and Table S3). We calculated the K2P (Kimura two-parameter) distance as a measure of genetic distance (Kimura, 1980) between and within the MOTUs using MEGA X (Kumar et al., 2018; Tables S4 and S5, Fig. S2). The newly acquired COI sequences are available online on the Barcode of Life Data System (BOLD) via their assigned IDs, provided in Table S2.

Environmental and climatic parameters

Freshwater ecological conditions along the latitudinal gradient from the mountainous Continental Germany to Mediterranean Greece differ significantly (Tockner et al., 2022). In an attempt to capture the variation in abiotic conditions, we measured on site flow velocity (P670; Dostmann electronic, Wertheim, Germany), pH (PHC201; Hach HQ40d multi, Loveland, CO, USA), oxygen content and saturation (LDO101; Hach HQ40d multi, Loveland, CO, USA), as well as conductivity (CDC401; Hach HQ40d multi, Loveland, CO, USA). We extracted the altitude of the sampling sites from Google Earth (http://earth.google.com/). For evaluating the quality of the sampling site sediment and water samples were taken. The water samples were analysed for their concentrations of nitrate (NO₃⁻), nitrite (NO₂⁻) and phosphate (PO₄³⁻) with a photometer (FinwellPro for ponds, Glen Iris, Victoria, Australia; MDE GmbH & Co. KG, Puchheim, Germany) and its according measuring kits. Carbonate hardness (HCO₃⁻) was measured using a colorimetric kit (Merck MColortests) (for the environmental data, see Table S6).

We further gained information on hydro-environmental watershed and river characteristics from Domisch, Amatulli & Jetz (2015). This dataset contains spatially continuous and freshwater-specific set of environmental variables for a standardized 1 km river network grid. Domisch, Amatulli & Jetz (2015) used upstream accumulation techniques to obtain watershed contributions for each pixel, providing information on land cover, climatic, river topographic and geological conditions of the years 1970 to 2000. For our modelling approach, we considered flow length (size of the catchment upstream of the sampling site as the sum of 1 km contributing grid cells), cultivated and managed vegetation (lc_avg_07; “cultivated land cover”), urban/built-up land cover (lc_avg_09; “urban land cover”), annual mean temperature (Bioclim1; “annu. mean temperature”) and annual precipitation (Bioclim12; “annu. precipitation”). These variables are computed for a 1 km grid cell of a terrestrial dataset, where the mean of each variables is weighted over the distance of the upstream catchment, starting from the corresponding sampling site (Domisch, Amatulli & Jetz, 2015). This allows us to account for environmental factors in the upstream watershed that have a significant influence on physicochemical conditions of the respective sampling site (for the climatic data, see Table S7 in the Supplementary).

Effect-based assessment of chemical contamination at sampling sites

It was shown that when taking human influences into account, niche differences within a single taxon can be predicted more accurately as seen for example in parakeet distributions (Strubbe et al., 2015). To determine the extent of contamination at each sampling site, we performed a set of in vitro assays to evaluate the baseline toxicity and endocrine activity on the sediment samples. Sediment acts as a sink for various environmental pollutants such as hydrophobic organic contaminants, accumulating and storing pollutants over time (Keiter et al., 2006). As a result, sediment can provide a record of historical contamination levels in a particular location. Specifically, endocrine disruption caused by environmental estrogens can be responsible for disrupting the reproduction and development of organisms, thus a screening for endocrine activity can evaluate the hormonal burden within the aquatic environment (Bergman et al., 2013). As benthic dwellers amphipods are exposed to sediment extensively and sediment analysis is a good indicator of actual exposure to contaminants (Nguyen, Muyssen & Janssen, 2012). In order to consider this pollution information in our analysis, we took sediment samples at each sampling site and stored them at 4 °C until processing. In the laboratory, we then used yeast reporter gene assays to screen for endocrine activity (YES) at the human estrogen receptor α (hERα; Routledge & Sumpter, 1996) and dioxin-like activity (YDS) at the aryl-hydrocarbon receptor agonists (AhR; Stalter et al., 2011) according to Giebner et al. (2018). For the analysis, we eluted 10 g of fine sediment (< 20 µm) in 50 mL methanol an placed the sample on an orbital shaker for 1 h at 220 rpm (GFL 3017; GFL Gesellschaft für Labortechnik GmbH, Burgwedel, Germany). Afterwards, we used an ultrasonic bath (Sonorex RK 52 H; Bandelin electronic, Berlin, Germany) for 10 min to completely detach remaining substances from the sediment. After centrifugation at 4400 rpm for 5 min the supernatant was decanted using a 1.5 µm glassfilter to withhold all suspended soils. The methanol was then removed with a rotary evaporator for 5 min at 56 °C (Heidolph Laborota 4000-efficient, vacubrand CVC 2000; Heidolph Instruments GmbH & Co. KG, Schwabach, Germany; VWR RC-10 Digital Chiller, VWR International GmbH, Darmstadt, Germany) and the remains were dissolved in 500 µL DMSO. Initially, we also tested for androgenic activities at the androgen receptor hAR (Sohoni & Sumpter, 1998) using a YAS assay, but could not find any significant loads in the sediment, so we did not include this in the further analysis. A microtox assay (labelled as “base. tox”) with the bacterium Aliivibrio fischeri was carried out to test for non-specific baseline toxicity of the sediment. A. fischeri is a Gram-negative marine bacterium that produces light through bioluminescence. This bacterium has a high sensitivity to toxic chemicals, making it a useful indicator of acute toxicity. The EC₅₀ value in the microtox assay is the concentration of the sediment that causes a 50% reduction in bioluminescence of the bacteria (i.e., the lower the EC₅₀ value, the more toxic the sample is considered to be). In the further calculations we subtracted the EC₅₀ value by 100 to have high values labelled as more toxic than low values. The microtox assay was in accordance to the ISO guideline 11348-3 (International Organization for Standardization, 1998) and adapted after Escher et al. (2008) for a 96 well plate format. All analysed data is available in the Tables S6 and S7.

PCA and niche-modelling

To compare the Grinnellian niche space of G. roeselii sensu lato and that of the individual MOTUs, we used all the environmental data we collected in combination with online available climatic data to evaluate niches for each MOTU. To do this, we first conducted a principal component analysis (PCA) using the function ‘dudi.pca’ within the R package “ecospat” (Di Cola et al., 2017) using the environmental, chemical and climatic factors (see “Environmental and climatic parameters” and “Effect-based assessment of chemical contamination at sampling sites”). All data was standardised prior performing the PCA so the mean value for all variables is at zero and has a unit variance value by setting “scale = TRUE” within the ‘dudi.pca’ function (Di Cola et al., 2017). Based on a PCA with all the environmental variables, we visualised the niches using the R package “ecospat” (Di Cola et al., 2017). The niche visualisation was based on the first two principal components (PCs) which explained 39% of the variance of the environmental variables. The G. roeselii niche space was then plotted underneath the environmental PCA as a representation of the overall environmental conditions that are suitable for G. roeselii. Overall, the G. roeselii niche space was more or less evenly distributed along all environmental axes. The first two PCs mainly represented a gradient in conductivity, carbonate hardness (HCO₃⁻), annual mean temperature, precipitation and flow length (PC1, explained = 24.1%) as well as a gradient in pH, altitude and mean flow velocity (PC2, explained = 14.9%; see Table S8).

To quantify the niche overlap between each MOTU, the probabilities were calculated where only the first MOTU (i.e., ‘unfilling’), both MOTUs (i.e., ‘stability’) or only the second MOTU (i.e., ‘expansion’) occurs. All MOTUs were compared with each other and displayed in a cross-table. We calculated Schoener’s D (as a measure for niche overlap; Broennimann et al., 2012) and Warren’s I (as a degree of niche similarity of two niches; Warren, Glor & Turelli, 2008) for each MOTU combination. Schoener’s D and Warren’s I range from 0 to 1, with 0 indicating complete dissimilarity and 1 indicating complete similarity. In addition, the p-value of these statistics were calculated using the equivalency test with the hypothesis of a niche overlap that is less equivalent/similar to a randomly generated niche (overlap.alternative = “lower”) and niche conservatism is assumed (expansion.alternative = “lower”, stability.alternative = “higher” and unfilling.alternative = “lower”). We repeated the same analysis with PCs 3 and 4 (see Figs. S3 and S4) that explain additional 11.6% and 9.7% respectively. PC3 and PC4 mainly describe the effect of agriculture, with cultivated land use and O₂ saturation (PC 3; explained variance = 11.6%) as well as Nitrate (NO₃⁻) and Nitrite (NO₂⁻) (PC 4; explained variance = 9.7%; see Table S8).

Distribution of the MOTUs sampled

Our COI barcoding approach confirmed the presence of (at least) five MOTUs in our study region (Fig. 2). Most of the sequences acquired at the sampling sites selected accordingly to Grabowski et al. (2017) have the same MOTU classification except for two sampling sites which are our sampling sites 19 and 11 (i.e., ID 13 and 11 respectively in Grabowski et al., 2017). Probably due to our sample size of n = 2, we could not detect co-occurring MOTU E or J here. In Slovenia and Germany we only found MOTU C (sites 1–9). MOTU G (sites 10–21) was found in many stream systems in Albania, as well as in lakes Ohrid and Prespa (Albania and Greece). We also found MOTU GG (site 13 and 15), but we have refrained from including it in the niche analysis, as this MOTU is nestled within MOTU G so we treated this as a MOTU G. In Greece, we found a clear north-south structuring, with MOTU A (sites 22–29) inhabiting the northern part and MOTUs K (sites 35–40) and L (sites 30–33 and 42) the central part of Greece. The occurrence of the MOTUs is therefore in accordance to Grabowski et al. (2017) with the exceptions for sampling site 19 and 11 where only MOTU G was found.

Characterizing the G. roeselii sensu lato niche space

We found G. roeselii sensu lato mainly in larger rivers, but also occasionally in first-order streams. Accordingly, environmental conditions were highly variable, where e.g., altitude varied from sea level up to 905 m, conductivity from 210 to 1100 µS/cm and pH from 7.53 to 9.76 (for a complete list of environmental conditions, see Tables S7 and S8). Pollution load of the sediment also varied strongly, from near-pristine rivers in Slovenia with low nitrate, nitrite and phosphate, to rivers in Greece, characterized by high levels of nitrate and nitrite (100 and 1.64 mg/L respectively; see Table S6). The niche space depicted by the first two PCs was mostly influenced by the conductivity, flow length and annual precipitation along PC1 and the altitude and pH along PC2 (Fig. 3). Pollution variables “baseline toxicity”, “dioxin-like activity” and “estrogenic activity” correlated positively along PC1 and PC2 with dioxin-like activity having a relatively strong impact (Table S8). Nitrate and nitrite concentrations correlated with the occurrence of cultivated land.

Niche overlap of the MOTU niches

We used the multidimensional PCA information to transform the direct niche space into niches for each MOTU (Fig. 4). The MOTU-specific niches always cover only part of the G. roeselii sensu lato niche. MOTUs C, G and A each have a broad niche, with G and A showing the highest niche overlap and similarity (D and I; Fig. 4). MOTU C and G, both with broad niches but distributed along different PCs, have some significant overlap (D = 0.20; p = 0.0495) and similarities (I = 0.26; p = 0.0495). MOTU L has the smallest niche of all sampled MOTUs, showing no overlap or similarity with niches of MOTU C, G and A (D and I; p < 0.01). The data acquired for MOTU L represents a niche that is only present within the niche of MOTU K (unfilling = 0.00). Additionally compared to MOTU K, there is only a small overlap (D = 0.16; n.s.) and the similarity is modest (I = 0.38; n.s.). Therefore the niche of MOTU L and K can’t be separated by statistical means. MOTU L and A on the other hand show significantly no overlap or similarity (D = 0.00; I = 0.02; both p < 0.01) and the niches of MOTU L and A are therefore significantly different from each other. The southernmost MOTU K and northern most MOTU C show a low overlap but the niches have some significant similarities (I = 0.26; p = 0.0396).

When considering further environmental variables (i.e., PC 3 and PC 4), we found a less clear differentiation. Here, only MOTU C differed significantly from MOTU A (I = 0.18; Fig. S4). The overlap and similarity (I and D respectively) are high for all MOTUs (mean I = 0.47; mean D = 0.35), indicating that all MOTUs were also found in agricultural dominated regions.

Niche differentiation with genetic distance

The linear regression of the K2P distance and niche overlap (Schoener’s D) did not show a significant correlation (Pearson’s correlation coefficient r = −0.2677; p = 0.4546; Tables S3 and S4, Fig. S2). Thus, we could not find a correlation of increasing niche differentiation and greater genetic distance between MOTUs.

Stable distribution pattern

Our samplings confirmed the molecular patterns of the G. roeselii complex sampled by Grabowski et al. (2017) 15 years ago. Moreover, our intensified sampling around the previous sites showed that on a small-scale (e.g., within a catchment) usually only one exclusive MOTU occurs. This is in concordance to most of the findings for Gammarus species in Europe (e.g., G. balcanicus: Mamos et al., 2016, 2021; G. fossarum: Wattier et al., 2020; G. roeselii: Grabowski et al., 2017). Nevertheless, our sample size was small and cannot fully exclude the possibility that co-occurence with rare MOTUs has been overlooked. Co-occurrence of MOTUs is possible as it was shown for other gammarids such as G. fossarum (Copilaş-Ciocianu & Petrusek, 2015; Bystřický et al., 2022) or G. pulex (Grabner et al., 2015). Assuming that MOTUs of a morpho-species use similar resources (e.g., they compete for the same resources, such as food and/or space), competitive exclusion would be expected (Scudo, Ziegler & Lotka, 1978; Tilman, 1982; Webb et al., 2002). If a shift in resource use is unsuccessful, competitive exclusion usually results in the local extinction of one of the species/MOTUs, as the more efficient or better-adapted competitor outcompetes and excludes the other from the area (Tilman, 1982) and a spatial separation of the competing species is evident (e.g., Niphargus hrabei and N. valachicus; Copilaș-Ciocianu et al., 2017; corophiid Amphipods; Borza, Arbačiauskas & Zettler, 2021). Therefore, finding more than one MOTU at the same site seems to be more of an exception than a rule.

MOTU specific niche differentiation

Our niche modelling approach revealed that some MOTUs within the G. roeselii species complex occupied different niches, partially confirming our second hypothesis. MOTU C, G and A all have a broad niche, suggesting that they have a high tolerance to environmental conditions and anthropogenic pollution (Gaston & Spicer, 2001). Interestingly, the wide range of MOTU C is not directly reflected in a larger niche, compared to MOTUs such as A and G. Nevertheless, MOTU C showed a clear differentiation from the other MOTUs (Schoener’s D ≤ 0.2). With its extended geographic distribution, MOTU C is a special case in the G. roeselii species complex. MOTU C is to be assumed the sister clade of MOTU A that separated around 4.5 million years ago (Grabowski et al., 2017). The relatively quick expansion of MOTU C in Central West Europe may have been promoted by human activity in combination with a warming climate (Csapo et al., 2020). The ability to cope with different abiotic conditions is represented by the broad Grinnellian niche, supporting the assumption that MOTU C is an ecological generalist (Csapo et al., 2020). Moreover, the differentiation from other MOTUs suggests that a niche shift has taken place, in favour of colder and larger river systems. This is not surprising, considering the colonization of the vast Danube river system (Sommerwerk et al., 2022) and that of northern river systems such as the Rhine catchment (Wantzen et al., 2022). The possibility to cope with a wide range of environmental conditions due to a broad generalistic niche may be one of the reasons for its widespread distribution in anthropogenically influenced rivers in Central Europe today (Csapo et al., 2020). We also found cases of large niche overlap for example for MOTUs G and A which occur in similar geographic regions. Although these two MOTUs have phylogenetically and geographically separated over 8 Million years ago (Grabowski et al., 2017), they still share similar Grinnellian niches, as commonly found in other allopatrically separated species (Wiens, 2004; Wiens & Graham, 2005). Kozak & Wiens (2006) discussed the role of niche conservatism and allopatric speciation, demonstrated in sister taxa of salamanders Desmognathus and Plethodon. Here, niche conservatism was found even though they only show small climatic differences between the closely related species (Kozak & Wiens, 2006). The strong overlap and highly similar niches of MOTUs A and G in addition to the clearly separated river systems suggest that they have conserved their niche requirements across geographic regions, indicating a degree of niche conservatism. Most likely their ancestral population was lacustrine, but became separated due to geological changes such as the forming of river channels and hydrological networks (Grabowski et al., 2017).

We found the smallest and most clearly separated niche for MOTU L. MOTU L and K belong to the southern members of the G. roeselii species complex, that phylogeographically split from the northern clades about 18 Million years ago (Grabowski et al., 2017). It is assumed that in Greece the southern MOTUs are spatially scattered and only occur in very localised and mostly isolated waterbodies (Grabowski et al., 2017). This is already an indication for possibly smaller niches, simply because amphipods can hardly spread across watersheds and are therefore severely limited in their dispersal capacity. This means that MOTU L colonizes the whole Pinios river system that is shaped by intensive use of water for agriculture, with 95% of available water resource used for crops and farmland (Skoulikidis, 2016). Such river systems are prone to enormous low flows or complete desiccation during summer, making it a system of high human selection pressure (Cooper et al., 2013; Skoulikidis, 2016; Karaouzas et al., 2018). It seems that MOTU L can tolerate this kind of anthropogenic pressure, but further research is needed to determine the specific adaptations or mechanisms that allow it to persist in these challenging conditions (see “Higher tolerance of generalist species in global change?”). In general, when interpreting the niches, it should be considered that the sample size was small. This is mainly due to the large geographical range, the complex analysis of environmental samples and the need for genetic characterization of all organisms. Nevertheless, we were able to sample a representative range of environmental conditions for all MOTUs especially when considering the limited and narrow distribution of most of the MOTUs (Grabowski et al., 2017). Only for the widely distributed MOTU C, which we were only able to collect selectively, we can assume an even greater tolerance to environmental conditions than already shown here.

Overall, we could not find evidence for a relationship between phylogenetic distance and niche differentiation, as shown with the non-significant correlation between K2P distance and Schoener’s D. Thus, our third hypothesis is not supported; it appears more as if the niche differentiation is case-specific, and in some cases the Grinnellian niches of the MOTUs differ in size and space. This was also the case for the amphipod genus Niphargus where different cryptic species showed ecological niche differences but the degree of differentiation was not linked to the genetic distances of the cryptic species (Fišer et al., 2015). The linkage between niche differentiation and phylogenetic relatedness could be more relevant on the clade level as pointed out by Wiens (2008) and Fišer et al. (2015).

In general, it remains to be considered that mitochondrial COI markers reflect geographical diversity well, but in some cases they tend to overestimate species levels in cryptic species and are incongruent with species delimitation by using nuclear markers (Mamos et al., 2021; Hupalo et al., 2022). For future ecological niche modelling of cryptic species it might therefore be useful to further back up the information from COI markers by the additionally use of nuclear markers.

Higher tolerance of generalist species in global change?

The Grinnellian niches represent the abiotic conditions and pollution loads under which we collected the individuals. Actual tolerances may differ. Our initial characterisation of the Grinnellian niche now provides the opportunity for follow-up questions. For example, it remains to be tested whether the small niche of MOTU L may be associated with lower tolerance to the abiotic factors considered here. Vice versa, it can be assumed that MOTU C, G and A tend to have larger tolerances. Such empirical tests of the environmental thresholds found here, would also show whether the species have certain Grinnellian niches only because of geographical barriers, i.e., it can be checked whether the approach used here incorrectly predicts too small Grinnellian niches for dispersal-limited/spatially structured species (e.g., MOTU L).

In recent years, Gammarus spp. have increasingly been used as test organisms for studying adaptation processes to various anthropogenic stressors (Schulz & Liess, 1999; Bundschuh & Schulz, 2011; Shahid et al., 2018a, 2018b; Siddique et al., 2020; Siddique, Shahid & Liess, 2021; Grethlein et al., 2022). So far, mostly intraspecific differences between populations in anthropogenically stressed regions were compared with pristine regions. The G. roeselii species complex provides an excellent framework to study adaptations in five closely related cryptic species that occur at both, the inter- and intraspecific levels.

Niche differences as an indicator of trait divergence?

The pronounced difference in the Grinnellian niches, especially of MOTU L and C, calls for investigation under an integrative framework to include a morphological/phenotypic characterization and focus on trait differentiation. Follow-up studies, should quantify a set of characteristics in MOTU L compared to MOTU A, C and G in order to see the role of ecological influence on the phenotype. If phenotypic differentiation is visible within the cryptic species complex, they should be more pronounced in MOTUs where we have observed niche differentiation (Whitlock, 1996; Fišer, Robinson & Malard, 2018). One might specifically expect MOTU L to have a different spectrum of phenotypic differentiation compared to MOTUs with a broad niche size (Sexton et al., 2017). Whether this phenotypic differentiation might then be more or less pronounced and in correlation with niche size differences remains to be seen. Additionally, if phenotypic differentiation is found, it should be discussed whether these differences also affect functional ecosystem processes (Bickford et al., 2007; Fišer, Robinson & Malard, 2018).