Effects of changes in isotopic baselines on the evaluation of food web structure using isotopic functional indices

Site description and sampling protocol

Five waterfalls (with vertical drops ranging from 18 to 72 m), from small streams to large rivers (widths ranging from 6 to 50 m), were studied in Southern Quebec, Canada (between 46−47°N and 72−73°W). Their catchment areas are situated on the Canadian Shield (corresponding to a metamorphic geological formation), making the running water weakly conductive and slightly acidic (ranging from 20 to 50 µS cm⁻¹with an average pH value around 6.3 ± 0.4 at investigated sites). To use changes in DIC-δ¹³C values brought by waterfall CO₂-outgassing as an integrative indicator of changes in algal isotopic baselines, each site was sampled at two locations immediately upstream and downstream of the waterfall (hereafter above- and below-waterfall), and the maximum distance between the two sampling points was 300 m. Paired sampling locations were also selected to have similar environmental conditions (water velocity, riverbed substrates, water depth, surrounding vegetation cover, canopy cover, etc.). Therefore, as habitat structures in each waterfall paired-sites were similar, food web structures of above- and below-waterfall invertebrate communities were also expected to be similar.

C-gas sampling and carbon stable isotope analysis of DIC

Water chemistry (partial pressure of CO₂: p CO₂, DIC-δ¹³C) was measured at each sampling site to characterize biogeochemical effects of waterfalls and quantify expected shifts in algal δ¹³C values. The selected sites were visited between 1 to 2 times in spring and summer (in early May and late June 2016). p CO₂ was measured using the headspace method (Campeau & Del Giorgio, 2014). 30 mL of water sample was collected from approximately 10 cm below the water surface, using a 60 mL polypropylene syringe, and 30 mL of ambient air was added into the syringe to create a 1:1 ratio (ambient air: water sample). Then, the syringe was vigorously shaken for 1.5 min to equilibrate the gases in the water and air fractions. 30 mL of the headspace was then injected into a 40 mL glass vial, prefilled with saturated NaCl solution (360 g L⁻¹ at 20 °C), through a butyl rubber septum. A second needle was used to evacuate the excess of NaCl solution. Vials were kept inverted for storage, and headspaces were analysed using a Shimadzu GC-8A Gas Chromatograph with flame ionization detector at University of Quebec at Montreal (Montreal, Canada). pCO₂ in water samples was then retrocalculated using the headspace ratio, water temperature and ambient air concentrations of CO₂ at studied sites. pCO₂ measurements were performed in duplicates for each sampling site. Supersaturation ratios (SR) were also calculated by dividing the gas water partial pressure by the atmospheric CO₂ concentration.

At each sampling site, 500 mL of water sample was also collected in duplicates in early May and late June 2016 to analyse DIC- δ¹³C. Water samples were filtered at 0.2 µm using nitrocellulose membrane filters and stored for a maximum of 72 h in the dark at 4 °C until analysis. 150 µL of phosphoric acid (H₃PO₄; 85%) was added into 12.5 mL amber borosilicate vials to ensure that all DIC content in the water sample would be converted into CO_2. Then, vials were flushed using Helium during 10 min to ensure a full evacuation of ambient air. 4 mL of water sample was injected in He-flushed vials through the rubber septa using fine needles, and vials were equilibrated at 20 °C for 18 h. DIC- δ¹³C was obtained using with a ThermoFinnigan Gas Bench II coupled to an Isotope Ratio Mass Spectrometer (IRMS), and results were expressed as the delta notation with Vienna Pee Dee Belemnite as the standard: δ¹³C (‰) = ([R_sample/R _standard] – 1) ×1000; where R = ¹³C/¹²C. Sample measurement replications from internal standards (C1 = −3.0‰, and C5 = −22.0‰) produced analytical errors (1σ) of ± 0.3‰ (n = 17).

Invertebrate sampling and carbon stable isotope

Each sampling station was sampled in early July 2016, and benthic invertebrates were collected in riffle sections using a kick-net (0.1 m², 600 µm mesh size). Equal sampling effort was applied to each habitat type within above- and below-waterfall sites. Invertebrate specimens were sorted immediately in the field into taxonomic groups and transported in the dark at 4 °C back to the laboratory 4–8 h later to be frozen at −20 °C until analysis. A small isotopic deviation can be observed using this method (see also Feuchtmayr & Grey, 2003; Wolf et al., 2016), but we assumed that this effect was the same for all samples. Invertebrates were identified at the genus level (Merritt & Cummins, 1996), and specimens were then classified into different feeding groups as herbivores, detritivores, and predators (Thorp & Covich, 2009; Merritt & Cummins, 1996; Electronic supplementary material S2). Samples were then dried at 60 °C for 72 h and ground into fine powder. Carbon (δ¹³C) and nitrogen (δ¹⁵N) stable isotopic ratios were then analysed using an Isotope Ratio Mass Spectrometer interfaced with an Elemental Analyser (EA-IRMS) at University of Quebec at Trois-Rivieres (Trois-Rivieres, Canada). Results were expressed according to the delta notation (see above). Sample measurement replications from three internal standards (STD67: δ¹³C = −37 ‰ and δ¹⁵N = 8.6‰, UTG40: δ¹³C = −26.2‰ and δ¹⁵N = −4.5‰, and BOB1 δ¹³C = −27‰ and δ¹⁵N = 11.6‰) produced analytical errors (1 σ) of ± 0.02‰ for δ¹³C values and ± 0.2‰ for δ¹⁵N values (n = 148).

Isotopic functional indices and data analysis

Non-metric Multidimensional Scaling (NMDS) was used to visualize dissimilarities among/within invertebrate communities at waterfall sites, and the Bray–Curtis index was used to measure dissimilarities of invertebrate communities based on presence/absence data. T-tests were also performed on trophic guilds δ¹⁵N values to compare their trophic positions in food webs at above- and below-waterfall sites.

Means of δ¹³C and δ¹⁵N values of all individuals for each species calculated at each sampling site were used to derive six IFIs following Layman et al. (2007): δ¹³C range (CR), δ¹⁵N range (NR), total area of the convex hull encompassing all the observations (TA), mean distance to centroid (CD), mean nearest neighbour distance (MNND) and standard deviation of nearest neighbour distance (SDNND). Layman’s IFIs can be grouped into isotopic functional richness (CR, NR and TA); isotopic functional divergence (CD); and isotopic functional evenness (NND and SDNND). As isotopic functional richness indices (CR, NR and TA) provide a quantitative indication of the extent of the isotopic niche space of the entire community and are calculated using the dispersion of species in the δ-space (Layman et al., 2007; Jackson et al., 2011), we hypothesized that those IFIs should be more sensitive than those based on the distribution of species in the δ-space (CD, MNND and SDNND; Appendix. 1). All indices were calculated using SIAR package for R (Parnell & Jackson, 2013). Principal component analysis (PCA) was also performed to display changes in structural properties of above- and below-waterfall invertebrate communities and provide an overview of relationships between IFIs and changes in DIC- δ¹³C values. All statistical analyses and plots were performed using the R 3.5.2 software (R Core Team, 2018).

Waterfalls, community structure and basal resources

NMDS (Fig. 2) and taxonomic list (Table 1) showed only little changes in taxonomic composition and suggested that compositional differences among communities were higher among waterfalls than within each paired site. Furthermore, no changes in δ¹⁵N values of aquatic consumers (Fig. 4) were reported suggesting that the analyzed organisms occupied similar trophic positions in the food webs (Cabana & Rasmussen, 1996)) and might therefore feed on similar diets above and below the waterfall sites. Therefore, in our study, waterfalls did not significantly impact the food web structure of invertebrate communities, and these results could strengthen previous findings showing the absence of major effect of waterfall on invertebrate community composition in four tropical rivers (Baker et al., 2016).

Waterfall CO₂-outgassing and associated shift in DIC- δ¹³C values should induce punctual changes in algal δ¹³C values and could also help to decipher the respective contribution of allochthonous and autochthonous carbon to lotic food webs. With few exceptions, δ¹³C values of trophic guilds were higher at below-waterfall sites than those of above-waterfall samples (Fig. 3), supporting the view that aquatic invertebrates mainly feed on in-stream algae (Tanentzap et al., 2017). However, differences among trophic guilds observed in our study might also suggest varying reliance on algae (Fig. 5). Surprisingly, large isotopic shifts were also observed for detritivores (but were in general smaller than those for herbivores), suggesting an important dependence on autochthonous sources for these organisms theoretically relying on detritus (Fig. 5; McNeely, Clinton & Erbe, 2006). Therefore, our study could support the hypothesis of the prevalence of autochthony in river and stream food webs (Brett et al., 2017; Tanentzap et al., 2017) but could also emphasize the issue of trait plasticity for inveterate leading to differences between actual and theoretical feeding habits.

Sensitivity of IFI to changes in isotopic baselines

IFIs have become increasingly used in aquatic ecology (Olsson et al., 2009; Abrantes, Barnett & Bouillon, 2014; Dézerald et al., 2018; Burdon, McIntosh & Harding, 2020), but IFI concept mainly relies on untested assumptions. In this study, we consider waterfall systems as a natural experimental set-up to quantify impacts of changes in isotopic baselines on ecological inferences from IFIs. We consider that waterfall CO₂-outgassing and associated shift in DIC- δ¹³C values should induce punctual changes in algal δ¹³C values and therefore help to understand how changes in isotopic baselines impact upon IFIs.

Notable differences in in the expected direction in IFI values were reported between above- and below-waterfall sites (Fig. 6), and changes were likely driven by shifts in DIC- δ¹³C values (Fig. 7B). Isotopic functional evenness and divergence indices (CD, SDNND and MNND) were only slightly impacted by changes in isotopic baselines (Fig. 6). Indeed, those indices were calculated based on the distribution of species in the δ-space (Layman et al., 2007), and ecological inferences were therefore not strongly impacted by changes in the vertical or horizontal distribution of specimens in the isotopic space. In contrast, isotopic functional richness indices (CR, NR and TA), calculated on species dispersion in the δ-space and providing a quantitative indication of the extent of isotopic niche space of the entire community, were strongly influenced by changes in algal δ¹³C values (Figs. 6 and 7). Moreover, differences in CR, NR and TA values between above- and below-waterfall sites often exceeded the range of changes previously reported in literature (Rigolet et al., 2015). These results could strengthen previous findings that these indices are very sensitive to changes in ranges of consumer δ¹³C and δ¹⁵N values (Brind’Amour & Dubois, 2013; Syväranta et al., 2013). Hence, isotopic functional richness indices (CR, NR and TA) might be frequently misinterpreted in river studies comparing food webs across sites and/or over time with fluctuating isotopic baselines.

Our study suggested that the reliability of IFI inferences can be strengthened by identifying changes in IFIs driven by variability in isotopic baselines. In that vein, different propositions have already been made in the literature to improve the understanding of food web structure and resource partitioning in consumers (e.g., Jabot et al., 2017). The most promising approach likely consists of increasing the number of isotopes studied (like hydrogen and sulphur: Doucett et al., 2007; Proulx & Hare, 2014) and including other types of data (such as gut content, fatty acids contents or compound specific isotope analysis). Indeed, the combination of these complementary proxies will provide new insights on actual energy pathways through food webs. By enabling a better understanding of trophic interactions in food webs, future IFI-based studies will contribute to better document food web structural properties.