DOI: 10.5061/dryad.z34tmpgdt
title: Lake food webs: Species invasion progressively disrupts the trophic structure of native food webs

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: <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.
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)