Lake food webs: Species invasion progressively disrupts the trophic structure of native food webs
Creators: name: Wainright, Charles
name: Muhlfeld, Clint
name: Elser, James
name: Bourret, Samuel
name: Devlin, Shawn
Publisher: name: Zenodo
Year: 2021
Subjects: bull trout
C and N stable isotopes
stable isotope ecology
lake ecosystems
aquatic invasive species

description: Species invasions can have substantial impacts on native species and ecosystems, with important consequences for biodiversity. How these disturbances drive changes in the trophic structure of native food webs through time is poorly understood. Here, we quantify trophic disruption in freshwater food webs to invasion by an apex fish predator, lake trout, using an extensive stable isotope dataset across a natural gradient of uninvaded and invaded lakes in the northern Rocky Mountains, USA. Lake trout invasion increased fish diet variability (trophic dispersion), displaced native fishes from their reference diets (trophic displacement), and reorganized macroinvertebrate communities, indicating strong food web disruption. Trophic dispersion was greatest 25-50 years after colonization and dissipated as food webs stabilized in later stages of invasion (&gt;50 years). For the native apex predator, bull trout, trophic dispersion preceded trophic displacement, leading to their functional loss in late-invasion food webs. Our results demonstrate how invasive species progressively disrupt native food webs via trophic dispersion and displacement, ultimately yielding biological communities strongly divergent from those in uninvaded ecosystems.
description: This repository holds scripts, datasets, and an Rmarkdown for "Species invasion progressively disrupts the trophic structure of native food webs."<br> The contents of this repository are named and briefly described below. The Rmarkdown provides a step-by-step explanation of the analyses in this project.<br> rmd_food_web_10012020_14.html (Rmarkdown html document explaining unified R script's analyses with code chunks and figures)<br> rmd_food_web_10012020_14.rmd (the R file to produce the html version of the Rmarkdown) The unified R script is annotated, calls the study's main isotope dataset, and includes code to reproduce analyses.<br> unified_food_web_script_08242020_29.R (study's main R script)<br> isotope_dataset_MT_lake_trout_food_web_11072020.csv (study's main dataset) Besides the main isotope dataset, the unified R script also calls three other datasets:<br> transect_data.csv (for NMDS ordinations)<br> lake_data.csv (for NMDS ordinations)<br> conversion_09302020_2.csv (for binomial regression) The lake trout diet MixSIAR script is separate from the unified script.<br> lake_trout_mixsiar_model_08272020_3.R (lake trout diet R script)<br> lake_trout_mix.csv (MixSIAR mixtures; this is called in the lake trout diet R script)<br> lake_trout_source2.csv (MixSIAR sources; this is called in the lake trout diet R script)<br> lake_trout_trophic_discrimination2.csv (MixSIAR discrimination factors; this is called in the lake trout diet R script) Two other scripts make multi-panel plots from the contents of the unified food web script.<br> trophic_disruption_superplot_11132020.R (combines trophic dispersion and displacement results into one figure)<br> script_superplot_overlap_conversion_11132020_3.R (combines diet overlap and conversion results into one figure) Funding provided by: U.S. Geological Survey<br>Crossref Funder Registry ID: http://dx.doi.org/10.13039/100000203<br>Award Number: Funding provided by: Flathead Lake Biological Station*<br>Crossref Funder Registry ID: <br>Award Number: Funding provided by: Flathead Lake Biological Station<br>Crossref Funder Registry ID:
description: <b><i>Study systems. </i></b>The study area consisted of nine natural lakes and one reservoir (collectively called "lakes'' throughout this article) west of the Continental Divide in northwestern Montana, USA (SI Table 2). These oligotrophic, dimictic lakes are in forested and undeveloped watersheds on U.S. public lands, like State and National Forests and Parks. Study lakes were classified into three categories based on their history of lake trout invasion: 1) reference, 2) mid-invasion (i.e., middle), and 3) late-invasion (SI Table 2). Reference lakes have a native fish assemblage and have no lake trout (current conversion = 0). Mid-invasion lakes have sympatric bull trout and lake trout populations and current (i.e., 2019 or most recent available) conversion between 0.4 and 0.8. Late-invasion lakes also have sympatric bull trout and lake trout populations but have current conversion values greater than 0.8. 'Conversion' is analogous to lake trout dominance. <b><i>Food web sampling. </i></b>All samples were collected between June and October in 2017, 2018, and 2019. Fish were collected with sinking and floating mono- and multifilament gill nets, littoral fyke nets, benthic hoop nets, hook and line, and backpack electrofishing concurrent with U.S. National Park Service and Montana Fish Wildlife and Parks fisheries surveys. Gill nets were 38 m long by 2 m deep panels of 10 to 100 mm bar mesh. Fyke nets had 8 m leads and 4 m hoop sections with one 75 mm vertical trapping pane, one 90 mm throat, and black 6 mm stretch mesh. Benthic hoop nets were 4 m long with two 90 mm throats and black 6 mm stretch mesh. Fyke and hoop nets were deployed in twelve-hour increments. Electrofishing was conducted in shallow water along lake shores using a Smith-Root LR-24 (Smith-Root, Inc. Vancouver, WA). Animal (fish) sampling was conducted by U.S. National Park Service and Montana Fish Wildlife and Parks management agencies during routine monitoring surveys in accordance with agency animal use and care protocols. Bull trout collections were authorized with a special collection permit (Section 10) issued by the U.S. Fish and Wildlife Service. All collected fish were identified to species and measured for length and weight (total length, mm; wet weight, g). A subsample of collected fish were biopsied for stable isotope (δ<sup>15</sup>N and δ<sup>13</sup>C) analysis. Only bull trout and lake trout presumed to be piscivorous (total length ≥ 200 mm) (31) were biopsied, which controlled for ontogenetic diet shifts (SI Figure 12, SI Table 9). A 4-mm diameter biopsy punch (Integra Miltex 336; Integra Life Sciences, Princeton, NJ, USA) was used to extract a sample of each fish's dorsal white muscle (32). Muscle samples were stored in 100% industrial ethanol (95% ethanol, 5% methanol) while afield and stored in a -10⁰C freezer for later processing. Littoral macroinvertebrates were collected by kick-netting with a 500-µm D-net at 0.5, 1, and 1.5 m depths at seven transects in each lake. Profundal macroinvertebrates were collected with a grab sampler dredge from depths exceeding twice the Secchi maximum at each lake and filtered through a 500-µm D-net. Macroinvertebrates were depurated in clean lake water, identified to family, and preserved in 100% industrial ethanol for isotope sample preparation. Bulk zooplankton were collected using repeated tows at 10 m depth of a 100-µm tow net in the pelagic epilimnion of each lake. Epilithic periphyton was collected from littoral rocks at 0.5 m depth of each lake using a brush. Periphyton and zooplankton were stored in filtered lake water, kept cold while afield, and stored in a -10⁰C freezer for later processing. <b><i>Laboratory analyses. </i></b>Fish muscle samples were further subsampled to generate a representative and comparable analytical dataset for each lake. The analytical subset of fish muscle tissue is as follows: 1) all available bull trout; 2) 10 lake trout; and 3) five of all other sampled fish species. The analytical subset was selected before sending samples to the stable isotope lab and controls for bias that could result from unequal sample sizes among species within a functional group and the covariance between species richness and invasion category (SI Table 8). All samples were dried in a 60⁰C oven for 72 hours and homogenized to a powder using a mortar and pestle (33). Dried samples were weighed (animal tissue: 1 mg ± 0.1 mg; periphyton: 10 mg ± 0.1 mg) and loaded into tin cups (5x9 mm, Costech 41077). Stable isotope cap combustion was conducted at University of California at Davis on a PDZ Europa ANCA-GSL interfaced to a PDZ 20-20 Europa Scientific isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK). Absolute isotopic ratios were expressed in standard delta "δ" notation relative to Vienna PeeDee Belamnite (δ<sup>13</sup>C) and atmospheric nitrogen (δ<sup>15</sup>N) (34). Stable isotope analysis using ratios of nitrogen (<sup>15</sup>N:<sup>14</sup>N; δ<sup>15</sup>N) and carbon (<sup>13</sup>C:<sup>12</sup>C; δ<sup>13</sup>C) is a quantitative way to infer two main parts of food web structure: trophic position and energy sources for secondary production (20). Nitrogen stable isotope concentration, δ<sup>15</sup>N, is used to infer trophic position because <sup>15</sup>N is enriched (approximately 3-4 ‰) in predators relative to their prey (20, 35). In contrast, carbon stable isotope concentration, δ<sup>13</sup>C, changes little between predators and prey (&lt;1 ‰) but benthic algae are <sup>13</sup>C enriched relative to phytoplankton (20, 36). Therefore, δ<sup>13</sup>C is used to determine the source(s) of organic carbon used in secondary production (20, 36). Absolute (i.e., empirical) δ<sup>15</sup>N and δ<sup>13</sup>C values vary as isotopic baselines change seasonally and among lakes (37, 38). To account for spatial and temporal isotopic variability among our study lakes, we corrected empirical δ<sup>15</sup>N and δ<sup>13</sup>C using the formula δX<sub>corrected </sub>= δX<sub>empirical </sub>− δX<sub>baseline</sub> described by Hobson et al. (39). We used our study system's most abundant mayflies (Ephemeroptera: Baetidae, Leptophlebiidae, Heptageniidae, and Ephemerellidae) as lake-specific primary consumer isotopic baselines (SI Table 2). We chose mayflies as isotopic baselines because they were abundant in our study lakes, mayflies overlapped each other in δ<sup>15</sup>N and δ<sup>13</sup>C (mean ± standard error), and mayfly δ<sup>15</sup>N and δ<sup>13</sup>C were less variable than other primary consumers or periphyton. We selected mayflies as our trophic baseline because they showed the lowest within-lake δ<sup>13</sup>C and δ<sup>15</sup>N variability; a common protocol for controlling lake effects in isotopic food web studies (40). To arrive at this baseline selection, we analyzed within-lake δ<sup>13</sup>C and δ<sup>15</sup>N variability selected six primary consumers (mayflies, zooplankton, caddisflies, hydrachnidia, snails, and profundal chironomids) and periphyton, hereafter called trophic baseline 'candidates.' These candidate taxa are ideal for trophic baseline corrections due to low tissue turnover (37), their abundance across our study lakes, and taxa specific overlap in δ<sup>13</sup>C and δ<sup>15</sup>N, which indicates these candidates had consistent δ<sup>13</sup>C and δ<sup>15</sup>N through space and time (35). We determined the within-lake isotopic variability for each baseline candidate by calculating the within-lake δ<sup>13</sup>C and δ<sup>15</sup>N standard error for each candidate taxa (SI Table 6) and among-lake δ<sup>13</sup>C and δ<sup>15</sup>N variability by averaging the within-lake δ<sup>13</sup>C and δ<sup>15</sup>N for each candidate (SI Table 7). Mayflies were most suited taxa as a trophic baseline because they had the lowest average within-lake δ<sup>13</sup>C and δ<sup>15</sup>N (SI Table 7) and mayfly consumer δ<sup>13</sup>C demonstrates remarkable within- and among- lake correlations between both basal resource δ<sup>13</sup>C and secondary consumers (41). <b><i>Data analyses and food web quantification</i></b><b>. </b>Fish species were assigned to five functional groups for analyses: 1) bull trout, 2) lake trout, 3) littoral forage fish, 4) generalist fish, and 5) pelagic forage fish (SI Table 3). Functional groups (i.e., littoral and pelagic forage fish and generalist fish) aggregated presumed mesopredator fishes based on habitat and feeding patterns (42, 43) and trophic position relative to lake trout and bull trout (44). Littoral forage fish, like redside shiners (<i>Richardsonius balteatus</i>), occupy nearshore habitat (43). Generalist fish, like cutthroat trout (<i>Oncorhynchus clarkii lewisi</i>), may occupy either nearshore or offshore habitats (43). Pelagic forage fish, such as kokanee (<i>Oncorhynchus nerka</i>), occupy offshore habitat (43). We used zooplankton and periphyton δ<sup>13</sup>C end members to make inferences on littoral vs. pelagic carbon utilization. Two-end member δ<sup>13</sup>C models are a common method by which comparative food web studies infer changes in fish δ<sup>13</sup>C because δ<sup>13</sup>C patterns are conserved through trophic transfers (20). Refer to SI Figures 6 and 7 for fish δ<sup>13</sup>C inferences (i.e., manuscript figures 2b and 2c) with end-member data presented. Bulk periphyton is an excellent endmember for determining littoral derived energy because the potential range and variability of littoral basal resource δ<sup>13</sup>C induced by PAR driven variability in rates of benthic primary production (BPPr) are conserved and integrated. PAR is attenuated in lake water as a function of depth within lakes, as well as a function of subsurface depth within the periphyton mat assemblages themselves (41). Many different functional groups and taxa may rely completely on periphyton autochthony while exhibiting substantial variation in δ<sup>13</sup>C induced by the depth-specific, light driven patterns in BPPr. It is absolutely essential to consider the relationship between bulk periphyton δ<sup>13</sup>C and BPPr when selecting the proper littoral endmember. Failure to do so will result in misinterpretation and poor measurements of trophic structure due to improper endmember selection and overestimation of pelagic carbon sedimentation (41). For instance, snails graze on surficial periphyton mats in high light conditions and consume the most <sup>13</sup>C enriched, most productive periphyton. In contrast, chironomids are found deeper in the periphyton mat at the exact same locale, where less actualized irradiance leads to less productive, less enriched δ<sup> 13</sup>C periphyton. In this case, both snail and chironomids were 100% reliant on periphyton but had divergent δ<sup>13</sup>C signatures and this divergence is typically interpreted as mixing of non-littoral sources (i.e., pelagic, terrestrial, and chemogenic). Poor endmember selection propagates the δ<sup>13</sup>C divergence/source mixing misunderstanding by strongly skewing consumer carbon use estimates towards non-littoral sources by as much as 30% (41). Without considering how benthic consumer partitioning of proximally shared resource(s) results in consumer isotopic divergence without source mixing, the spatiotemporal patterns in food web structure, littoral-pelagic coupling, and the isotopic findings presented in this manuscript regarding the mechanisms of species invasion and the trophic disruption hypothesis would either missed or misinterpreted. All data analyses were completed in R version 3.5.1 (45) and RStudio version 1.1.456 (46). All figures were produced with R package "ggplot2" version 3.1.0 (47). Linear mixed effects models from R package "lme4" version 1.1-19 (48) were used to examine magnitude and direction of changes in δ<sup>13</sup>C and δ<sup>15</sup>N to infer trophic displacement between lake trout invasion categories. Model fit and residual normality were confirmed using residual plots. R packages "lmerTest" version 3.1-1 (49) and "effects" version 4.1-4 (50) were used to test for statistical significance of δ<sup>13</sup>C and δ<sup>15</sup>N changes among invasion categories. Baseline-correcting lakes and using a mixed-effects models accounted for 'lake effects' from among-lake isotopic differences (SI Figure 10 and 11). Bayesian standard ellipse area (SEA.b) was calculated to quantify isotopic niche area and infer trophic dispersion using R package "SIBER" version 2.1.4 (51). SIBER's Markov chain Monte Carlo (mcmc) was provided baseline-corrected δ<sup>13</sup>C and δ<sup>15</sup>N for each fish group and uninformative prior probabilities for population mean and variance (µ and σ<sup>2</sup>). Then SIBER's mcmc then produced a posterior distribution for population mean and variance. Next, SIBER generated a distribution of isotope ellipses (51) based on the posterior distributions of population mean and variance. Finally, SIBER was used to calculate SEA.b for posterior-derived ellipses (51). As diet specificity increases, ellipse area decreases so SEA.b is a probabilistic approach to inferring diet specificity (51). Niche overlap probabilities were modelled with R package "nicheROVER" version 1.0 (52, 53). Pairwise comparison of isotopic niche overlap describes directionality of overlap, which is useful for examining the likelihood of competitive exclusion (52). Asymmetric overlap, where one species is likely to be in another species' isotopic niche, but the opposite is not likely, can suggest competition (52). Conversely, symmetric overlap, where both species are likely to exist in each other's isotopic niche, can suggest resource partitioning (52). R package "MixSIAR" version 3.1.10 (54, 55) was used to model lake trout diet in mid- and late-invasion lakes. MixSIAR's mcmc was provided uninformative prior probabilities for prey consumption, baseline-corrected lake trout δ<sup>13</sup>C and δ<sup>15</sup>N, baseline-corrected prey fish mean and standard deviation δ<sup>13</sup>C and δ<sup>15</sup>N, and Post's (35) trophic discrimination factors. MixSIAR's mcmc then produced a posterior probability distribution of proportions of each prey category, bull trout and other fish. Model convergence was evaluated with traceplots. Standardized gill net survey data and binomial linear regression were used to substitute space for time in our study area. Equation 1 (SI Table 5) was used to calculate each lake's 2019 (or most recent) conversion (adapted from Catford et al. (7)). Binomial generalized linear modeling (GLM) was used to regress conversion based on historical gill net survey data (SI Table 5, Equation 2). GLM fit was evaluated with residual plots. Model r-squared was calculated using R package "rsq" version 2.0 (56). R package "vegan" version 2.5-3 (57) was used to calculate non-metric dimensional scaling (NMDS) ordinations to evaluate macroinvertebrate community similarity within and among lake trout invasion timesteps. Macroinvertebrate communities were ordinated by sampling transect because transect is the lowest aggregation at which sampling design allowed for sample independence. Next, NMDS scores were calculated in vegan to quantify macroinvertebrate community similarity each transect, NMDS scores were grouped by invasion category (i.e., reference, mid-, and late-invasion) and NMDS scores were plotted with 95% confidence interval ellipses. Finally, community similarity among invasion categories was tested using permutational multivariate analysis of variance (permANOVA). Since reservoirs experiencing annual water drawdowns have different littoral macroinvertebrate communities than natural lakes (58), Hungry Horse Reservoir was excluded from these ordinations.