Modeling the trophic impacts of invasive zooplankton in a highly invaded river

The lower Columbia River (Washington and Oregon, USA) has been heavily invaded by a large number of planktonic organisms including the invasive copepod Pseudodiaptomus forbesi and the planktonic juveniles of the invasive clam, Corbicula fluminea. In order to assess the ecological impacts of these highly abundant invaders, we developed a multivariate auto-regressive (MAR) model of food web dynamics based upon a 12-year time-series of plankton community and environmental data from the Columbia River. Our model results indicate that plankton communities in the lower Columbia River are strongly impacted by the copepod P. forbesi at multiple trophic levels. We observed different ecological effects across different life stages of P. forbesi, with nauplii negatively impacting ciliates and autotrophs, and copepodite stages negatively impacting Daphnia and cyclopoid copepods. Although juvenile C. fluminea were highly abundant in the summer and autumn of each year, our best fit MAR model did not show significant C. fluminea impacts. Our results illustrate the strong ecological impact that some zooplankton invaders may cause within rivers and estuarine systems, and highlight the need for further research on the feeding ecology of the planktonic life-stage of C. fluminea. Overall, our study demonstrates the manner in which long-term, high resolution data sets can be used to better understand the ecological impacts of invasive species among complex and highly dynamic communities.


Introduction
Aquatic invasive species are an increasingly common aspect of ecological change in marine and freshwater systems. Such invasions often originate from the global transport of zooplankton in the ballast water of commercial shipping vessels [1][2][3], although the overland transport of recreational boats, movement of fishing gear, and natural vectors such as the movement of aquatic birds play a role as well [4]. Aquatic invaders have in some instances catalyzed wholescale ecosystem shifts, as in the case of Dreissenid zebra and quagga mussels in the North American Great Lakes [5][6][7] and the invasion of the Black Sea by the ctenophore Mnemiopsis leidyi [8,9]. Aquatic invaders have also been shown to inflict heavy damages to power generation facilities, imperil drinking water supplies, interfere with recreation and tourism, and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The clam, Corbicula fluminea, is one of the most widely distributed invasive bivalves across the northern hemisphere, and considered a species of ecological concern in many watersheds across North America and Western Europe [31,42]. Dispersal of C. fluminea primarily occurs via the release of planktonic juveniles that resemble miniature adults with well-developed shell, adductor muscles, foot, gills and digestive system [43]. These planktonic juveniles are suspended by turbulent water currents and transported downstream before settling to the benthos where they reach maturity in 3-6 months [43]. Adult C. fluminea exhibit relatively nonselective filter feeding across a wide range of microplankton taxa [44], with some evidence for selection against cyanobacteria [45]. Despite the extensive body of literature devoted to the adult stage of C. fluminea [43][44][45][46], the feeding habits of the planktonic juvenile stage are essentially unknown. C. fluminea is one of the most numerous planktonic organisms in the CR during summer and autumn [30], but it is unclear if this juvenile stage interacts with other members of the plankton community in any significant manner. We hypothesize that young C. fluminea exhibit similarly non-selective feeding, but we are uncertain if the juveniles feed at rates sufficient to impact standing stocks of microplankton prey.
In order to better understand the trophic relationships and community impacts of P. forbesi and juvenile C. fluminea, we produced a multivariate auto-regressive (MAR) model of plankton community dynamics using a 12-year series of monthly plankton samples collected from the lower CR. Our MAR model is based upon a set of regression models (one for each species contained in the dataset) which estimates the strength and directionality of interspecific interactions among species. MAR modeling has been previously employed to elucidate food-web structure for a number of different plankton communities [47][48][49][50][51][52][53]. For example, a MAR model generated from a 60-year time series of plankton populations in Lake Baikal (Russia) was instrumental in understanding the causes and consequences of a long-term increase in cladoceran biomass in the lake, and elucidating the relationship between increased temperatures and shifts in community composition [48]. Likewise, MAR analysis of a 33-year time series of plankton community abundances from Lake Washington (Washington, USA) was used to validate conceptual food web models, and to identify previously uncharacterized trophic interactions [47].
Our study examines the role of invasive zooplankton in the trophic web of plankton communities of the lower CR (Washington and Oregon, USA). Our main objectives were to evaluate the ecological impacts of P. forbesi and juvenile C. fluminea, and to elucidate the key interactions among both native and invasive members of the community. Furthermore, our analysis aims to characterize the trophic habits of the juvenile planktonic life stage of C. fluminea, which is very much understudied relative to later benthic life stages.

Study site
The Columbia River (CR) Basin drains an area of 669,300 km 2 which includes portions of seven U.S. states and two Canadian provinces [54]. The CR extends for 1,954 km and discharges on average 224 billion m 3 of annual outflow into the Pacific Ocean [54]. The hydrology of the river is heavily influenced by the timing of spring and early summer snowmelt, although more than 200 impoundments on the river have greatly reduced seasonal variation in flow relative to historical patterns [55]. Large commercial fisheries for anadromous fishes such as steelhead trout, chinook, coho, chum, and sockeye salmon exist on the CR, with multiple fish stocks protected under the U.S. Endangered Species Act [56].
We collected plankton samples and environmental data from a single fixed location in the CR, as part of a larger ongoing program of plankton monitoring across the lower river and estuary (Fig 1). This site is located at a pier in Vancouver, WA, USA (45.6222˚N, 122.6772˚W), approximately 170 river km upstream from the mouth of the CR and 64 km downstream from the Bonneville dam, which is the furthest downstream impoundment on the river. Hydrological conditions at this site are characterized by high flow with no thermal stratification, and water depths ranging seasonally from 8-11 m [30]. This site is upstream of any saltwater intrusion but remains tidally influenced.

Sample collection
We collected a suite of plankton samples and environmental measurements once every month across an uninterrupted 12-year period spanning January 2005-December 2016, yielding a total of 144 continuous monthly samples. On each sampling date, we collected three replicate zooplankton samples via vertical net tows from a depth of 1 meter above river bottom to the surface using a 73-μm mesh, 0.5-m diameter ring net with attached flowmeter (General Oceanics). Zooplankton samples were fixed in a 10% buffered formalin solution in the field for later taxonomic processing. Water samples for microplankton community analysis were collected via triplicate bucket samples of surface water and preserved in 5% acid Lugol's solution. Surface water samples were considered as representative of the water column due to strong vertical mixing at this location. We measured salinity, temperature, and water clarity on each sampling date using a YSI 85 temperature/salinity probe (YSI Incorporated) and a Secchi disk. Temperature and salinity data were collected at two-meter intervals from surface to bottom. We additionally collected triplicate samples of surface water for chlorophyll a analysis. These surface water samples were kept chilled in 70 mL opaque bottles until subsequent filtration and fluorometric measurement of chlorophyll a in the laboratory via the acidification method [57] using a Turner model 10-AU fluorometer (Turner Designs). As water samples were collected from a public location and did not involve vertebrate or otherwise regulated species, no special permits were associated with this project.

Taxonomic processing of samples
A minimum of 200 non-naupliar organisms were subsampled from aliquots of each zooplankton sample using a Stempel pipette, which were examined using a Leica MZ6 stereomicroscope (Leica Microsystems) at 40X magnification. Specimens were identified to the lowest possible taxonomic rank using Thorp and Covich [35], with most rotifers and microcrustaceans identified to the genus or species level. We converted counts of individual taxa to density (individuals m -3 ) by multiplying each count by the ratio of sample volume to subsample volume, and then dividing by the total volume of water filtered. Two replicate samples were processed for each sampling date.
Microplankton samples (i.e. planktonic protists, eukaryotic algae, cyanobacteria) were taxonomically processed via 1-10 ml subsamples of the Lugol's preserved water samples. Subsamples were settled overnight in Utermohl chambers, and the chambers examined using an Olympus CK-40 inverted microscope at 200X (numerical aperture = 0.75). All individuals were identified to genus (and species when possible) using Wehr et al. [58] and sized using an ocular micrometer in order to calculate biovolume and carbon biomass based on geometric shape [59,60].

Multivariate autoregressive (MAR) modeling of ecological interactions
The MAR model is a system of p linear equations describing the abundances for each species in the community, with p equal to the number of taxa in the model and q equal to the number of environmental covariates. In matrix form, the MAR model is written as follows: where x t is the p × 1 vector of log abundances for each of the p species at time t, B is a p × p interaction matrix whose elements b ij (the interaction coefficients) describe the effect of the density of species j on the per capita growth rate of species i, a is the p × 1 vector of a values (the intrinsic rates of increase) for each species, C is the p × q matrix whose elements c ij describe the effect of covariate j on species i, and u t-1 is the q × 1 vector of covariate values at time t-1. The vector of process errors w t is assumed to be drawn from a multivariate normal distribution with a mean of 0 and covariance matrix S. See Ives et al. [61] and Hampton et al. [49] for more complete overviews of the MAR approach.
Our complete time series contained more than 100 mesozooplankton and microplankton taxa, which is intractable for the purposes of statistical modeling. We therefore aggregated our data into a set of functional and/or higher-level taxonomic groups based upon our ecological questions of interest (Table 1). Prior to model input, all data were normalized within season and standardized across variables and all species abundances were log(x+1) transformed [61]. Invasive taxa (P. forbesi and C. fluminea) were modeled as extrinsic forcings on the native plankton community through the inclusion of these taxa as covariates in the MAR-but on the same time lag (t-1) as the native species variables [61]. We did not include the invasive B. coregoni in the model, as it occurs irregularly and at typically low abundances in our system. In order to disentangle interspecific interactions from broad patterns of seasonal succession, abundance data were normalized (i.e. "de-seasoned") by subtracting the raw abundance values from the monthly mean abundance for each taxon. Standardization across variables was then achieved through conversion of these de-seasoned values to dimensionless z-scores.
When data are prepared for MAR input, linear interpolation is sometimes used to fill missing values, provided that the gaps are short and infrequent [49]. This process was used to fill one instance of missing microplankton data (March 2016) which was lost due to human error. Model validation was conducted through visual inspection of residuals, and evaluation of parameter estimate convergence times. The robustness of model results was further assessed through iterative permutation of the model structure.

Phenological patterns of the system
Water temperatures ranged between 2.8-23.5˚C with average minimum temperatures of 4.5˚C typically occurring in January of each year and average maximum temperatures of 21.2˚C occurring in August (Fig 2A). River outflow, as measured at the nearby upstream Bonneville Dam (Fig 2B), ranged between 2.38-13.  (Fig 2C). Chlorophyll a concentration ranged between 0.01-23.3 μg/L, with a mean annual peak of 11.8 μg L -1 occurring in each spring (Fig 2D). Total microplankton abundance peaked in the spring of each year, typically reaching a highly variable maximum in April (mean of 5590 individuals mL -1 ), followed by more or less continuous decline until the following spring ( Fig 2E). The microplankton were largely comprised of single celled photosynthetic organisms such as diatoms, cyanobacteria, and green algae, and to a much lesser extent mixotrophic/heterotrophic ciliates and dinoflagellates (Table 1).
Total zooplankton abundance exhibited two distinct seasonal peaks-one in the spring and one in late summer / early autumn of each year (Fig 2F). The smaller spring zooplankton peak was largely comprised of large rotifers of the genera Asplanchna and Brachionus, and to a lesser extent native cladocerans and cyclopoid copepods. Note that the mesh size of our zooplankton

PLOS ONE
collection net (73-μm) may have allowed smaller rotifer species to pass through the net, especially at the youngest life stages. The larger fall zooplankton peak was almost entirely comprised of the invasive copepod Pseudodiaptomus forbesi and juveniles of the invasive Asian clam, Corbicula fluminea. During these two seasonal zooplankton peaks, total abundance typically reached or exceeded 2,500 individuals m -3 , and between peaks, as few as 2 individuals m -3 were observed, with the smallest abundances typically observed during winter months. The zooplankton community experienced an annual cycle of oscillation between states of extreme numerical dominance by invasive taxa and states comprised almost exclusively by native taxa (Fig 3D). Three invasive taxa were observed in total: the calanoid copepod Pseudodiaptomus forbesi (Fig 3A), larvae of the Asian clam Corbicula fluminea (Fig 3B), and the cladoceran Bosmina (Eubosmina) coregoni ( Fig 3C). P. forbesi and C. fluminea were consistently observed in far greater abundances than B. coregoni (Table 1, Fig 3). The onset of the period dominated by invasive taxa typically occurred in June or July, peaked in August or September, and was followed by a transition to a relatively depauperate winter community in November.

PLOS ONE
During periods of peak invader abundance, these taxa typically comprised more than 90% of total zooplankton abundance (Fig 3D), with values greater than 99% observed on multiple occasions.

Ecological impacts of invasive zooplankton taxa
Our best-fit MAR model contained significant interactions among native members of the community, as well as invasive members acting on native members of the community (Fig 4;  Table 2). Additionally, our model contained several significant interactions between water temperature and native zooplankton. The model contained one significant negative interaction between native plankton taxa, and four significant negative interactions between invasive and native taxa (Fig 4). P. forbesi nauplii (larval stage) and copepodites (juveniles and adults)

PLOS ONE
were treated as separate taxa in the model due to their high degree of morphological and ecological dissimilarity [36,37]. Model results support this treatment, as interactions associated with P. forbesi nauplii were distinct from those of copepodites. P. forbesi nauplii showed a strong negative effect on autotrophic microplankton (diatoms, flagellates, green algae, and cyanobacteria) and ciliates. P. forbesi copepodites were associated with strong negative effects on native Daphnia and native cyclopoid copepods. The single significant negative interaction between two native zooplankton groups contained in our model was a strong negative effect of Daphnia on native bosminid cladocerans. Our model also contained several positive interactions between members of the microplankton community (Fig 4). Autotrophic microplankton exhibited a strong positive effect on both Asplanchna rotifers and ciliates, and ciliates exhibited a strong positive effect on autotrophic microplankton.

Discussion
Prior studies have shown that the Columbia River and estuary are highly invaded by several species of zooplankton [26,27,30,32,65]. Among these invaders, the calanoid copepod P. forbesi is clearly the most abundant taxon, with community composition approaching a monoculture of this species in the autumn of some years [30]. While prior studies have shown that P. forbesi is both widespread and seasonally abundant in the CR (as well as the San Francisco Estuary), the ecological impacts of this invader have largely remained an unresolved question [26,28,30,41,66].
Our investigation has shown a clear pattern of community-level effects arising from the invasion of Pseudodiaptomus forbesi. Our results indicate that the various life stages of P. forbesi exhibit negative effects on native cyclopoid copepods, Daphnia, ciliates, and autotrophs (diatoms, flagellates, green algae, and cyanobacteria). These results accord with laboratorybased studies which showed that P. forbesi is capable of feeding on diatoms, ciliates, flagellates,  Maximum likelihood parameter estimates and 95% confidence intervals for significant MAR terms. Rows with the same predictor and response variable represent density dependent population growth. and dinoflagellates, and may exhibit a strong preference for diatoms and ciliates [40,41]. We therefore interpret the negative effect of P. forbesi on autotrophic microplankton and ciliates in our MAR model as a direct effect of grazing. P. forbesi are not known to predate upon larger metazoans, and thus we interpret the negative effects of P. forbesi on native cyclopoid copepods and Daphnia as a consequence of competition for microplankton prey resources. It remains unclear what functional traits are responsible for this apparent competitive advantage over native members of the community, although laboratory studies indicate that lower rates of fish predation relative to native Columbia River zooplankton may play some role [67].
Our study highlights the need for a better understanding of the role that juvenile C. fluminea play in planktonic food webs. While C. fluminea spend only a brief portion of their life in this life stage, they are typically among the most abundant zooplankton in the CR for several months of each year (Fig 4). Our MAR model did not show significant interactions between juvenile C. fluminea and native members of the community, although it did contain some non-significant interactions that were high in magnitude. These results may be interpreted in several different ways. On one hand, our model may indicate that juvenile C. fluminea feed at low rates, and thus exert little trophic influence on CR plankton communities. Alternatively, the grazing effects of juvenile C. fluminea might be more taxon-specific than we can detect in our relatively course grouping of microplankton species. These hypotheses might be discriminated through laboratory-based studies that directly examine the feeding habits of C. fluminea juveniles, or DNA and/or microscopy based-analysis of gut contents. In future years, the inclusion of additional data from the ongoing CR time-series may also reduce uncertainty around our estimated parameters and allow for more precise characterization of some ecological interactions not resolved in the current iteration. Nonetheless, we can presently conclude that juvenile C. fluminea do not exhibit trophic impacts that match the magnitude and breadth of the invasive P. forbesi.
We did not attempt to estimate the ecological impacts of the invasive cladoceran Bosmina (Eubosmina) coregoni in our MAR model, as it was uncommonly observed in the lower CR. The species is a recent arrival to the Pacific coast of North America, having been first detected in the CR in late 2006 [30,68]. It is unclear if the sporadic appearance of B. coregoni in the lower CR represents the early phases of a nascent invasion, the erratic population dynamics of a failed introduction, or the downstream advection of individuals from an upstream source population. The fact that B. coregoni is more typically associated with lake environments, rather than fast-flowing rivers [68,69], lends some support to the hypothesis that these individuals have been carried downstream from an established population at one or more of the many upstream impoundments of the CR.
Our MAR model also identified several interactions between native members of the plankton community, including a strong negative effect of native Daphnia on native bosminid cladocerans. This result is consistent with previous findings that under certain circumstances, Daphnia may strongly suppress the growth of Bosmina populations by more effectively grazing shared food resources [70,71]. Our model also identifies a strong positive effect of autotrophic plankton (diatoms, flagellates, green algae, and cyanobacteria) on Asplanchna. Asplanchna are relatively large predatory rotifers that typically feed upon smaller rotifers [72,73] and protozoans [74], which in turn graze upon a wide range of autotrophic microplankton and bacterioplankton [35,75,76]. As small-bodied rotifers and protozoans were not included in our MAR model, the apparent positive effect of autotrophic plankton on Asplanchna may be an indirect positive effect via the (unobserved) smaller rotifer and protozoan species that Asplanchna prey upon, rather than a direct positive effect from grazing. Our model also showed a biologically plausible positive effect of autotrophic plankton on ciliates, as many species of ciliated protozoa (most notably oligotrichs) are known to graze heavily on phytoplankton [77]. The meaning of the reciprocal positive interactions between ciliate and autotrophic plankton is unclear, and may represent trophic dynamics within the microplankton community that requires analysis at finer taxonomic or temporal scales. One possible interpretation is that ciliates may be selectively consuming one group of autotrophic taxa (e.g. green algae) that is then reducing competition among all autotrophic microplankton for nutrients or light, thus resulting in an overall increase of autotrophic microplankton.
Our best-fit MAR model contained a small number of significant interactions relative to the number of estimated parameters. This is a typical characteristic of MAR models that are constructed from ecological datasets [48,49,51], and does not indicate that CR plankton dynamics are governed by few interactions. Rather, our results highlight the interspecific interactions which are sufficiently strong and consistent that they can be isolated from our highly dynamic time-series data. Our model contained an additional number of strong (but non-significant) interactions which were eliminated after bootstrap-based estimation of confidence intervals. This result was anticipated given the relatively short duration of our time-series, relative to the multi-decadal datasets from which MAR models of plankton food webs are typically constructed [48,50,52,53]. Our best-fit model also showed very few significant effects of temperature, which was not unexpected given the wide confidence intervals surrounding many of our parameter estimates. We nevertheless chose to retain temperature as a variable in our MAR model because its inclusion improved model performance, but we would caution against over-interpretation of the lack of significant temperature effects retained in our final MAR model. Our group continues to collect monthly samples from the lower CR, and this additional data may substantially reduce the uncertainty surrounding some of our parameter estimates in future iterations of our model.
There are important limitations to the questions informed by our model results. Foremost is the fact that our model does not capture the dynamics of the entire Columbia River ecosystem, rather the focus is on the interspecific interactions within the CR plankton community. For example, we have not attempted to estimate the effects of predation by planktivorous fish or larger invertebrates. This choice was driven by both a desire to maintain a specific focus on the effects of P. forbesi and C. fluminea on native members of the plankton community, and a lack of predator data on the appropriate temporal scale for this modeling approach. It is difficult to predict the potential strength of top-down impacts of planktivory from fishes and larger invertebrates (such as adult C. fluminea) as the effect would depend acutely on both the abundance and timing of co-located predator and prey given the seasonal booms and busts in zooplankton abundance (Fig 2). Several studies have examined specific fish-zooplankton interactions in the Columbia River [67,78,79], but consistent monitoring of fish populations has been restricted to only a few species of game fishes [e.g. 80,81], thus making it difficult to rigorously infer community-wide changes in predation pressure across the period of study. Similarly, nutrient concentrations in the river likely play an important bottom-up role in regulating the growth of many of the modeled microplankton species, but were not addressed in these models.
An additional point to bear in mind is that our approach is based upon examining patterns of covariance, from which we attempt to infer causal mechanisms based upon the relevant ecological literature. It is possible that some of the patterns that we have observed arise from a shared correlation with an unmeasured variable, rather than a direct causal relationship. In regards to purely temporal correlations, we have controlled for seasonality by normalizing and standardizing our raw data prior to model input, but hidden correlations from other unmeasured variables are possible. Finally, the taxonomic groupings that we employed, especially in regards to microplankton prey, may have obscured any highly selective feeding habits that did not strongly affect overall group abundance.
In conclusion, our results illustrate the strong ecological impact that zooplankton invaders may cause within rivers and estuarine systems. Our model results indicate that plankton communities in the lower CR are strongly affected by the invasive copepod Pseudodiaptomus forbesi at multiple trophic levels. We observed different ecological impacts across different life stages of P. forbesi, with nauplii negatively impacting ciliates and autotrophs, and copepodite stages negatively impacting Daphnia and cyclopoid copepods. Our model does not resolve any ecological interactions between juveniles of C. fluminea and the rest of the CR zooplankton community, and highlights the need for additional research on the planktonic phase of this widely distributed and abundant invader. Overall, our study demonstrates the manner in which long-term, high resolution data sets can be used to better understand the ecological impacts of invasive species among complex and highly dynamic communities.