The study area comprised ten tributaries (Strahler streams order of 4 to 6) of the North Saskatchewan River basin in Alberta, Western Canada (Fig. 1). Three reaches (WMD2, BMD2, and BMD3) are located in the highly developed and populated Edmonton Metro Area, Alberta, while the other seven reaches (COW1, BAP2, ROS2, POP1, BEA1, SMO1, and VER4) are located in less developed and lower populated rural and/or natural areas. Of the latter group, COW1, BAP2, ROS2, and POP1 are tributaries that originate in the montane and foothills regions of the basin and feature colder mean summer water temperatures. In contrast, BEA1, SMO1, and VER4 are warmer water streams that originate in the prairies (Fig. 1). It is important to note that, to investigate potential competition for and/or partitioning of dietary resources between F. virilis and native species in this study, we assume that dietary resources are limited in all 10 reaches.

Figure 1. 

Locations of ten study reaches in the North Saskatchewan River basin. Reaches where F. virilis are absent are represented by empty white circles. Reaches occupied by F. virilis are represented by circles filled with a black solid circle. Unique reach identification codes are located near each reach’s location marker.

We sampled F. virilis during the summer (June-August) of 2020 in each tributary along a 200–300 m reach using modified minnow traps with 5.7 cm diameter openings (Mangan et al. 2009; De Palma-Dow et al. 2020). Individual traps were tied to nylon rope 3 m apart in sets of six to form a trapline (Rosewarne et al. 2014). Individual traps in each trapline were baited using salmon-based cat food-filled perforated film canisters (one canister Purina Friskies® per trap) (Mangan et al. 2009). Each trapline was then affixed to the riverbank by a loop and rebar stake on the upstream end. A total of four traplines were deployed at each reach for a total of 24 traps per reach. Because crayfish are nocturnal and most active at night (Styrishave et al. 2007), baited traplines were left overnight to increase the chance of capture. The morning following deployment, traplines were retrieved from the water. All captured crayfish were enumerated and humanely euthanized using a 15-minute ice bath and pithing in accordance with animal handling and ethics regulations (Animal Use Protocol No.: AUP00003578). Whole specimens were frozen at -20 °C prior to sample processing and analysis. F. virilis relative abundance was estimated using catch per unit effort (CPUE) and calculated for each reach as the mean number of individuals captured in a single 24-hour overnight baited trapline survey and was reported in units of individuals per trapline (Zale et al. 2012) (Suppl. material 1: table SS1).

Fish and benthic invertebrates were sampled within the same week as crayfish were sampled in each tributary to minimize effects of any changes in flow conditions over time. Each 300 m reach was subdivided into six 50 m transects within which fish were sampled via backpack electrofishing in a sweeping systematic pattern. All fish captured were identified to species, fork length measured, and released back into the river from which they were caught, except for individuals of target fish species. Target native fish species included invertivore/herbivore secondary consumers; Longnose Dace (Rhinichthys cataractae [Valenciennes, 1842]), Lake Chub (Couesius plumbeus [Agassiz, 1850]) and Trout Perch (Percopsis omiscomaycus [Walbaum, 1792]); native detritivorous consumers Longnose Sucker (Catostomus Catostomus [Forster, 1773]) and White Sucker (Catostomus commersonii [Lacepede, 1803)]); and one native piscivorous species: Burbot (Lota lota [Linnaeus, 1758]). Target species were selected due to their ubiquity in the basin so that direct species comparisons could be made across all sampled sites. Individuals of target fish species were humanely euthanized via single-blow blunt force trauma followed by pithing to ensure death (Research License 20-3812 RL). Whole specimens were frozen at -20 °C prior to sample processing and analysis. All fish sampling and euthanasia was conducted in accordance with the animal handling and ethics regulations (AUP No.: AUP00003578), and under a valid Research License issued by the Government of Alberta (RL# 20-3812).

Benthic macroinvertebrates were collected for isotopic baseline calculation. At all ten reaches, benthic macroinvertebrate communities were sampled using a triangular kick net (400 µm mesh) in a single zig-zag pattern, sweeping over erosional zones (riffles and runs) a minimum of 3 m of where traplines were set and standardized to three-minute sampling effort. Collected material was placed into enamel pans where invertebrates were separated from other material using forceps and wash water bottles. All samples were preserved in 70% ethanol (Hobson et al. 1997). Benthic invertebrate samples were identified and enumerated to the family or genus when further identification was possible and necessary to resolve functional trait discrepancies. All identification of aquatic insect taxa and their functional feeding traits were determined following the work of Merritt et al. (2019).

During dissections, a sample of dorsal muscle tissue was collected from each fish specimen. Also recorded were specimen wet weight (g), fork length (mm), total length (mm) and sex (m/f/juv) based on presence of female or male gonads. Additionally, the stomach contents of L. lota (a larger predator species) were inspected for evidence of F. virilis consumption. A sample of abdominal (tail) muscle tissue was collected from each of the 52 crayfish specimens. Weight (g), carapace length (mm), total length (mm), and sex (m/f) were also recorded.

Prior to stable isotope analysis, all fish and crayfish tissue samples were preserved at -20 °C, while all benthic macroinvertebrate samples were preserved whole, in 70% ethanol. Fish dorsal muscle and crayfish abdominal muscle tissue samples were freeze dried at -55 °C and 0.015 Bar for 24 hours in a LABCONO® FreeZone 1 Liter Benchtop Freeze Dry System (Labconco 2021) to constant weight. Benthic invertebrate samples (separated by taxonomic family) were dried whole to constant weight at 60 °C for 24 hours in a Precision® Compact Gravity Convection Oven (Thermo Scientific 2009). Once dried, each sample was ground into a homogenous powder and weighed into a 6X8 mm tin capsule to 0.4000–0.4999 mg using the UMX2 Ultra-microbalance (Mettler Toledo 2004). All samples were analyzed for δ¹⁵N and δ¹³C ratios at the Natural Resources Analytical Lab in Alberta, Canada using ThermoScientific DeltaV Advantage isotope ratio mass spectrometer (IRMS) coupled to the ThermoScientific FlashSmart Organic Elemental Analyzer and ConfloIV. All values are reported in delta notation (‰) relative to international standards: Pee Dee Belemnite (vPDB) for δ¹³C and atmospheric N₂, Vienna Air (VAIR) for δ¹⁵N. Calibration was performed using seven-point normalization which includes six certified reference material and two international reference standard material LSVEC and IAEA-N2 with known values of δ¹³C = -46.6‰ and δ¹⁵N = 20.3 ‰, respectively. In-house pea grain standard with δ¹³C = 8.94‰ and δ¹⁵N = 3.234‰ were analyzed every 12 samples to monitor drift and ensure accuracy. Accuracy and precision are in relative standard deviation and were better than ± 0.20‰ δ¹³C and δ¹⁵N. Accuracy was within error for secondary reference material. Blanks and empty capsules were included at the beginning of each run cycle to monitor for contamination.

Prior to analysis, raw isotope data were inspected for carbonate contamination using multiple linear regression analysis to compare δ¹³C and percent carbon of samples of each reach (Jardine et al. 2003). Next, we investigated carbon depletion caused by lipid richness. Samples with C:N ratios > 4 were considered lipid rich and in need of correction (Post et al. 2007). We detected no lipid richness. Further, because ethanol fixation is known to have a measurable effect on stable isotope carbon values (e.g. Ventura and Jeppesen 2009; Hogsden and McHugh 2017; Blechinger et al. 2024) we corrected for the effect of preservation in 70% ethanol on our benthic invertebrate δ¹³C signatures by applying the mass balance approach as described by Ventura and Jeppesen, (2009) for sample sets where the original C:N is unknown. Specifically, we calculated each corrected δ¹³C value using the equation:

δ¹³C_EtOHCorr = δ¹³C_EtOHfixed – [D ((C:N_EtOHmean – C:N_EtOHfixed) / C:N_EtOHmean)]

Where δ¹³C_EtOHfixed is the original ethanol contaminated δ¹³C signature of the given sample, D is equal to 4.65 – the difference coefficient for ethanol preserved freshwater benthic invertebrate samples that was determined by Ventura and Jeppesen (2009), C:N_EtOHfixed is the C:N ratio of the given sample, C:N_EtOHmean is the mean C:N value of all samples, and δ¹³C_EtOHCorr is the final δ¹³C signature of the given sample after ethanol preservation correction. No ethanol preservation correction was applied to δ¹⁵N values as the balance of literature indicates this effect is negligible (e.g. Hobson et al. 1997; Hogsden and McHugh 2017; Hajisafarali et al. 2023). Finally, C.-commersonii and C.-catostomus samples were combined into a single group called Catostomus spp. because the sample size of C.-catostomus (n = 6 over all sampled reaches) was insufficient for subsequent stable isotope analysis. This combination was based on the two species sharing a close taxonomic lineage (same genus) and similar life history, morphology, and diet (Scott and Crossman 1973).

To account for spatial isotopic variability among reaches and to compare stable isotope metrics between reaches, primary consumer benthic macroinvertebrate samples from each reach were used to calculate isotopic baselines. We used primary consumer benthic invertebrates for baseline source estimates rather than long-lived primary consumers such as clams and snails because these long-lived organisms were unavailable. Benthic invertebrates have been demonstrated as reasonable indicators of baseline δ¹³C and δ¹⁵N values in previous studies (e.g. Anderson and Cabana 2007; Jackson and Britton 2014). Baseline δ¹⁵N values were calculated following the suggestions of Anderson and Cabana (2007) to convert raw δ¹⁵N and δ¹³C values into trophic position (TP) and corrected δ¹³C (δ¹³C_Corr), respectively. The family Elmidae was the most nitrogen depleted benthic invertebrate family and present in 60% of sampled reaches. Where Elmidae did not occur, the next collector family with the lowest δ¹⁵N value was selected and corrected to Elmidae (Anderson and Cabana 2007). We used δ¹³C and δ¹⁵N values from the most nitrogen depleted benthic invertebrate families within each reach to calculate the mean δ¹³C and δ¹⁵N values, pooled across all reaches. The pooled δ¹³C and δ¹⁵N baseline values were then used to calculate TP and δ¹³C_Corr.

Trophic position was calculated for each consumer’s muscle tissue sample using the pooled baseline δ¹⁵N value and individual raw δ¹⁵N values by substituting them into the single source trophic position model described by Post (2002):

TP_con = λ + (δ¹⁵N_con - δ¹⁵N_base)/Δ_n

Where TP_con is the trophic position of the consumer, λ is the trophic position of baseline organisms (λ = 2 for herbivorous benthic invertebrates), δ¹⁵N_con is the isotopic nitrogen value of the consumer, δ¹⁵N_base is the calculated pooled baseline δ¹⁵N value, and Δ_n is the trophic enrichment factor (TEF) equal to 3.4 ± 1‰) which is the applicable fractionation value first determined by Post (2002). 3.4‰ is used as the fixed nitrogen TEF for studies conducted on wild populations (Dionne et al. 2016), where trophic TEFs specific to species and diet (McCutchan et al. 2003) are not available.

We corrected δ¹³C consumer values based on the pooled δ¹³C benthic invertebrate baseline value as was used for calculation of trophic position using the equation described by (Olsson et al. 2009).

δ¹³C_corr = (δ¹³C_con – δ¹³Cµ_baseline)/CR_baseline

Where δ¹³C_corr is the basal isotopic corrected δ¹³C value of the consumer, δ¹³C_con is the raw δ¹³C value of the consumer, δ¹³Cµ_baseline is the pooled δ¹³Cvalue, and CR_baseline is the range of source δ¹³C values across all reaches. Trophic position and corrected δ¹³C values were used in all subsequent statistical analyses. Summary statistics of trophic position and δ¹³Ccorr are provided by species (Table 1) and by reach (Suppl. material 1: table SS1).

To determine the size and position of each species’ basin-wide and realized trophic niches, maximum likelihood fitted small sample size corrected standard ellipse area (SEA_c) was calculated and used as a measure of each species core trophic niche width using the R SIBER package (Jackson and Parnell 2020 Package “SIBER” ). SEA_c represents approximately 40% of the spread of the data and is ideal for calculating the core trophic niche of a species when working with small sample sizes (n < 30) (Jackson et al. 2011). To determine the basin-wide trophic niche of each species, we pooled samples by species across all reaches before plotting and calculating the niche width of each species. By pooling individuals of the same species over all study reaches, differential resource availability and interspecific interactions such as competition and predation of the multiple reaches are evened (Baltensperger et al. 2015). Therefore, the standard ellipse area of conspecific individuals that are pooled over all study reaches is representative of that species’ basin-wide niche width.

In contrast, when the species’ standard ellipse areas of individual reaches are plotted separately, they represent the realized niche of the species present in that reach. Here, we assume that all individuals within a reach are subjected to similar interspecific interactions and resource availability. Therefore, we plotted the standard ellipse areas within each reach to inspect the realized niche widths of each species. Niche widths were reported in units of ‰² area.

To detect if dietary resources were being consumed by both native fish and F. virilis, the basin-wide and realized niche widths were inspected for presence and degree of overlap between F. virilis and native fishes. The degree of overlap was calculated as a proportion using the R SIBER package. Proportional overlap values were then calculated as a proportion of the non-overlapping area of the two ellipses using the following equation:

P_overlap = [V_overlap / (V_ellipse2 + V_ellipse1 – V_overlap)]/100

Where P_overlap is the unit-less proportion overlap of the two trophic niches being compared; V_overlap is the ‰² area value of overlap of the two species’ trophic niches being compared; and V_ellipse2 and V_ellipse1 are the calculated trophic niche area of species 1 and species 2, respectively. The final proportional overlap was reported as a percentage between 0% and 100% with an overlap of 0% indicating completely unique ellipses and an overlap of 100% indicating complete overlap.

The following analyses were conducted on three native secondary consumer fish species (NSCFS) (R. cataractae, C. plumbeus, and Catostomus spp.) only. Percopsis omiscomaycus and L. lota were excluded from the following analyses due to insufficient total sample size (n_{P. omiscomaycus} = 3, n_{L. lota} = 11).

To detect potential trophic impacts of F. virilis sympatry on NSCFS in each reach, we calculated the Bayesian estimate of realized standard ellipse area (SEA_B) using the R SIBER Package. We used a null prior distribution to estimate the SEA_B niche widths with 95% probability intervals for each NSCFS population over 20,000 iterative runs of the Bayesian bivariate distribution model. We then compared the SEA_B niche widths of populations that were sympatric with F. virilis and those that were not. SEA_B niche widths were considered significantly different from each other when the 95% probability intervals around the means being compared did not overlap (Jackson et al. 2011; Pettitt-Wade et al. 2015).

In addition to differences in SEA_B, the δ¹³C_Corr carbon range of NSCFS were calculated and compared between F. virilis sympatric populations and allopatric populations to determine if the richness of consumed dietary resources was reduced (narrowed carbon range) when sympatric with F. virilis. Carbon range was calculated as the difference between the greatest individual δ¹³C_Corr value and smallest δ¹³C_Corr value and was expressed in units of Δ‰ (Layman et al. 2007).

To determine if F. virilis sympatry and/or trophic niche overlap may be related to reduced body condition of native secondary consumer fish, the relative weight (Wr) fish condition metric was calculated as described by Wege and Anderson (1978). Wr was calculated using the R FSA Package (Ogle 2019, Package “FSAdata.”), using equations derived by Bister et al. (2000) and Giannetto et al. (2011, 2012). The standard weight intercept and slope values for riffle daces (genus: Rhinichthys) and brook chub were used to calculate the relative weights of R. cataractae and C. plumbeus, respectively. The relative weight intercept and slope values of brook chub were used because the relative weight equation for C. plumbeus does not yet exist and brook chub and lake chub share similar life history, taxonomy, and morphology (Giannetto et al. 2012). The relative weight intercept and slope values of C.-commersonii were used to calculate relative weights of Catostomus spp.. Relative weight is expressed as a percentage of the previously determined standard weight for an individual of that species and fork length (Ogle 2018). For example, a relative weight of 100 indicates that the individual is the exact expected weight for a typical individual of its size and species from a reference population. Relative weights of < 100 or > 100 indicate underweight or overweight individual for its size and species, respectively (Wege and Anderson 1978; Ogle 2018). Mean relative weight and standard deviation of each NSCFS in F. virilis occupied and unoccupied reaches were calculated and a two-sided t-test was used to determine if there were significant differences between the mean relative weight of each species dependent on the presence of F. virilis. The two-sided t-test is sufficiently accurate to determine significant difference in means with sample sizes n ≥ 5 (de Winter 2013).

Percopsis omiscomaycus possessed the smallest basin-wide niche (SEA_c) width of all fish species with an area of 0.051‰² (Table 2). The largest basin-wide niche width of all fish species belonged to R. cataractae with an area of 0.610‰² (Table 2). The basin-wide niche width of F. virilis was found to be 0.067‰², having the second smallest basin-wide niche of all sampled species (Table 2). Faxonius-virilis’ basin-wide niche occupied a moderate trophic position and showed a δ¹³C_corr component nearly triple that of the TP component, which was reflective of the species omnivory (Fig. 2). The basin-wide niche of F. virilis overlapped with three out of five native fish species: L. lota, C. plumbeus, and R. cataractae (Fig. 2). Two of these were moderate overlaps with the basin-wide niche of L. lota overlapping 31.2% with that of F. virilis and the basin-wide niche of C. plumbeus overlapping 23.8% with that of F. virilis (Fig. 2, Table 2). The basin-wide niche of R. cataractae was nearly independent of that of F. virilis with an overlap value of 0.14% (Fig. 2, Table 2).

Figure 2. 

Corrected isotopic carbon (δ¹³Ccorr) and trophic position (TP) biplots showing each species’ core basin-wide isotopic niche width. Isotopic niche widths are expressed in ‰² and were calculated using small sample size corrected standard ellipse area (SEAc) which contains 1 SD around the mean or approximately 40% of the data for each species. Isotopic niches are labeled with the corresponding species’ shorthand name, ellipse color, and marker type. Black open circles = F. virilis; red crossed circles = L. lota; orange open diamonds = C. plumbeus; blue open triangles = R. cataractae; green open squares = Catostomus spp.; and gold hourglasses = P. omiscomaycus. Plotting of core isotopic niches were done using the SIBER R package.

Core realized niches (SEA_C) were plotted for each species in the five reaches where crayfish were present. Realized niches were mostly segregated in isotopic space (Fig. 3, Table 2). Out of a total of seven potential overlap events with native fish species, the realized niches of F. virilis overlapped with those of native fish only once (Fig. 3E). This overlap occurred in reach WMD2 with a minor overlap of 2.57% with L. lota (Table 2).

Figure 3. 

Corrected isotopic carbon (δ¹³Ccorr) and trophic position (TP) biplots showing the core realized isotopic niche width of each species within each reach where northern crayfish were found to be present. Panel letters indicate the specific reach as follows: (A) BEA1, (B), BMD2, (C) BMD3, (D) VER4, and (E) WMD2. Isotopic niche widths are expressed in ‰² and were calculated using small sample size corrected standard ellipse area (SEAc) which contains 1 SD around the mean or approximately 40% of the data for each species. Isotopic niches are labeled with the corresponding species’ shorthand name, ellipse color, and marker type. Black open circles = F. virilis; red crossed circles = L. lota; orange open diamonds = C. plumbeus; blue open triangles = R. cataractae; and green open squares = Catostomus spp.. Plotting of core isotopic niches was done using the SIBER R package.

The mean Bayesian estimated core realized niche width area (SEA_B) of C. plumbeus and R. cataractae in F. virilis sympatric reaches were statistically similar (95% probability intervals overlapping) to those of their conspecifics in F. virilis allopatric reaches (Fig. 4A, B). However, the realized niche of Catostomus spp that were in sympatry with F. virilis were significantly smaller in two reaches (reaches BEA1 & VER4) compared with conspecifics found in two F. virilis absent reaches (POP1 & SMO1) (Fig. 4C). Additionally, all three mean SEA_B estimates of Catostomus spp were lowered in F. virilis present reaches compared to those where F. virilis was absent.

Figure 4. 

Density plots of realized isotopic niche widths (SEA_B ‰²) of the three secondary consumer fish species ((A) C. plumbeus [n = 35], (B) R. cataractae [n = 65], and (C) Catostomus spp. [n = 33]) compared where F. virilis are present vs. absent. Black dots represent the bootstrapped mean SEA_B areas. Blue crosses represent the small sample size corrected standard ellipse area (SEAc). Boxes around means indicate the 95%, 75%, and 50% probability intervals of the SEA_B area. Lower case letters indicate significant differences between mean SEA_B values where different letters indicate significant differences with 95% confidence and like letters indicate statistically similar mean SEA_B values. Unique reach codes appear below their respective bar.

For all three NSCFS, carbon ranges were neither consistently broader nor narrower in F. virilis occupied and absent reaches, indicating no detectable effect of F. virilis on the richness of diet consumed by NSCFS (Fig. 5).

Figure 5. 

Range plot comparing reach specific δ¹³Ccorr ranges (Δ‰) which are reflective of dietary source richness of the three secondary consumer fish species (R. cataractae, C. plumbeus, and Catostomus spp.) compared between where F. virilis are present (grey bars) vs. absent (white bars). Unique reach codes appear on the y-axis for each respective bar.

Finally, the mean relative weights of C. plumbeus, R. cataractae, and Catostomus spp. were 79.56 ± 5.11%, 79.23 ± 9.66%, and 80.59 ± 7.35%, respectively. There was also no significant difference in mean relative weight between NSCFS sympatric and allopatric populations of C. plumbeus (t-test p-value = 0.1772), R. cataractae (t-test p-value = 0.8038), or Catostomus spp. (t-test p-value = 0.7582) (Fig. 6).

Figure 6. 

Violin plots comparing mean body condition (as described by the relative weight condition metric [Wr] and reported in %) of R. cataractae (n = 43), Catostomus spp. (n = 9), and C. plumbeus (n = 13) over all reaches in which crayfish are present (grey) against all reaches where crayfish are absent (white). Significant difference between means is represented by an asterisk (*).

Overlap of basin-wide niches between F. virilis and three native fishes suggested that F. virilis have the potential to consume the same dietary resources as native fish species. Overlap at the basin scale seen between F. virilis and C. plumbeus and R. cataractae is consistent with diet studies showing that these fish species, like crayfish, are known benthic feeders who readily consume benthic macroinvertebrates, macrophytes, and/or benthic detritus (Scott and Crossman 1973; Brazo et al. 1978). Additionally, the overlap between L. lota and F. virilis is reflective of the carnivorous diets of these species. Lota lota are known carnivorous fish which feed on benthic invertebrates such as mayfly nymphs and crayfish when small (51- 305 mm) (Scott and Crossman 1973). The overlap of the basin-wide niche of F. virilis with the lower TP portion of that of L. lota suggests that F. virilis are consuming similar benthic macroinvertebrates as small L. lota in the system – further reflecting F. virils’ omnivory.

In contrast with the basin-wide niche analyses and our original hypothesis, within individual reaches, the realized niche of F. virilis was largely segregated from those of native fishes. Lack of overlap of realized niches suggested that while crayfish do consume the same resources as native fish at the basin scale, F. virilis and native fishes are not consuming a significant amount of the same dietary resources when in sympatry. This lack of overlap suggests that F. virilis is not competing for dietary resources with our study’s native fish species. A possible explanation for this unexpected lack of realized trophic overlap is dietary plasticity. Specifically, F. virilis may be using plasticity in food resource selection to avoid competition for dietary resources with sympatric native fishes, such as C. plumbeus and R. cataractae. Omnivory and the ability of crayfish to exercise dietary plasticity are well documented in the literature (e.g. Momot 1995; Veselý et al. 2020). While F. virilis show a preference for consuming animal tissues (Momot 1995), they are able to successfully subsist on less preferred dietary resources by readily consuming macrophytes (Momot 1995), macroinvertebrates (Hanson et al. 1990), and detritus. The ability of Faxonius genus crayfish, specifically Faxonius rusticus and Faxonius limosus, to be highly plastic in their diet has been seen in response to environmental stimuli such as seasonal changes in dietary resource availability and over-invasions – when a formerly successful invader is supplanted by a new one (Tran and Manning 2019; Linzmaier et al. 2020). During times of high animal tissue abundance crayfish preferentially consumed animal tissue whereas crayfish increased their consumption of diatoms and detritus significantly during times of low animal tissue abundance (Tran and Manning 2019). Veselý et al. found similar patterns of dietary plasticity through time and among populations with different dietary resources (Veselý et al. 2020). The use of dietary plasticity to avoid trophic overlaps with native species is also shown in studies of crayfish invasion. In mesocosm experiments, dietary plasticity facilitated niche differentiation between native and non-native crayfishes (Jackson et al. 2014). Further, invasive crayfish seem to use dietary plasticity to occupy a different trophic niche than sympatric invasive fish species (Jackson and Britton 2014).

Our trophic metric and body condition results are consistent with the segregation of F. virilis’ realized trophic niche from those of native fishes. In all but two cases, the niche widths (SEA_B) and body condition of sympatric NSCFS populations were statistically similar to conspecific allopatric populations. The two cases where niche widths were significantly reduced in the presence of F. virilis were seen in Catostomus spp. which had no overlap with either the basin-wide or realized niches of F. virilis. Therefore, we conclude that the observed reduction in niche widths is likely not associated with F. virilis sympatry. Additionally, the overall similarity of realized SEA_B niche widths, carbon ranges, and body condition of F. virilis sympatric and allopatric native fish populations indicates that F. virilis have had no detectable detrimental trophic effects on sympatric NSCFS. This finding is consistent with the results of our first objective that indicated F. virilis and native fishes are consuming different dietary resources and not participating in resource competition when in sympatry.

The ability of F. virilis to use dietary plasticity to occupy a trophic niche that is unoccupied by native species could facilitate the species’ establishment in currently unoccupied areas of the basin. Our study indicates that F. virilis are currently not competing with common native fishes R. cataractae, C. plumbeus, or Catostomus spp.. However, these species are generalist and generally robust (Scott and Crossman 1973). Currently unoccupied tributaries with different species assemblages containing rare, specialist, and/or sensitive fish species may be vulnerable to F. virilis in different ways. If introduced, F. virilis could exert negative effects on sensitive species by way of indirect and/or direct competition, predation, habitat modification, etc. Further, of all Canadian provinces, excluding the Maritimes and Territories, Alberta has the lowest native fish richness (Scott and Crossman 1973). This is due to dispersion barriers for routes from glacial refugia after the last (Late Wisconsian) glaciation (Nelson and Paetz 1992). As biotic resistance to invasion is positively correlated with biodiversity (Elton 2020), species-poor North Saskatchewan River tributaries may be especially vulnerable to new invasions by F. virilis. However, the relationship between biodiversity and biotic resistance is highly nuanced depending on the invasive species, native community, and abiotic factors and cannot be assumed as positive in all cases (Levine and D’Antonio 1999; Lockwood et al. 2013). To prevent further movement and potential impacts of F. virilis in native fish, watershed managers should continue to implement and practice measures preventing further expansion of F. virilis within the North Saskatchewan River basin and other Alberta watersheds.

Our study has two potential limitations. First, in order to evaluate potential competition between F. virilis and native fishes, we assume that dietary resources are limited in all sampled reaches. However, we did not explicitly quantify resource availability. As limited resources are a requirement for competition to occur, it is possible that, rather than using dietary plasticity, F. virilis may be co-occurring without competition due to ample resources in the sampled reaches. However, we contend that if the latter was the case, we would expect there to be more instances and greater percentages of overlap between realized niches of F. virilis and native fishes. As our results stand, we feel that the most plausible explanation is that F. virilis is using dietary plasticity to exploit a different trophic niche than native fishes. Second, we used ethanol to preserve benthic invertebrate samples prior to stable isotope analysis and baseline calculation. While this was done to preserve samples collected in the more remote western sites where freezing was not a feasible preservation method and applied the mass balance approach to correct for the effect of fixation on carbon signatures, we recognize that this could still have introduced some error into our results. However, comparison of our results with and without the mass balance correction indicate that the patterns gleaned where consistent and likely are robust to ethanol preservation-induced error.

Our study evaluated the trophic effects of F. virilis on three common and generally robust native species. While we found little evidence that F. virilis is competing with these native fishes for dietary resources, our study does not exclude the possibility that F. virilis may be competing for resources with other North Saskatchewan River basin fish species. Further investigation should be made into the trophic effects of F. virilis on native rare and sensitive fish species, as they could be more vulnerable to F. virilis presence than the species we studied. For example, crayfish have been shown to compete with benthic carnivorous fish species for spatial resources (Reynolds 2011). Crayfish have been shown to force juvenile L. lota to leave preferred shelter habitats which can make juveniles more vulnerable to predators (Hirsch and Fischer 2008), while sculpins are displaced from shelters and spend increased time fleeing in the presence of crayfish which resulted in reduced growth rates and lowered body condition (Light 2005).

Our results also do not rule out the possibility that F. virilis exert direct negative effects on native fish by way of predation. For example, instream experiments have shown that crayfish actively prey upon adult benthic darter species (Thomas and Taylor 2013) as well as the eggs and fry of threatened fish species (Fitzsimons et al. 2002). Considering predation in the other direction, our stomach content analysis revealed that F. virilis are being preyed upon and consumed by at least one piscivorous fish species in the North Saskatchewan River basin: L. lota. Lota-lota has been documented to prey upon crayfish as a natural prey item in their native range (Jacobs et al. 2010).

The impacts of F. virilis in the North Saskatchewan River basin may not be limited to fishes. F. virilis have been known to change the species assemblages of benthic macroinvertebrates drastically and decimate native snail and clam biomass (e.g. Hanson et al. 1990; Rodríguez et al. 2005). Future studies would do well to investigate F. virilis’ behavioral interactions with juvenile piscivorous fish; if North Saskatchewan River basin sculpin species are being preyed upon by F. virilis; the mercury concentrations of F. virilis and their contribution to bioaccumulation in benthic predatory fish such as L. lota; and/or the potential effects of F. virilis on benthic invertebrate communities using stable isotope mixing models and diversity indices.

Lastly, while our study investigates the potential competition for trophic resources between NSCFS and invasive F. virilis, we acknowledge that physical and physiological resources such as habitat and preferred temperature can be limiting factors for crayfishes and could impact the interpretation of our results. All five of our sites where F. virilis is present possess a similar, moderate degree of habitat complexity and are centrally located within the North Saskatchewan River basin and therefore have similar summer water temperature regimes (Van Mierlo et al. 2022). Therefore, we contend that the sites that we evaluated are likely similar enough in habitat and preferred temperature to have a negligible influence on the patterns seen in this study. However, investigation of how limited physical or physiological resources interact with and/or affect the potential competition for trophic resources between NSCFS and F. virilis is an interesting area of study that should be pursued in the future.

In conclusion, overlap of F. virilis’ basin-wide niche with those of native fishes indicated that F. virilis have the potential to consume the same resources as and/or compete with native fishes. However, segregation of realized niches showed a lack of resource competition within communities of the North Saskatchewan River basin. Our results suggest that rather than participate in resource competition, F. virilis may be using dietary plasticity to exploit a slightly different trophic niche than those occupied by native fishes and in doing so, avoid competition for dietary resources through resource partitioning. While F. virilis were not found to negatively affect the common, generalist fish species in this study, dietary plasticity may facilitate the invasion of F. virilis in currently unoccupied tributaries. Watershed managers should therefore continue to prevent F. virilis introductions into currently unoccupied tributaries to prevent potential negative effects on sensitive native fish species.