Zooplankton invasion on a grand scale: insights from a 20-yr time series across 38 Northeast Pacific estuaries

We present the first comprehensive analysis of the Pacific Northwest estuaries (PNWE) zooplankton time series, which encompasses 38 estuaries distributed across more than 1000 km of the North American Pacific Coast. With observations spanning more than 20 yr, we here examine biogeographic trends among zooplankton communities, patterns of biological invasion across the region, and environmental correlates with dominant native and invasive taxa. Our results show that some estuaries across the region are invaded by multiple zooplankton species and that the geographic extent of invasion is far greater than previously reported for at least five species of copepods: Pseudodiaptomus inopinus, Pseudodiaptomus forbesi, Oithona davisae, Limnoithona sinensis, Sinocalanus doerrii, and the cladoceran Bosmina coregoni. Some of these species appear to be rapidly spreading across the region, while others have occupied a relatively static geographic range for decades. The copepod, P. inopinus, is by far the most abundant and geographically widespread of these invaders, comprising more than 90% of all zooplankton abundance at some sites. We propose that the geographic distribution of these invaders is strongly constrained by geomorphic characteristics that define the salinity and mixing regimes in these estuaries, reflecting the strong role that physical forces play in structuring estuarine zooplank-


Estuarine zooplankton invasions
Zooplankton are a taxonomically diverse assemblage of predominantly microscopic organisms that play a crucial role in freshwater, estuarine, and marine ecosystems. Freshwater and estuarine zooplankton populations tend to be structured across landscapes in a manner determined by the hydrological connectivity of the water bodies that they inhabit. Some species of zooplankton achieve long-range dispersal via stress-resistant eggs carried by floodwaters, aquatic birds, or air currents (C aceres andSoluk 2002, Frisch et al. 2007), but varied anthropogenic activities now facilitate the long-range transport of many zooplankton species.
Of the numerous zooplankton invasions that have been documented in recent decades (Bollens et al. 2002, Dexter andBollens 2019), most have been associated with the discharge of ballast water from commercial shipping vessels (Carlton and Geller 1993, Ruiz et al. 2000, Grosholz 2002. After an initial period of establishment, newly arrived estuarine species may undergo secondary spread to upstream locations (e.g., freshwater-tolerant species such as the copepod Pseudodiaptomus forbesi) or downstream to coastal waters (e.g., salt-tolerant species such as the copepod Oithona davisae; Svetlichny and Hubareva 2014, Emerson et al. 2015). However, the mechanisms and patterns of secondary spread across estuarine systems are not well-understood, because regular monitoring of zooplankton has occurred at only a small number of estuaries, making it difficult to understand invasion processes across regional scales.
For more than two decades, our group has surveyed zooplankton communities across the Pacific Northwest coast of North America, from southern British Columbia, Canada, to northern California, USA, with recent expansion into central and southern California. These surveys have been conducted every four years  by sampling along the upstream salinity gradient at 20-30 estuaries. The surveyed estuaries span a wide range of environmental conditions and land use, from small rural watersheds to large urbanized water bodies (Lee and Brown 2009).
A least 12 nonindigenous zooplankton species have been documented across this region, but most observations have been from two-well studied bodies of water: the San Francisco Estuary and the lower Columbia River. Note that we employ the term "nonindigenous" in reference to any nonnative species, but reserve the term "invasive" for species that have been associated with negative ecological impacts or rapid population expansion (see Davis and Thompson 2000 for a discussion of this terminology). This list of nonindigenous species includes nine copepods (e.g., Acartiella sinensis, Limnoithona sinensis, Limnoithona tetraspina, Oithona davisae, Pseudodiaptomus forbesi, Pseudodiaptomus inopinus, Pseudodiaptomus marinus, Sinocalanus doerrii, and Tortanus dextrilobatus), the cladoceran Bosmina coregoni, and larvae of the clams Corbicula fluminea and Potamocorbula (Carlton et al. 1990, Orsi and Walter 1991, Cohen and Carlton 1998, Orsi and Ohtsuka 1999, Gifford et al. 2007, Winder and Jassby 2011, Bollens et al. 2012, Smits et al. 2013, Hassett et al. 2017). Among these species, 50% occur at both locations. Zebra mussels (Dreissena polymorpha) and quagga mussels (Dreissena bugensis) are not yet present at either site, but quagga mussels have recently arrived in the western United States and may be poised for rapid spread across the region (Wong and Gerstenberger 2011).
We present the first detailed analysis of the PNWE zooplankton time series, in which we examine biogeographic trends among Pacific Northwest zooplankton communities, patterns of biological invasion across the region, and environmental correlates of dominant native and invasive taxa. We also examine the salinity and stratification conditions occupied by a wide range of native and nonnative taxa, based on a previously developed habitat occupancy model for the regionally widespread invader P. inopinus . We also examine whether zooplankton communities exhibit sharp biogeographic breaks or slowly intergrade across broad distances, and whether we can identify geographic hotspots of introduction across the region.

Sample collection
The Pacific Northwest estuaries (PNWE) zooplankton time series comprises estuarine zooplankton community surveys conducted across the Pacific Northwest coast of North America. This program includes 38 estuaries and tributary rivers across more than 1000 km of coastline, with several recent sites extending into southern California (Fig. 1). We here present data from six surveys undertaken from 1996 to 2016. These surveys were conducted in collaboration with a variety of regional agencies and stakeholders every four years, with some interannual variation in the geographic extent of sampling (Table 1). This variation arises from sample design trade-offs between greater spatial coverage and reduced temporal resolution at the periphery of our study region, as well as the evolving research objectives of our regional collaborators.
In each survey year, we collected samples across an approximately 6-week window centered on the month of September. Zooplankton samples were collected from each water body along a longitudinal transect spanning a salinity gradient from 0 to 15 practical salinity units (PSU). Samples were collected from a small (4 m) boat at approximately mid-channel for each station via vertical net tows from just above the bottom to the surface using a 0.5 m diameter net. Zooplankton were collected using a 250-lm mesh net in 1996-2004, but a 75-lm mesh net was deployed from 2008 onwards to more effectively sample smaller taxa. Samples were preserved in the field in a 10% buffered formalin solution for later taxonomic processing.
In order to sample across a range of salinity conditions, target salinities at each station were set at 0, 1, 3, 5, and 15 PSU (measured at the bottom). In rare instances, we were unable to sample all target stations due to issues of channel depth or abruptness of the salinity gradient. Surveys of freshwater tributaries were conducted at 5 equally spaced stations along a~10-km transect. In all statistical analyses, we used measured salinity values rather than target salinity values.
Surface and bottom salinity and temperature, as well as water depth, were recorded at each station. Salinity and temperature values were measured using a YSI probe (Yellow Springs Instruments), and water depth was recorded using either an electronic depth finder or a marked and weighted line (depending on the year and boat). We calculated a simple stratification index as the difference between surface and bottom values at each station for salinity and temperature. This index serves as a rough indicator of mixing, but not an estimate of rate change across depth, or a measurement of the depth or shape of the thermocline or halocline. At the majority of our sites, however, the difference between surface and bottom salinity was <1 PSU, suggesting a generally well-mixed water column.
In the laboratory, zooplankton samples were washed through a 75-lm screen to remove formalin preservative, and then subsampled using a Hensen-Stempel pipette to obtain aliquots for taxonomic processing under a stereomicroscope. Plankton taxa were enumerated and adult individuals identified to genus, and in most cases, species. Each sample was processed until 100 specimens of the single most abundant taxon were enumerated. Larval forms were generally identified to higher taxonomic levels such as order (e.g., Calanoida), suborder (e.g., Balanomorpha), or class (e.g., Bivalvia). More than 200 zooplankton taxa were identified in total, which were then aggregated into 37 taxonomic groups to reduce the number of rare species across the dataset and to ensure taxonomic uniformity across sampling years (Table 2).

Statistical analysis
We conducted a multivariate statistical analysis of community data via a suite of nonparametric methods according to the methodology proposed by Field et al. (1982) and Clarke (1993). Zooplankton community structure was visualized via nonmetric multidimensional scaling (NMDS) of sample data (Kruskal 1964). Prior to ordination, abundance data were transformed to relative abundance within each sample (following Clarke 1993, McCune andGrace 2002), with taxa present in <3% of samples excluded from analysis. Samples were ordinated in a 2-dimensional space using a Bray-Curtis dissimilarity matrix, and the fit of the final configuration was assessed using the Dexter et al. (2018a) NMDS stress test, which evaluates ordination stress against a distribution of stress values generated from constrained permutations of the data. The associated Shepard diagram (Shepard 1980) was also visually inspected for outlying points of poor fit on the ordination space.
Individual samples were grouped according to similarity of community composition via hierarchical agglomerative clustering. Clustering was conducted using the flexible beta algorithm with beta at 0.6 (Milligan 1989) upon the Bray-Curtis dissimilarity matrix employed for NMDS ordination. A suite of summary statistics was calculated for each cluster, and clusters were visualized as an interpretive overlay on the NMDS ordination. Environmental factors (e.g., salinity, temperature, and stratification) were plotted as vectors on the NMDS ordination using the envfit function in the vegan package for R, which produces a goodness-of-fit statistic (r 2 ) for each environmental factor across the ordination space. The significance of each association was assessed using a Monte Carlo randomization test with 1000 permutations.
Finally, we characterized the salinity and stratification conditions under which the most abundant native and invasive taxa were observed. Pairwise comparisons were conducted for areas of overlap between all invasive taxa, as well as between the 10 most abundant native taxa and the regionally significant invader, P. inopinus. We delineated areas in this 2dimensional salinity-stratification space as the conditions under which at least one sample was at least 50% comprised of the given taxon. This high threshold was employed to exclude marginal habitats from the analysis. Unless otherwise specified, all species abundances were standardized to relative abundances (0-1 scale) within each sample. This standardization was performed to reduce the effects of variation in sampling gear (e.g., net mesh size) and personnel across years, and to more readily compare community composition across sites with markedly different levels of overall abundance (Field et al. 1982, Clarke 1993. All statistical analyses were conducted in R (version 3.5.1, R Core Team 2016) with figures generated using the colorblind-safe palette Viridis in ggplot (Wickham 2016). The complete R-code is provided in Data S1.

Overview
A total of 584 samples were collected from 38 estuaries and associated tributary rivers ( Table 1). Most of the estuaries were connected directly to the Pacific Ocean or the Salish Sea, but four tributary rivers draining into the Columbia River estuary were also surveyed (Grays, Cowlitz, Youngs, and Willamette rivers). Zooplankton taxa richness varied widely across samples, ranging from 1 to 18 taxa collected per station. The copepods Oithona davisae, Acartia spp., and Eurytemora affinis, and the larvae of bivalves and Cirripedia were the most abundant taxa overall (Fig. 2). For some of these taxa, high abundances were observed across a relatively small number of sites. When considered in terms of relative abundance across sites, the copepods Pseudodiaptomus inopinus, Pseudodiaptomus forbesi, Pseudobradya sp., and Eurytemora spp., and various Harpacticoida were the most abundant taxa (Fig. 3).
Bottom and surface salinity ranged from 0 to 34 PSU across the entire dataset, with sites ❖ www.esajournals.org ranging from unstratified (no difference between top and bottom salinity) to strongly stratified (34 PSU difference from top to bottom). Bottom temperatures ranged from 5.5°to 22.4°C and tended to be strongly correlated with surface temperatures, which ranged from 5.6°to 22.7°C. Across the entire dataset, differences between top and bottom temperatures ranged from 0°to 7.5°C.
Mixing index values of the sampled estuaries broadly reflected the coastal geomorphology of the study area (Emmett et al. 2000, Hickey and Banas 2003, Lee and Brown 2009. Sites located in northern Oregon and southern Washington (approximately 43°-47°N) tended to exhibit long and well-mixed zones of transition from fresh to marine waters, consistent with their origins as low-elevation drowned river valleys (e.g., the Tillamook River; Table 1). In contrast, sites at the northern and southern end of our study region tended to show stronger vertical stratification and short salt wedges at the mouth of the estuary (e.g., the Chetco River; Table 1). This pattern reflects the steep elevation gradient of the southern Oregon and northern California coast, and the glacial fjord origins of our northernmost estuaries (Emmett et al. 2000, Hickey andBanas 2003). These hydrological differences were apparent in the field, where sharp gradients of salinity prevented us from precisely sampling the target maximum salinity of~15 PSU at most northern and southern estuaries, but not estuaries in the central region (Fig. 4A). Hereafter, we refer to the central latitudinal range between 43°a nd 47°N (southern Washington and northern Oregon) as the "Central Coast" subregion, with sites located above 47°N (northern Washington and British Columbia) as the "North Coast" subregion, and with sites located below 43°N (southern Oregon and California) as the "South Coast" subregion.

Community-level patterns
NMDS ordination of the community data yielded an evenly distributed cloud of points without obvious breakpoints in community structure (Fig. 5). Higher-dimensional NMDS solutions (not shown) showed little change in the relative position of samples along the first two axes. Visual inspection of the associated Shepherd plot for outlying points (Field et al. 1982, Clarke 1993) did not show specific regions of poor fit on the ordination space. Results from the Dexter et al. (2018a) NMDS stress test (1-sample Z-test; P < 0.05) based on 1000 constrained permutations of the data matrix indicated that the stress value of 0.25 achieved by the 2-dimensional ordination remained within acceptable limits for interpretation. Hierarchical agglomerative clustering of community samples yielded three groups which could be distinguished by several biological and physical characteristics (Fig. 6). Division of the data into a greater number of clusters yielded groups that were difficult to distinguish from each other. Group 1 samples typically originated from higher salinity stations, with a mean salinity of 13.4 PSU compared with 4.0 and 2.6 PSU for groups 2 and 3, respectively. Group 1 encompassed samples from across the entire geographic range of the study, while groups 2 and 3 were predominantly found in Oregon and Washington. Samples from group 3 were associated with slightly warmer waters than groups 1 and 2, although all groups spanned a wide range of temperatures (Fig. 6). The three community groups are visualized on the NMDS ordination, with the associated environmental factors shown as vectors (Fig. 5). Note that samples are continuously distributed across ordination space and that samples positioned at group borders are likely intermediate in values.
Salinity, temperature, and stratification showed strong patterns of correlation across the NMDS ordination ( Fig. 5; Table 3). These associations are plotted as vectors on the NMDS ordination with vector length scaled to the r 2 value of the environmental variable.  divides groups 2 and 3) are associated with warmer water temperatures. In contrast, salinity and stratification both increase toward the lower right of this ordination space (i.e., the axis which divides group 1 from groups 2 and 3). Depth was uncorrelated with community structure in our analysis, likely because most samples were collected from small and shallow estuaries where depth could rapidly fluctuate with tidal stage and rainfall.
Community groups were also clearly differentiated based on their taxonomic composition. Although some taxa were common across all samples, each group could be readily characterized by its most abundant members (Table 4). Group 1 contained a high proportion of more marine-oriented taxa (e.g., Acartia spp., Cirripedia, polychaetes), reflecting the higher salinities typical of this group. Groups 2 and 3 were numerically dominated by brackish and freshwater copepods, but showed striking differences in species composition. Samples in group 3 were numerically dominated by the invasive copepod P. inopinus, comprising 83% of the average abundance. P. inopinus comprised only 4.3% of average abundance in group 2, with the harpacticoid copepod Pseudobradya sp. being the most abundant member of this group.

Nonindigenous zooplankton taxa
We observed a large increase in the number of nonindigenous zooplankton species across the 20-yr period of study, with most taxa found exclusively in the central coast subregion (Fig. 7). At the outset of the PNWE program, we observed only the invasive copepod P. inopinus among our surveyed sites. We then detected an increasing number of nonindigenous species at nearly every subsequent survey, such that by 2016, six nonindigenous species were detected across multiple sites: the copepods Pseudodiaptomus forbesi, Pseudodiaptomus inopinus, Oithona davisae, Sinocalanus doerrii, Limnoithona sinensis, and the cladoceran Bosmina coregoni (Fig. 8). With the exception of the copepod O. davisae, these species were found only in the central coast. This increase in the number of detected species does not appear to be driven by the 2008onwards deployment of a finer (75-lm) mesh net, as the small-bodied O. davisae were collected in large numbers prior to 2008, and our first detection of B. coregoni (which are large enough to be captured in both sizes of net) occurred 8 yr after the switch to a finer mesh.
The copepods P. forbesi, P. inopinus, and O. davisae have been previously characterized as invasive on the U.S. West Coast (Cordell and Morrison 1996, Cordell et al. 2008, and each comprised a major fraction of total zooplankton abundance across multiple estuaries that we sampled. P. inopinus was the single most abundant taxon in terms of relative abundance within individual samples (Fig. 3) and was observed in 200 samples across 14 estuaries. When present, P. inopinus typically comprised more than 50% of total zooplankton abundance, and in some cases more than 90% (Fig. 4B). We observed P. inopinus predominantly within the central coast subregion (Fig. 4B) although a small number of P. inopinus individuals were occasionally observed at more distant sites.
The congeneric copepod P. forbesi was detected in 12 samples collected across 3 sites: the Willamette River, the Youngs River, and the Grays River. Each of these rivers is a tributary of the lower Columbia River, where large populations of P. forbesi are established (Cordell et al. 2008, Fig. 6. Boxplots of selected biological and environmental variables across each of the three identified community groups.  Emerson et al. 2015). In all cases, P. forbesi were collected across a narrow range of salinities (0-7.3 PSU), with most specimens observed at <1 PSU. Although P. forbesi was detected at only a few locations, it was the most abundant taxon at several of the sampling stations where it occurred. The third invasive copepod, O. davisae, was detected in 56 samples collected across 10 different estuaries (the Russian, Samish, Eureka Slough, Nanaimo, Coos, Noyo, Elk, Skagit, Tijuana, and Fraser river estuaries). O. davisae were observed at temperatures ranging from 10°t o 22°C, and across a wide range of salinities, with roughly equivalent levels of abundance observed across 0.1-34 PSU. O. davisae were highly abundant at a small number of sites, but were more typically found in low numbers. O. davisae have been detected at an increasing number of sites in recent years and appear to be undergoing further spread across the region.
In contrast, the copepods S. doerrii and L. sinensis and the cladoceran B. coregoni were detected only in relatively low abundances. S. doerrii was detected in 8 samples across 3 different sites: the Chehalis River, the Youngs River, and the Russian River. It was detected across a relatively wide latitudinal range (38°-47°N) at salinities of 1-7 PSU. This represents a considerable expansion beyond the previously documented range of S. doerrii in North America (Sytsma et al. 2004, Cordell et al. 2007, Bollens et al. 2012), but in all cases, S. doerrii were observed in low abundances (<3% of total zooplankton).
Finally, the nonindigenous cladoceran B. coregoni was detected in five samples across four sites (the Chehalis, the Cowlitz, the Willamette, and the Tillamook rivers), two of which are tributaries of the Columbia River (the Cowlitz and the Willamette rivers). B. coregoni were primarily found at salinities of~0 PSU, and never in salinities greater than~3 PSU. In all cases, B. coregoni comprised only a small fraction of total zooplankton abundance (<3% of total zooplankton).

Environmental correlates and physicochemical niches of major invaders
Water bodies invaded by P. inopinus were characterized as having relatively low salinity (averaged across all stations) and little to no stratification (Fig. 9). As previous studies identified both of these factors as key predictors of P. inopinus spread ), we examined overlap in this salinity-stratification niche between P. inopinus and other native and invasive zooplankton. With respect to other dominant invasive zooplankters (P. forbesi and O. davisae), each exhibited a large degree of overlap with P. inopinus in terms of this salinity and stratification space, but not with each other (Fig. 10). In general, O. davisae were more likely to numerically dominate the zooplankton community in high salinity conditions, P. inopinus were more likely to dominate in brackish conditions, and P. forbesi Note: Taxa are ranked in descending order by the mean relative abundance of that taxon within each group. tended to dominate in freshwater conditions. With respect to native zooplankton, we found that nearly all examined taxa occurred under the salinity-stratification conditions in which P. inopinus typically dominated, although some occupied a much wider range of conditions (Fig. 11).
We identified several presently uninvaded estuaries which may, in terms of this salinity-stratification niche, be viable habitat for P. inopinus (Fig. 9). Most of these estuaries lie in the Puget Sound region of Washington state (Duwamish, Dungeness, Samish, Skagit, and Snohomish rivers) and have previously been characterized as potentially viable habitat for P. inopinus ). However, we found that the Russian River and the Eureka Slough  (both in central California) fall within this salinity-stratification niche as well.

Biogeographic patterns
Estuaries are highly dynamic environments in which physical and chemical conditions rapidly shift across time and space. Salinity and temperature are typically considered the primary environmental factors that determine estuarine zooplankton community structure (Collins and Williams 1982, Soetaert and Van Rijswijk 1993, Tackx et al. 2004, Marques et al. 2006, Graham and Bollens 2010, Bollens et al. 2011, Breckenridge et al. 2015, with water column stratification, eutrophication, and other anthropogenic disturbances exerting further influence (Laprise andDodson 1994, Marques et al. 2006). Our study confirms that salinity is a primary correlate of community differentiation within and across estuaries of the U.S. Pacific Northwest. We found that moderate-to-high salinity stations hundreds of kilometers apart tended to be far more similar to each other than to fresher stations situated slightly upstream in the same estuary and vice versa.
In contrast, temperature was a poor predictor of structure among the zooplankton communities that we studied, with most species being found across a wide range of temperatures. We observed that sites highly invaded by P. inopinus tended to be slightly warmer than uninvaded sites, but it is unclear whether this association reflects a causal relationship, or whether temperature merely serves as a proxy for other traits that differentiate invaded from non-invaded sites (e.g., flushing rates). A strong temporal correlation has been observed between water temperature and the seasonal abundance of the closely related P. forbesi in the Columbia River (Dexter et al. 2015), and it is interesting to see a similar (but weaker) correlation across space with P. inopinus.
Like salinity (but unlike temperature), stratification appears to be a primary factor differentiating the zooplankton communities that we surveyed. We found highly invaded estuaries to show low average salinity and little vertical stratification of the water column (i.e., well-mixed waters). These mixing regimes reflect local coastal geomorphology, with the low-elevation drowned river valleys of northern Oregon and southern Washington showing much higher rates of invasion than other regions of the coast. We also found that average salinity and stratification conditions across the entire estuary were a far better predictor of occupancy for P. inopinus than conditions at individual sampling stations. These Fig. 9. Mean salinity and degree of thermal stratification (averaged across stations) for each study site. Sites with similar conditions to invaded estuaries but which remain uninvaded (shown inside circle) are predominantly situated in the Puget Sound area. Fig. 10. Visualization of the salinity and stratification conditions under which the three highly abundant zooplankton invaders (Pseudodiaptomus inopinus, Pseudodiaptomus forbesi, and Oithona davisae) comprise more than 50% of total zooplankton abundance.
findings accord with studies of rocky intertidal communities along the U.S. West Coast, which showed that local coastal features tended to delineate biogeographic breaks, and that latitudinal gradients in temperature explained biogeographic patterns only at continental scales (Blanchette et al. 2008).

Accumulation of nonindigenous species
During our two-decade period of study, we observed a steady increase in the number of nonindigenous zooplankton among the estuaries that we surveyed. In 1996, we detected only a single invader, the copepod P. inopinus. During the next two decades, a new species was documented approximately every four years. We do not attribute this increase to the gradual inclusion of new sites across survey years, because new arrivals were primarily detected at consistently monitored sites in the central study region, not at peripheral sites.
This pattern of increase raises several questions. Should we expect this relatively constant increase in the number of nonindigenous species to continue in the coming years, or is the system nearing a saturation point? Why have new species arrived as a steady stream and not a single wave of introductions? What effect do these newly arrived species have on the establishment and persistence of other previously observed species? Are nonindigenous species facilitating the further introduction of new species (i.e., invasional meltdown; Simberloff and Holle 1999) or does competition limit establishment and coexistence among new arrivals?
Although our data do not allow us to directly examine ecological interactions, we observe several lines of evidence pointing toward strong competition among some invaders. For example, we can identify several instances in which a newly arrived invasive species appears to have displaced a previously established (and highly abundant) invader. The invasive copepod P. inopinus was one of the most abundant species of zooplankton in the Columbia River Estuary for at least a decade, but disappeared from the system shortly after the establishment of P. forbesi (Cordell et al. 1992, Sytsma et al. 2004. It is also noteworthy that despite the fact that the invasive copepods O. davisae, P. forbesi, and P. Fig. 11. Visualization of the salinity and stratification conditions under which Pseudodiaptomus inopinus comprises more than 50% of total zooplankton abundance (yellow polygons) relative to the same criteria for each of the 10 most abundant zooplankton taxa (green polygons), with areas of overlap shown in yellow-green.
inopinus dominate zooplankton communities under overlapping geographic ranges and environmental conditions, these three species rarely co-occur. These examples suggest that competition may play a strong role in constraining the establishment and spread of nonindigenous species in the Pacific Northwest. It remains unclear whether competition or hydrological conditions play the stronger role in constraining newly arrived species, or whether a combination of both factors ultimately delineates the extent of invasion spread.

Patterns of dispersal
Our surveys reveal several themes regarding the origin and dispersal of zooplankton invaders among Pacific Northwest estuaries. First, each of the six nonindigenous species that we detected is native to southeast Asian estuaries (Orsi et al. 1983, Orsi and Walter 1991, Uye and Sano 1998, Wang et al. 2007, Cordell 2012. Secondly, nearly all of these species have been collected from the ballast water of commercial cargo vessels entering Pacific Northwest waters (Cordell et al. 2009, Lawrence and. Finally, three of these species (P. forbesi, O. davisae, and S. doerrii) became established in the San Francisco Estuary prior to their appearance elsewhere in North America (Orsi et al. 1983, Ferrari and Orsi 1984, Orsi and Walter 1991. These patterns suggest an invasion corridor between southeast Asian estuaries and the Pacific coast of North America that is primarily driven by international shipping with secondary spread facilitated in part by domestic shipping-and strongly corroborates findings from ballast water surveys of ships entering Pacific Northwest harbors (Cordell et al. 2009, Lawrence and. While ballast water is the most probable vector for the initial introduction and secondary spread of nonindigenous zooplankton into larger estuaries, it is unlikely the cause of secondary spread to smaller estuaries across the region, because most of the estuaries harboring nonindigenous species have little or no commercial shipping (Pacific Maritime Association 2019). Despite the lack of a shipping vector, P. inopinus and O. davisae have rapidly colonized smaller rivers and estuaries across the region. O. davisae were absent from all surveyed sites in 1996, but occurred in at least 7 estuaries by 2000 (Fig. 8). We detected O. davisae at estuaries ranging from southern Canada to the United States-Mexico border. Similarly, P. inopinus spread across more than 500 km of coastline within approximately one decade (Cordell et al. 1992, Cordell andMorrison 1996).
Because zooplankton can disperse via varied natural and anthropogenic means (Havel and Shurin 2004), the absence of commercial shipping among smaller estuaries suggests that other transport vectors such as oyster aquaculture (Mckindsey et al. 2007, Minchin 2007, Dumbauld et al. 2009) or overland transport of recreational boats (Havel and Stelzleni-Schwent 2001) may have facilitated the spread of nonindigenous plankton species. Tides and coastal currents may also be effective means of dispersal for species that can tolerate marine waters (Christy and Stancyk 1982), while inland or salt-intolerant species may be transported on the plumage of aquatic birds (Frisch et al. 2007) or as wind-carried eggs (C aceres and Soluk 2002). In a previous study, a genetic-based reconstruction of the North American invasion of P. inopinus found that colonization and subsequent spread likely occurred in an erratic manner across the west coast, largely independent of commercial shipping routes (Dexter et al. 2018b).
The dispersal route of the nonindigenous cladoceran, B. coregoni, may be of particular interest in that (unlike other recently arrived zooplankton) B. coregoni are native to both Europe and Asia (M€ uller 1985, Demelo and Hebert 1994, Wang et al. 2007, Guijun et al. 2012. B. coregoni have been present in the North American Great Lakes since the 1960s (Deevey and Deevey 1971, Balcer et al. 1984, Lieder 1991, but it was not until the late 2000s that B. coregoni were detected at several sites across the Pacific Northwest (Smits et al. 2013, Dexter et al. 2015, Lee et al. 2016. Each of these sites is regularly monitored, and B. coregoni are easily distinguished from native Bosmina; thus, these detections likely represent a reasonable estimate for the time of population establishment. This raises the question of whether these populations originated via westward expansion from the Great Lakes, or from ballast water introductions from native populations in Asia. As zebra and quagga mussels continue westward expansion across the North American continent (Wong and Gerstenberger 2011), any information regarding dispersal of other taxa from the Great Lakes region to the Pacific Northwest would be invaluable in assessing the risk of these high-impact invaders being introduced to Pacific Northwest waters.

Mechanisms and impacts of P. inopinus invasion
Although multiple invasive species can dominate zooplankton communities across the Pacific Northwest, P. inopinus is unique in terms of numerical abundance and the sharp geographic boundaries defining its range. Based on a series of habitat occupancy models, Cordell et al. (2010) predicted that salinity and stratification would be the primary factors constraining the spread of P. inopinus and that a number of uninvaded estuaries in the Puget Sound region were at risk of invasion. Our results support both conclusions and furthermore identify several uninvaded estuaries on the California coast as viable habitat for P. inopinus.
Despite the apparent suitability of multiple uninvaded sites, the geographic range of P. inopinus has remained static for more than 20 yr. This cannot be attributed to the absence of transport vectors: P. inopinus have been detected in the ballast water of numerous ships entering the Puget Sound, which contains several estuaries apparently suitable for its colonization (Cordell et al. 2009). Puget Sound estuaries may possess physical or biological characteristics that prevent the establishment of P. inopinus-a question which merits further investigation. In contrast, we have recently detected a small number of P. inopinus at several of the California sites which we have identified as potentially viable habitat. We propose that P. inopinus may yet undergo further range expansion in northern California, with the Russian or the Rogue rivers serving as sources for southern expansion.
It is unclear whether highly abundant invaders such as P. inopinus have displaced native fauna, or whether they have simply occupied open niches in Pacific Northwest estuaries. Previous authors have noted that estuaries across the North American west coast tend to be relatively species-poor-hypothesized to be a consequence of the young geological age of these estuaries (Jacobs et al. 2004, Dumbauld et al. 2009). Several native taxa were consistently found in the brackish conditions where P. inopinus occur (Fig. 11), but few of those species are known to tolerate saline waters and their presence may be attributed to advection from upstream (fresher) reaches of the estuary. A notable exception is the brackish-tolerant copepod Eurytemora affinis, of which a unique subclade exists in the Pacific Northwest (Lee 2000), but which tended to occur only in low abundances across our survey sites. Now that we have a clearer understanding of the magnitude of zooplankton invasion across Pacific Northwest estuaries, we suggest that future studies more closely examine the ecological dynamics that underlie this set of invasionsboth in terms of the dynamics between native and invasive members of the community, and in respect to the dynamics between different invaders as well.

CONCLUSIONS
Across more than 1000 km of coastline, we found consistent patterns of zooplankton community differentiation within and across estuaries, primarily driven along axes of salinity, temperature, and stratification. Some estuaries in the Pacific Northwest region are highly invaded by multiple zooplankton species, while others appear resistant to new taxa. We observed six nonindigenous species of copepods and cladocera distributed among our sites, of which three (P. inopinus, P. forbesi, O. davisae) were sufficiently widespread and abundant to be considered invasive. In most cases, the geographic extent of all nonindigenous species extended far beyond previous literature reports. We observed a continual accumulation of nonindigenous species across our two-decade period of study. The three most widespread and abundant nonindigenous species rarely co-occurred and appeared to partition the environment across a gradient of salinity. The copepod P. inopinus was by far the most abundant and geographically widespread of these species, comprising more than 90% of all zooplankton abundance at some sites. The geographic range occupied by P. inopinus has remained stable for more than two decades, but some species (such as O. davisae) appear to be rapidly spreading across the region and can be locally abundant. We propose that the geographic distribution of these nonindigenous species may be largely determined by geomorphic characteristics that define the salinity and mixing regimes in estuaries-a reflection of the dominant role that these physical forces play in structuring estuarine zooplankton communities.

ACKNOWLEDGMENTS
The funding for this research was provided by a U.S. Environmental Protection Agency STAR grant (#FP91780901-0) awarded to E. Dexter and S. Bollens, and a series of grants from Washington Sea Grant/ National Oceanic and Atmospheric Administration, most recently #R/HCE/PD-5 to S. Bollens, J. Cordell, and G. Rollwagen-Bollens. We thank Ben Bolam, Sean Nolan, Vanessa Rose, and Josh Emerson for assistance in the field, and Olga Kalata for assistance with sample processing. Additional comments on this manuscript were provided by Stephanie Hampton, Stephen Katz, and S everine Vuilleumier as members of E. Dexter's Ph.D. committee.