Location Matters: Passive and Active Factors Affect the Vertical Distribution of Olympia Oyster (Ostrea lurida) Larvae

Dispersal, retention, and population connectivity are impacted by current regime and the behaviors that drive larval distribution, so understanding both is key to informing restoration of native species like the Olympia oyster (Ostrea lurida) across its range in western North America. This study explores the relationships between several factors (temperature, [chl a], larval size, tidal stage, and estimated current speed) and Olympia oyster larval vertical distributions in Fidalgo Bay (48.4828, − 122.5811), a shallow, tidally flushed bay in the Salish Sea. Olympia oyster larvae collected from four depths over the tidal cycle from July 11–14, 2017, were ~ 20% deeper near slack tide and shallower during the faster parts of both ebb and flood, with a threshold for this transition around an estimated 25 cm s−1. This pattern does not suggest tidally timed migrations as has been shown in another population of Olympia oysters, nor can this pattern be totally explained by passive processes. Larvae did not cluster at depths with specific temperatures or [chl a] but there was a difference in larval size between surface and bottom waters, with older, larger larvae more common at the bottom. Fidalgo Bay does not exhibit two-way flow or strong vertical shear, so vertical distribution of larvae likely has little effect on transport in this system but might in other similarly shallow habitat areas with higher stratification that are target restoration sites in the Salish Sea. These results add to the growing number of studies that show location-specific differences in larval vertical distribution and behavior within taxa and underscore the importance of integrating local hydrodynamics into predictions of bivalve larval transport.


Introduction
Larval behavior and the physical environment interact to affect larval dispersal, retention, and population connectivity. In particular, bivalve larvae are capable of actively controlling their depth through vertical swimming and sinking behaviors, thus determining which horizontal currents they are carried in and whether they are retained near their natal habitats (Hidu and Haskin 1978;Shanks and Brink 2005;Knights et al. 2006). While passive processes can explain vertical shifts in the distribution of bivalve larvae depending on the hydrodynamic system (Weinstock et al. 2018), biophysical transport modeling that incorporates both local hydrodynamics and larval behaviors suggests that behavioral modifications to vertical distribution can significantly affect transport and dispersal patterns (for example, with eastern oyster larvae in Chesapeake Bay; Dekshenieks et al. 1996;North et al. 2008). Thus, incorporating data on both local hydrodynamics and larval vertical distributions is key for predicting larval transport and dispersal potential in bivalves.
Understanding the relationship between larval behaviors and the physical environment of the Olympia oyster (Ostrea lurida) is critically important because this species is currently the focus of regional scale efforts to restore self-sustaining populations across their native range from British Columbia, Canada to Baja, México (Dinnel et al. 2009;Blake and Bradbury 2012;Wasson et al. 2016). Every summer Olympia oysters release brooded larvae, which develop Communicated by Judy Grassle Electronic supplementary material The online version of this article (https://doi.org/10.1007/s12237-020-00771-8) contains supplementary material, which is available to authorized users. further in the water column for one to several weeks (reviewed by Pritchard et al. 2015). Their depth distribution behaviors are likely driven by a variety of larval responses to water column conditions, including temperature extremes (Daigle and Metaxas 2011) and food patches (Raby et al. 1994;Sameoto and Metaxas 2008). Bivalve larvae are also known to undergo ontogenetic vertical migration during which the photonegative and geopositive late-stage larvae (Young 1995) distribute in deeper waters than do early-stage larvae (Carriker 1951;Mann 1988;Baker and Mann 2003;Poirier et al. 2019). Differences in behavior of early-and late-stage larvae can increase dispersal of newly released larvae and increase retention of late-stage, competent larvae near settlement habitat depending on local current regimes (reviewed by Forward and Tankersley 2001).
Larval behaviors, however, may be context dependent, and the physical differences between geographically distinct areas might influence the suite of environmental cues that drive the vertical distribution behaviors of larvae (López-Duarte et al. 2011;Miller and Morgan 2013a). For example, a single study on Olympia oyster larval distributions in Coos Bay estuary (Oregon) suggested that Olympia oysters undergo tidallytimed vertical migrations to remain in surface waters during flood tide and bottom waters during ebb tide (Peteiro and Shanks 2015). Tidally-timed migrations may enhance retention within Coos Bay estuary which is a relatively large (5 0 km 2 ) and deep coastal estuary connected directly to the outer Pacific coast and has prominent freshwater influence that creates weak halocline conditions to which Olympia oyster larvae may respond (Peteiro and Shanks 2015). In contrast, in Fidalgo Bay, WA-a relatively small (~6 km 2 ) and shallow bay in the central Salish Sea with little freshwater influence-our initial observations suggest that Olympia oyster larvae do not perform tidally timed vertical migrations, but instead distribute in surface waters during both ebb and flood (unpublished data from 2014 and 2015; Hatch et al. 2018). Moreover, larvae from geographically separate populations can show local behavioral adaptation (Manuel et al. 1996) and phenotypic plasticity (Miller and Morgan 2013b) that may further drive site-dependent behaviors.
The purpose of this study was to gain a better understanding of the relationship between Olympia oyster larvae and their behavioral responses to their environment to both inform restoration activities and increase understanding of bivalve larvae transport and dispersal. We focused our study in Fidalgo Bay, WA, because it is a Washington State priority restoration site for Olympia oysters (Blake and Bradbury 2012) and, due to collaborative restoration activities by state, tribal, and non-profit groups, it now has one of the largest Olympia oyster populations in the region at an estimated 4.8 million (Dinnel 2016). The large population made it possible to sample larvae in the water column during the reproductive season and observe larval distribution patterns in relation to other factors. As a shallow bay with little freshwater input (Sutherland et al. 2011;Banas et al. 2015), Fidalgo Bay is similar to other priority restoration areas (Blake and Bradbury 2012) and historical Olympia oyster bed habitats (Hatch and Wyllie-Echeverria 2016). Thus, results from this study in Fidalgo Bay can inform predictions about larval transport and dispersal in other parts of the Salish Sea. These predictions can aid restoration efforts by informing whether certain sites will be larval sources or sinks and where to focus habitat restoration activities to encourage recruitment for growing self-sustaining populations.
Because larval dispersal is a physical and biological process, we investigated the hydrodynamics of tidal currents and other factors that might influence vertical distribution of Olympia oyster larvae in the bay. Our two main goals were to (1) characterize current velocities that move through the prominent main channel to determine whether Fidalgo Bay exhibits vertical shear or a two-layer flow and (2) explore the effects of current speed, tidal direction, temperature, chlorophyll-a concentration ([chl a]), and larval size (a proxy for larval age) on the depth distribution patterns of Olympia oyster larvae in the main channel of the bay. Based on our preliminary observations that larvae were more abundant in surface waters during ebb and flood tide (unpublished data, Hatch et al. 2018), we hypothesized that larvae would not exhibit tidally-timed vertical migrations, but would instead exhibit ontogenetic vertical migrations with more abundant early-stage larvae in surface waters to enhance export and late-stage larvae in bottom waters to enhance recruitment and retention. These results can be used to predict larval dispersal to aid Olympia oyster restoration activities and, more broadly, apply to knowledge about bivalve larval interaction with local environments.

Study Location
Fidalgo Bay is located in the central Salish Sea adjacent to the city of Anacortes, Washington. It is relatively shallow, has little freshwater input from nonpoint source run-off, small creeks, and outfalls, and is open at its northern end to the tidally forced Guemes Channel (Fig. 1). The city of Anacortes converted a historic trestle into a walking path with riprap reinforcement, which narrows and channelizes the bay. Based on restoration surveys, the most densely populated Olympia oyster beds are on the east side of the bay along the main channel near the riprap and trestle (Dinnel 2016;Fig. 1). Olympia oysters might also inhabit the entire bottom of the main channel, but this has not been surveyed. Adjacent to the main channel and south of the trestle are extensive mudflats that become mostly exposed during low tides of the bay's semi-diurnal tidal cycle. Within the main channel, we sampled larvae by boat at a location that was at~-3.5-m tidal elevation (labeled "Boat" on Fig. 1) and collected current velocity profiles using an acoustic Doppler current profiler at a location that was at~-4-m tidal elevation (labeled "ADCP" on

Experimental Design
We designed a sampling schedule that enabled us to use mixed effects modeling to investigate the effects of both depth-specific and depth-averaged water column conditions (estimated current velocity, temperature, and [chl a]) on larval depth distribution patterns and the potential for ontogenetic vertical migration. We simultaneously collected larvae, temperature data, and chorophyll-a samples during flood and ebb tides on 4 days from July 11-14, 2017, via the field sampling methods described below. Current velocities were measured during a separate 3-day (July 25-28, 2017) deployment of a bottom-mounted acoustic Doppler current profiler (ADCP) instrument and then used to estimate currents during the week of larval sampling.

Field Sampling
We measured larval abundance, temperature, salinity, and [chl a] from four depths at one location by boat each day from July 11-14, 2017 ( Fig. 1) to match the estimated peak reproductive timeframe for Olympia oysters (Dinnel 2016). We selected 4 days based on NOAA (2017) tidal predictions with enough tidal exchange to sample during both ebb and flood tide during daylight hours because Olympia oyster larvae are phototactic (personal observations). Each day, we completed 11 sampling events during which we collected samples from four depths in the water column: surface (0.5-m below surface), bottom (0.5-m above seafloor), and two mid-depth samples that evenly split the depth between surface and bottom samples. Before each sampling event, we measured the water depth with a sounding line to calculate the four sampling depths. It took 15-20 min to collect all four depth-specific samples during each sampling event. By the end of the four consecutive sampling days, we had performed 44 sampling events and collected a total of 176 individual samples.
To collect each larval sample, we used a modified bilge pump to filter 100 L of water from each targeted depth through a 102-μm mesh plankton net to ensure retention of Olympia oyster larvae, which are released from brooding at a shell length (SL) of~180-μm (reviewed by Pritchard et al. 2015). Each sample was stored on ice while in the field and then preserved in 70% ethanol until processing during the following 8 weeks; there was no visible shell degradation. At the end of filtering each 100-L sample, we collected 60 mL of bulk seawater from the pump for measurement of [chl a]. We filtered the 60 mL of seawater through a 0.7-μm pore size glass microfiber filter (Whatman™ GF/F); the foil-wrapped filters were held on ice in the field and then stored at − 80°C for later extraction. We extracted chl a using 90% acetone for 24 h and read fluorescence of each sample with a Turner Trilogy Fluorometer. We collected temperature and salinity measurements at the same times and depths as our larval pump sampling with a Hach Environmental Company HydroLab DS5 water quality multiprobe instrument.
We measured current velocities in the main channel of Fidalgo Bay with a Nortek 1-MHz Aquadopp ADCP fixed to a steel bottom frame (Fig. 1). Velocity was recorded in 0.3-m vertical bins every 60 s. The ADCP and frame were set center-channel (4-m depth with an estimated mean tidal range of 1.5 m and diurnal tidal range of 2.5 m). We carefully chose a location along the same path of water flow in the main channel as the larval sampling site (97-m apart), but with firmer sediment to provide stability for the instrument and reduce the likelihood of obstruction of the sensors by fine sediments. These efforts were successful; pitch (− 0.03°± 0.05) and roll (− 1.7°± 0.08) remained small and stable during the deployment.
Velocity measurements with the ADCP did not coincide with our larval sampling period due to instrument failure. Therefore, we chose a second timeframe from July 25-28, 2017, to redeploy the instrument that best represented the larval sampling period based on similar tide schedule (Fig. 2b). Current velocity data collected during the second deployment were used to estimate current velocities during the larval sampling period of July 11-14. These estimations work for this narrow and highly constrained system because sea surface pressure is a strongly dominant force of current velocities here (Fig. 2a). At this site, the majority of the signal of depth-averaged currents can be predicted from tidal pressure; therefore, once currents are determined accurately for one tidal cycle, they can be estimated for another tidal cycle with similar amplitudes (Fig. 2b). We compared depth-averaged alongisobath current velocities (1-to 4-m above seafloor) measured in Fidalgo Bay to sea surface pressure measurements in the adjacent Guemes channel (48.5092, − 122.6871) collected 90 min later from a Seabird SBE-39 and found a significant correlation ( Fig. 2a; Pearson's cor = 0.84, p << 0.001, adjusted R 2 = 0.70, p << 0.00), further suggesting that the current system within Fidalgo Bay is predominantly driven by sea surface pressure. Additionally, local weather conditions during the two measuring periods were generally similar ( Fig. 3; The Weather Company 2019); wind-driven pressure gradients did not cause currents to differ substantially between sampling timeframes. Thus, we estimated current velocities during the larval sampling timeframe with measured current velocities at corresponding tidal stage and sea surface pressures. The~0.5m difference in depth between the ADCP sampling location and the discrete water sampling location was also taken into account in adjusting for depth-specific and depth-averaged estimates. Finally, current vectors shown in Fig. 4 were rotated from directional earth coordinates (east-north-up) into a depth-averaged axis of principal variance (coordinates: along-channel, across-channel, up), calculated using empirical orthogonal function analysis.
We hand-sorted Olympia oyster larvae from each sample using an Olympus SZ-ST stereoscope fit with polarized lens filters. Olympia oyster larvae were distinguished from the other bivalves in our sample using several methods. First, we narrowed down all the potential local species of bivalve larvae based on reproductive season (Loosanoff et al. 1966;pers. comm. Julie Barber). We distinguished Olympia oyster larvae from the larvae of common bivalves in Fidalgo Bay (pers. comm. Julie Barber) by size and shell shape using reference Olympia oyster larvae and morphological identification keys (Loosanoff et al. 1966;Shanks 2001). The most abundant bivalve larvae in the samples were clams, which are easily distinguishable from oyster larvae by shell shape. The Pacific oyster (Crassostrea gigas) is the only other oyster in the region and we could easily distinguish the small proportion of oyster larvae that were C. gigas due to their skewed umbo (Loosanoff et al. 1966;Shanks 2001). Once we had sorted all Olympia oyster larvae, we measured shell lengths parallel to the hinge of each larva using a stereomicroscope equipped with a camera and the ImageJ software (Leica MC170 HD and Leica Application Suite, Leica, Wetzlar, Germany). Olympia oyster larvae were then divided into three size classes: 180-210 μm were newly released larvae, 211-260 μm were developing larvae, and > 260 μm were latestage larvae close to settlement. We chose these size classes based on reference materials (Loosanoff et al. 1966;Shanks 2001) and personal observations of development of lab-reared Olympia oyster larvae from the Fidalgo Bay population.

Data Analysis
We analyzed our data with mixed effects models to investigate the influence of depth-specific and depth-averaged environmental parameters on larval depth distribution patterns and whether Olympia oyster larvae exhibited ontogenetic vertical migration. Mixed effects models enabled us to account for the inherent violations of independence in our sampling design (i.e., we collected all samples from one fixed location over time) and the hierarchical structure of our sampling design (i.e., we collected four depth-specific samples per sampling event and 11 sampling events per sampling day) (Zuur et al. 2009). By accounting for these sources of error within our models, we were able to include all 4 days of depth-specific samples in each analysis. We represented the depth distribution of larvae as a model response variable in two distinct ways: (1) depth-specific larval abundance per sample (100 L) and (2) larval weighted mean depth (WMD) normalized to water column height for each sampling event, which was calculated from the four depth-specific samples as follows: is the proportion of the total number of Olympia oyster larvae in the ith depth-specific sample. S i is the sampling depth (m) of the ith depth-specific sample. H is the water column height (m) at the time of sampling.
Each larval WMD measurement was a ratio from 0 to 1; low values represent relatively shallower larval WMD and high values represent relatively deeper larval WMD.
For each model, we began by including all predictor variables of interest, then determined the most appropriate random component if needed, and then selected the most parsimonious final model by comparing nested models with log likelihood tests, the Akaike information criterion (AIC), and Bayesian information criterion (BIC). We calculated the conditional and marginal R 2 values for each final model, which indicate variance explained by the entire model and by the fixed factors alone, respectively ("MuMIn" package in R version 3.5.0, Nakagawa and Schielzeth 2013; Bartoń 2018). We validated the fit of each final model by visually inspecting residual and QQ plots. Additionally, we inspected data for temporal autocorrelation structures using autocorrelation function plots.
First, we fit a linear mixed effect model (LME) to investigate the fixed effects of depth (m), estimated current velocity (m s −1 ), temperature (°C), and [chl a] (μg L −1 ) on the abundance of larvae per 100-L sample ("nlme" package, Pinheiro et al. 2018). We included a random intercept of "sampling event" nested within "sampling day" to account for this inherent non-independence of our sampling design. We built the Fig. 2 Acoustic Doppler current profiler depth-averaged (1-to 4-m above seafloor) along-isobath current velocities (cm s −1 ) measured July 25-28 in Fidalgo Bay (48.4828, − 122.5811) plotted as a function of sea surface pressure (Decibars) measured 90-min later in Guemes Channel (48.5092, − 122.6871) with a Seabird SBE-39 (a). Sea surface pressure data were collected in 30-s intervals, so were averaged per min to match the 60-s measurement interval of the current velocity data. A 10-min moving average was applied to each dataset before correlation was analyzed and plotted. Linear relationship (Pearson's cor = 0.84, p << 0.001; adjusted R 2 = 0.70, p << 0.001) is shown by the gray dashed line. Sea surface pressure (decibars) during larval sampling (July 11-14; black) and current sampling (July 25-28; gray) are shown in (b) model to allow for unique variance structures by water column height to correct for an observed violation of variance homogeneity.
Second, we fit a generalized additive mixed model (GAMM) and a linear mixed effect model (LME). Preliminary plots suggested a non-linear relationship between larval WMD and estimated depth-averaged current velocity, so we fit a GAMM with a thin plate regression spline for the estimated depth-averaged velocity parameter to investigate the significance of this non-linear relationship ("mgcv" package; Wood 2003Wood , 2011Wood , 2017. We then fit an LME to investigate the fixed effects of estimated depth-averaged absolute current speed (cm s −1 ), current direction (ebb or flood), depthaveraged [chl a] (μg L −1 ), and proportion of late-stage larvae (> 260 μm) on larval WMD ("nlme" package; Pinheiro et al. 2018). We included a random intercept for "sampling event" to account for this violation of independence in the sampling design. A log likelihood ratio test using residual maximum likelihood estimation indicated that including sampling day as an additional part of the random component did not improve the model, so we removed this term to reduce model complexity. We performed a simple contrast from the full LME model to compare larval WMD on ebb versus flood tide after ensuring equal variance using Levene's test ("car" and "gmodels" packages; Fox and Weisberg 2011;Warnes et al. 2018).
Finally, to investigate whether the larvae appeared to undergo ontogenetic vertical migration, we first compared the WMD of the three size classes of larvae (180-210 μm, 210-260 μm, and > 260 μm) with a one-way ANOVA after applying Levene's test. Then, we applied a Kolmogorov-Smirnov test (KS test) to determine if the size distributions of Olympia oyster larvae were significantly different in the surface versus the bottom samples.

Results
Dominant tidal currents, along-isobath, flowed along the main channel with negligible cross-channel flow, across-isobath, and little variation with depth (Fig. 4a, b). Positive alongisobath currents (Fig. 4a) flowed shoreward toward the trestle and inner bay, which is in a southeast direction, and negative along-isobath currents (Fig. 4a) flowed seaward, which is in Fig. 3 Pressure (mbars), wind (m s −1 ), and air temperature (°C) plotted in 5-min intervals from July 13 at 00:00 to July 14 at 23:57 shown by black dots and from July 25 at 00:00 to July 28 at 23:57 shown by gray dots. Data are from a weather station on the south eastern end of Guemes Island, Anacortes, WA (Fig. 1,  48.533, − 122.585; The Weather Company 2019) the northwest direction (Fig. 1). Positive, across-isobath currents flowed southwest and negative across-isobath currents flowed northeast, and there was negligible across-isobath flow in either direction (Fig. 4b). Along-isobath currents reached 65 cm s −1 during maximum ebb and flood tides and were consistently < 50 cm s −1 during low tidal exchanges, which corresponded with the tides during each larval sampling period (Fig. 4a). There was little variation in along-and acrossisobath velocities from the surface to the seafloor at each sampling timepoint (Fig. 4a, b), although along-isobath surface water (just~0.4-m depth) appeared to maintain a seaward flow (negative current speeds) during some slack and flood tides (Fig. 4a). In addition, WSW wind influenced the current velocity of surface water to about 0.4-m depth on July 27 (Fig.  4b), which corresponded to the relatively higher wind speeds that day (1.8 to 6.6-m s −1 , 4.2 ± 1.6-m s −1 from 6:00-21:00) compared to the rest of that week (0-4.5 m s −1 , 1.1 ± 0.7 m s −1 ; x ± SE; Fig. 3) and the week of larval sampling (0-3.8-m s −1 , 1.2 ± 0.7 m s −1 ; x ± SE; Fig. 3). Vertical velocities remained slightly positive throughout the deployment with frequent fluctuations of upward and downward movement (1.6 ± 0.99 cm s −1 ; x ± SE; Fig. 4c). A 5-min moving average was applied to each dataset as displayed. In a and b, positive and negative current speeds correspond with flood and ebb tides, respectively. In c, positive and negative values correspond with upward and downward water movement, respectively. Black rectangles in a frame low tidal exchanges that were used to estimate current velocities during larval sampling from July 11-14. Black arrows in c point to periods of the low slack tide where we observed larvae deeper in the water column. See Supplementary Fig. 6 for grayscale plot Overall, we collected five times fewer larvae (mean 0.58 ± 0.48 larvae L −1 , range 0.02-3.17 larvae L −1 ) than we have in previous years in Fidalgo Bay (0.2-15.7 larvae L −1 ; Hatch et al. 2018). We used an LME model to determine whether depth (m), estimated depth-specific current velocity (m s −1 ), depth-specific temperature (°C), and depth-specific [chl a] (μg L −1 ) influenced larval abundance. Of these factors, only depth significantly affected larval abundance (LME t 131 = − 8.07, p << 0.001, marginal R 2 = 24.8; Tables 1 and 2; see Supplementary Fig. 1 for validation plots). Depth profile plots support these model results: depth distribution of larvae changed over the tidal cycle each day but did not appear to be related to depth-specific temperature or [chl a] (Fig. 5, see Supplementary Figs. 2-4 for depth profile plots of each sampling day). During field sampling, temperatures ranged from about 12 to 19°C (16.6 ± 1.3°C; x ± SE) with little temperature stratification each day. The maximum variation in temperature from surface to bottom was 1.1°C during ebb tide and 4.1°C during a few flood tide sampling events. The [chl a] ranged from 4.4 to 38.6 μg L −1 (19.1 ± 8.1 μg L −1 ; x ± SE), but there were no evident chl a maxima (Fig. 5, Supplementary Figs. 2-4).
Next, we investigated effects of estimated depth-averaged absolute current speed (m s −1 ), tidal direction (ebb vs. flood), and depth-averaged [chl a] (μg L −1 ) on larval WMD visually and using a GAMM and an LME. Depth profile plots show that larval vertical distribution patterns throughout the tidal cycle were similar all four sampling days ( Fig. 5; Supplementary Figs. 2-4). Specifically, larvae were shallower during ebb tide, deeper as the currents slowed around slack tide, and shallower when currents increased during flood tide ( Fig. 5; Supplementary Figs. 2-4). Our modeling results supported this visual pattern, indicating that larval WMD and estimated depth-averaged current velocity had a significant non-linear relationship (GAMM, F edf of 4 = 11.7, p << 0.001; Fig. 6a; See Supplementary Fig. 5 for corresponding model validation plots). More specifically, the model indicated that estimated absolute depth-averaged current speed, not current direction (ebb or flood), significantly affected larval WMD (LME; Tables 3 and 4, Fig. 6b, c). Larval WMD appeared to be~20% shallower during ebb and flood tides and deeper around slack tide, with a threshold for this transition2 5 cm s −1 (Fig. 6a). In addition, larval WMD was significantly shallower when depth-averaged [chl a] (μg L −1 ) in the water column was higher (LME, Tables 3 and 4, Fig. 7).
We also examined whether larvae were undergoing ontogenetic vertical migration by including the proportion of latestage larvae (> 260 μm) in the LME model, comparing the depth distribution of larval size classes, and comparing larval size distributions at the surface and bottom. Results showed that every size class of larvae exhibited a similar WMD pattern in relation to estimated current velocity (Fig. 6a), and this observation was confirmed by the LME model selection process indicating that proportion of late-stage larvae (SL > 260 μm) did not significantly affect larval WMD (Tables 3  and 4). Of the total Olympia oyster larvae sampled (n = 10,124), 48% were 180-210 μm, 43% were 210-260 μm, and 7% were > 260 μm, but there was no significant variation  in WMD among these three size classes of larvae (one-way ANOVA, F 2,129 = 1.81, p = 0.17). A η 2 of 0.03 suggested that size class to size class differences accounted for just 3% of the total variability in larval WMD. However, the size distribution of Olympia oyster larvae in the surface and bottom samples were significantly different (two-sample K-S test; D = 0.14, p value < 0.001). The larval size distribution was skewed, with a higher frequency of smaller larvae in the surface samples than the bottom samples (Fig. 8).

Discussion
Olympia oyster larvae in the Fidalgo Bay main channel were consistently more abundant near the surface during ebb, deeper during slack, then more abundant near the surface again during flood. We observed similar patterns of higher abundances of Olympia oyster larvae in surface than in bottom waters during ebb and flood in Fidalgo Bay in 2014 (unpublished data) and 2015 (Hatch et al. 2018). Our results do not suggest selective tidal stream transport (Forward and Tankersley 2001) yet cannot be totally explained by passive processes either (Weinstock et al. 2018). The pattern of Olympia oyster larval depth distribution over the tidal cycle was best explained by estimated tidal current speed (LME marginal R 2 = 42.2%; Table 4, Fig. 6b), and we hypothesize it is driven by a combination of passive and active larval movement.
Other studies, including one study on Olympia oyster larvae (Peteiro and Shanks 2015), have inferred that bivalve larvae are capable of active vertical distribution based on observed shifts in larval abundance with depth over the tidal (μg L −1 ), and white diamonds represent estimated current velocity (cm s −1 ). Negative and positive values of estimated current velocity indicate water moving out of the bay (ebb tide) and into the bay (flood tide), respectively. Horizontal dotted lines represent sea surface adjusted so the bottom of each plot represents the seafloor Fig. 6 Larval WMD as a function of estimated depth-averaged current velocity (cm s −1 ) fit with generalized additive modeling thin plate regression spline smoother for all larval size classes combined and shading representing 95% confidence intervals (a). Larval WMD versus (b) estimated absolute tidal current speed (cm s −1 ) fitted by linear mixed effect modeling and (c) ebb and flood tidal direction. Linear mixed effects modeling indicated a significant linear relationship between larval WMD and estimated current speed (LME, T = −4.87, p << 0.001) and no significant effect of tidal direction (LME). A simple pairwise contrast confirmed no significant difference between ebb and flood larval WMD (t 1 = 0.80, p = 0.43) Table 3 Model selection process using log likelihood ratio tests and comparing Akaike information criterion (AIC) and Bayesian information criterion (BIC) to determine which fixed effects explain variance for larval WMD including data from all sampling days. Fixed effects include average [chl a], estimated depth-averaged absolute current speed (CurrentSpeed), current direction (Tide), and proportion of late-stage larvae (PropLate). Also shown are degrees of freedom (DF) and log likelihood values (LogLik). The p value (p) is likelihood of calculating the log likelihood ratio (L.Ratio) if the null hypothesis (full model has a better fit than reduced model) is true cycle and have suggested that these behavioral changes lead to enhanced transport during flood tides and retention during ebb tides (e.g., Carriker 1951;Mann et al.1988;Forward and Tankersley 2001;Peteiro and Shanks 2015). However, observed tidal shifts in larval depth might also be explained by passive processes. Weinstock et al. (2018) showed via metaanalysis that local tidal range and slope, proxies for tidal vertical advection, explain a significant amount of variation in results among published studies that indicate a tidal shift in bivalve larval vertical distribution. Further, they found the magnitude of estimated tidal vertical advection corresponded with the degree of change in larval vertical position from ebb to flood tide (Weinstock et al. 2018). Mean variations in vertical velocity in Fidalgo Bay were slightly positive (1.6 ± 0.99 cm s −1 ) and exceeded our lab-measured maximum swimming and sinking speeds of Olympia oyster larvae of 0.4 and 0.8 cm s −1 , respectively (unpublished data). Nevertheless, there was no indication of downward movement of water during slack tide when our observed larval WMD deepened, so passive vertical advection of larvae does not fully explain the patterns in larval distribution that we observed. Moreover, Knights et al. (2006) similarly observed the highest density of mussel larvae nearest the seafloor when tidal currents were slowest at slack tides. In a follow-up modeling study, James et al. (2019) "reverse engineered" this pattern and showed that, in the Irish Sea, the patterns of shifts in bivalve larval distributions observed by Knights et al. (2006) could only be explained by active swimming behavior changes performed by larvae throughout the tidal cycle. Even if passive processes are at play, we should not discount active movement entirely; Olympia oyster larvae actively swim upward in the laboratory (unpublished observations) Fig. 7 Relationship between larval WMD and depth-averaged [chl a] (μg L −1 ) (LME, T = − 3.8, p << 0.001, Table 3). Thick black line shows average model fit of fixed depth-averaged [chl a] modeling component of the LME (i.e., relationship between larval WMD and depth-averaged [chl a] for all sampling events across four sampling days). Grayscale lines represent the model fit of this relationship for each sampling event and some biophysical models that incorporate local hydrodynamics reveal that eastern oyster larvae (Crassostrea virginica) in Chesapeake Bay, Mobile Bay, Eastern Mississippi Sound, and Delaware Bay do not distribute like passive particles and can affect dispersal patterns through active vertical movement (North et al. 2008;Kim et al. 2010;Narváez et al. 2012). The vertical distribution of Olympia oyster larvae in Fidalgo Bay might be driven by a combination of passive and active larval movement.
Our study found the larval WMD was significantly shallower when there was higher average [chl a] in the water column (Fig. 7), but higher depth-specific larval abundances in the surface during ebb and flood tide were not related to temperature or [chl a]. We suggest that [chl a] may have been a weak predictor of overall larval WMD (estimated R 2~1 3%; Table 4) due to light attenuation from a high density of phytoplankton during sampling. Like other bivalves (Bayne 1964;Barile et al. 1994), Olympia oysters are highly phototactic (unpublished observations) and so might distribute shallower when higher abundances of phytoplankton block light from penetrating deeply. We did not measure light intensity directly, but consistently shallow Secchi depths (0.75-1 m) suggest that even in the shallow waters of Fidalgo Bay light might be an important predictor of larval depth distribution. However, larvae did not remain near the surface during the entire tidal cycle as we would expect if they were actively and persistently swimming upward due to a phototactic response, so light does not fully explain the observed larval vertical distribution pattern.
Oyster larvae can also behaviorally respond to hydrographic conditions, which might have enabled Olympia oyster larvae in Fidalgo Bay to influence their depth distribution during our sampling period. In flume and grid-stirred flow tank experiments, eyed eastern oyster larvae (Crassostrea virginica) actively respond to hydrographic cues related to turbulence by increasing upward swimming speeds or by rapidly diving (Finelli and Wethey 2003;Fuchs et al. 2013;Fuchs et al. 2015;Wheeler et al. 2015). In a racetrack flume study, Finelli and Wethey (2003) estimated that when larvae dive at 0.8 cm s −1 , which is also the maximum downward velocity we measured for late-stage Olympia oyster larvae in the lab (unpublished data), they exert enough propulsive force to potentially control their depth in horizontal current speeds up to 17-52 cm s −1 . This estimation was based on two assumptions: (1) that larvae can control their vertical position when Rouse numbers, which are ratios of sinking velocity over shear velocity multiplied by von Karmen's constant, are > 0.75 (Gross et al. 1992) and (2) that shear velocity is 0.05-0.15 times freestream velocity over smooth bottoms (Finelli and Wethey 2003). Recognizing that Fidalgo Bay is not a controlled bottom boundary layer as in the laboratory experiment, we can cautiously use the estimations by Finelli and Wethey (2003) to hypothesize that Olympia oyster larvae in Fidalgo Bay could have been actively responding to hydrographic conditions and influencing their depth because horizontal current speeds were < 50 cm s −1 during our larval sampling. This hypothesis is supported by Peteiro and Shanks (2015) who observed Olympia oyster larvae in the Coos Bay estuary were significantly deeper during ebb than flood tides when currents speeds were less than 50 cm s −1 . In our study, larvae were deeper near slack tide and shallower during the faster parts of both ebb and flood (Figs. 5 and 6a), with a threshold for this transition at~25 cm s −1 . If larvae in Fidalgo Bay were actively influencing their depth, our results suggest they might have been actively sinking during slower current speeds or behaviorally responding to hydrographic cues related to turbulence by swimming upward faster during higher current speeds, similar to the behavior of eyed eastern oyster larvae in lab experiments (Wheeler et al. 2013;Fuchs et al. 2015). Interestingly, previous studies on upward swimming and diving behaviors in oyster larvae focused on only larvae that were competent or nearly competent to settle (Finelli and Wethey 2003;Wheeler et al. 2013;Fuchs et al. 2015), but we observed a consistent change in vertical distribution depending on tidal currents for all sizes of Olympia oyster larvae (Fig. 6a). Fuchs et al. (2015) noted that some early-stage snail larvae behaviorally respond to turbulence and wave accelerations, but the implications of behavioral responses to flow by early larval stages have not been thoroughly explored (but see DiBacco 2011 andWheeler et al. 2016).
Although Olympia oyster larvae in Fidalgo Bay do not exhibit tidal vertical migration, they might undergo ontogenetic changes in vertical distribution. While modeling suggested that all sizes of larvae exhibited a similar vertical distribution pattern (GAMM; Fig. 6a) and did not identify larval  Tables 3 and  4), we might have failed to detect an effect of larval size because of the dominant influence of predicted tidal current speed on larval WMD and the relatively low abundance of late-stage larvae (~7%) in our samples. The size distribution of larvae in the surface skewed significantly smaller than in the bottom (K-S test; D = 0.14, p value < 0.001; Fig. 8), which suggests that larger, late-stage larvae distribute deeper due to passive sinking related to buoyancy or that larvae undergo ontogenetic vertical migrations that are detectable even in this relatively shallow bay.
The distribution patterns of Olympia oyster larvae in Fidalgo Bay differed from those observed in Olympia oyster populations in Coos Bay, OR (Peteiro and Shanks 2015). The two populations are genetically distinct (Silliman 2019), so might have evolved different behavioral strategies to optimize larval survival. While larvae within Coos Bay estuary must be retained to avoid wastage into the open coast (Peteiro and Shanks 2015), larvae in the Salish Sea have opportunities to encounter potential habitat even if they are transported far from their natal population, as suggested by widespread historic distribution of Olympia oyster beds (Blake and Bradbury 2012; Hatch and Wyllie-Echeverria 2016) and population genetic studies that indicate little genetic differentiation (i.e., high larval exchange) at the scale of Salish Sea sub-basins (Stick 2011;Silliman 2019). Nevertheless, despite the fact that Olympia oyster larvae in Fidalgo Bay accumulate near the surface during both ebb and flood tide, larval vertical distribution patterns probably have little effect on transport of Fidalgo Bay populations because the bay does not exhibit two-way flow or strong vertical shear due to its lack of freshwater input and high channelization. However, we might expect that vertical distribution and ontogenetic migration will have a stronger influence on Olympia oyster larval transport and dispersal from other restoration sites and potential habitat areas that are similarly shallow bays but exhibit sheared flow or stratification due to geomorphology and freshwater input (Blake and Bradbury 2012;Hatch and Wyllie-Echeverria 2016). These areas might support source populations for the species, as larvae that distribute in the surface during ebb and flood tide in these habitats will likely have a net seaward transport and disperse from their natal populations. In contrast, larvae are probably mostly retained and settle within Fidalgo Bay based on results from this study and observed settlement patterns (Dinnel 2016) and so is likely a sink population for the central Salish Sea.
Differences in larval distributions between populations have been shown for several taxa and might be due to differing physical conditions that influence passive movement and/or trigger differing larval behavioral responses (Morgan et al. 2014) and could be driven by phenotypic plasticity or genetic differences to enhance larval survivorship given local conditions. Phenotypic plasticity in larval behavior has also been observed in several species and might be common for species with large ranges occupying diverse habitats (Morgan and Anastasia 2008). For example, certain fiddler crab species match their tidal vertical migrations with local tidal rhythms to enhance larval export (Morgan and Anastasia 2008), and shore crab larvae of adjacent coastal and estuarine populations exhibit different diel and vertical migration behaviors to enhance larval transport to offshore nursery areas given different spawning locations (Miller and Morgan 2013b). Giant sea scallop larvae from distinct populations display significantly different vertical migration patterns in the same laboratory conditions (Manuel et al. 1996). If the contrasting Olympia oyster larval distributions in Fidalgo Bay and Coos Bay were driven by differing physical conditions, depth is the most likely contributing factor. Sampling sites in both bays exhibited channelized tidal currents with comparable speeds, but Fidalgo Bay is about three times shallower than Coos Bay (3.5-m average depth during sampling compared to 12 m), which might have altered passive movement from differing slopes and tidal vertical advection (Weinstock et al. 2018) or active movement based on bivalve larval response to pressure (Bayne 1964;Mann and Wolf 1983;Young 1995).
Results of this study underscore that location matters when considering larval distribution patterns; local physical characteristics can drive both passive and active larval movement so neither can be discounted when predicting larval transport and dispersal potential. Nevertheless, while this study broadly highlights the importance of understanding local environmental and hydrodynamic conditions to best predict larval transport patterns, results from this study can also give us a first indication of the key factors that impact vertical distribution of Olympia oyster larvae and can inform larval transport and dispersal patterns from common habitats in the Salish Sea.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.