Tidal shifts in the vertical distribution of bivalve larvae: Vertical advection vs. active behavior

Some previous studies, though not all, have reported that vertical distributions of bivalve larvae shift upward in the water column on flood tides and downward on ebb tides. When observed, such shifts have been interpreted as reflecting the vertical migration of larvae (i.e., an active behavioral process). This interpretation assumes that tidally driven vertical advection is insignificant compared to larval swimming speeds, but that assumption may not be valid in all regions. We demonstrated that the vertical distribution of mussel (Mytilus edulis) larvae shifted upward on flood tides on all surveyed dates at sites in the eastern Gulf of Maine. We then used a high‐resolution coastal circulation model to predict the vertical distribution of passive particles at one site on four sample dates. The modeled distribution shifts explained 99% of the variance in the larval vertical distribution. Finally, we conducted a multipart meta‐analysis of published papers on tidal shifts to assess whether variation in the results among studies could be explained by tidal vertical advection via two proxies—tidal range and shelf slope. For well‐replicated studies that could be analyzed with a paired t‐test, tidal range explained 78% of the variance in the test statistic. A multiple regression approach applied to distributional shifts from a broader set of studies indicated significant effects of tidal range and interactions between tidal range, slope, and slope heterogeneity. Our results demonstrate the importance of considering passive processes when evaluating how the vertical distribution of weakly swimming meroplankton varies between tidal phases.

The dispersal of most sedentary and sessile marine invertebrates occurs during a relatively brief larval phase in which oceanographic processes exert substantial control over where larvae are transported (Cowen and Sponaugle 1980). Numerous laboratory, field, and modeling studies extending back > 50 yr have explored the extent to which larvae may participate in defining that dispersal path (e.g., Bayne 1964;Tilburg et al. 2012;Criales et al. 2015). A few taxa possess strongly swimming larvae that may be capable of actively resisting horizontal advection (Chia et al. 1984;Young 1995;Kingsford et al. 2002). However, the slower swimming larvae of most taxa are restricted to moving vertically in the water column and altering their dispersal trajectory by migrating between water masses that differ in velocity direction or magnitude (Young and Chia 1987;Kingsford et al. 2002).
Much of the vertical migration work on taxa with weakly swimming larvae focuses on bivalves. Despite their limited swimming ability (generally on the order of 10 −4 m s −1 for sustained swimming, though speeds of 10 −3 m s −1 are possible over shorter distances; Chia et al. 1984;Metaxas 2001), observed vertical distribution shifts have typically been interpreted as reflecting behavioral movement (Carriker 1951;Verwey 1966;Harding et al. 1986;Gregg 2002;Knights et al. 2006;Peteiro and Shanks 2015;see Korringa 1952). This interpretation rests on either implicit or explicit assumptions that the vertical water advection (commonly assumed to be on the order of 10 −4 -10 −5 m s −1 ) is negligible compared to larval swimming speeds (Chia et al. 1984;Metaxas 2001;Lloyd et al. 2012). Yet this assumption has remained largely untested, in part because vertical velocities are low and hence often difficult to directly measure in the field (Arthur 1965;Halpern and Freitag 1987;Schott and Leaman 1991).
Although vertical velocities can be generated through a variety of processes (Young and Chia 1987), they have received less attention than the stronger horizontal components of currents. Our focus is on the vertical advection created when horizontal tidal currents are deflected by sloping bottom topography (a subset of case 4 in Young and Chia 1987; essentially a tidally driven case of topographically induced upwelling or downwelling- Janowitz and Pietrafesa 1982). Tidal currents flowing over an upward sloping shelf will generally produce positive vertical velocity components on a flood tide and negative vertical velocities on an ebb tide (Manasrah et al. 2006). Bivalve larval distributions typically shift in the same direction as the likely prevailing vertical velocities-upwards on flood tides and downward on ebb tides (Carriker 1951(Carriker , 1961Kunkle 1957;Knights et al. 2006;Peteiro and Shanks 2015). Consequently, while shifts in vertical distribution have been interpreted as reflecting active behavior, the direction of those shifts is also consistent with expectations of passive transport.
Past work has explored several aspects of passive larval transport (within a tidal phase, Gregg 2002; in higher velocity upwelling and downwelling currents, Genin et al. 2005; via the vertical mixing associated with very high horizontal velocities, Peteiro and Shanks 2015), though without directly examining the effects of tidal vertical advection. A related mechanism (tidal straining, in which the tidally phased breakdown of stratification leads to greater horizontal mixing) has been proposed for more restricted conditions (requiring relatively high freshwater input and estuarine circulation; Simpson et al. 1990). Others have acknowledged possible roles of vertical advection without empirically addressing those effects . We argue that the passive larval transport commonly occurs due to tidal vertical advection and that vertical currents of sufficient magnitude to affect the distribution of weakly swimming larvae are more prevalent (at least in many regions) than previously suspected.
We explored the effect of tidal vertical advection on the vertical distribution of bivalve larvae through a combination of field sampling, simulation modeling, and a meta-analysis. This approach reflected an evolution of our thinking about the problem, and each component was devised in response to the results from the preceding component. First, we quantified the vertical distribution of mussel (Mytilus edulis) larvae on flood vs. ebb tides at five sites in the Gulf of Maine (GoM) on eight different sample dates. Second, for a subset of those samples, we modeled concurrent changes in the vertical distribution of passive particles between ebb and flood tidal phases using a high-resolution coastal circulation model originally developed for modeling larval transport (Conlon et al. 2018) and compared the observed vertical distributions from field samples with the distribution of our modeled particles. To evaluate the relevance of our empirical and modeling results to other bivalve studies, we conducted a three-part metaanalysis in which we analyzed the shifts in vertical distribution for most published bivalve studies and used proxies for tidal vertical velocity components (tidal range, shelf slope, and the heterogeneity of shelf slope) to evaluate the possibility of passive effects. To establish the validity of the proxies used in our meta-analysis, we also used our circulation model to explore the extent to which the local vertical velocities within the model domain could be predicted by variation in shelf slope and tidal range.

Materials and methods
Vertical distribution of mussel larvae in the Gulf of Maine Study sites and species We sampled at five different sites in Western and Machias Bays in the northern GoM during the summers of 2012-2014 ( Fig. 1). The GoM is a relatively shallow, semi-enclosed gulf with 3-7 m semidiurnal tides, with the greatest tidal range occurring just northeast of our study region in the Bay of Fundy (Fig. 1). Circulation in both sampled bays is driven primarily by tidal exchange, with relatively little freshwater input from river discharge and only minor contributions from variation in wind (Conlon et al. 2018). Density gradients primarily reflect temperature, rather than salinity differences (Tilburg et al. 2012;Conlon et al. 2018). The larger scale coastal circulation is dominated by the Eastern Maine Coastal Current (EMCC), which can act as an ecological barrier to the onshore/offshore transport of bivalve larvae (Tilburg et al. 2012;Yund et al. 2015). Our sampling was conducted well inshore of the EMCC frontal boundary and likely sampled larvae produced in the same or neighboring bays.
All Mytilus spp. larvae in our samples were assigned to M. edulis. However, our study sites were located just southwest of the southern range boundary for the congeneric Mytilus trossulus (Rawson et al. 2007) and we cannot completely exclude the possibility that some larvae from this congener may have been present. Unfortunately, our use of electron microscopy to identify mussel larvae rendered samples unsuitable for subsequent molecular analysis.

Sampling techniques
Larval vertical sampling took place in July and August of 2012-2014, shortly after local M. edulis populations spawned, during daytime ebb and flood tides in 20-30 m of water. To collect larvae, 100 L of seawater were pumped (via an intake positioned at each sampling depth) through a 50 μm plankton net and samples were preserved in ethanol. This collection volume typically yielded hundreds to thousands of bivalve larvae at depths where they were expected to be abundant. Samples from 2012 and 2013 were collected from five water depths (0 m, 3 m, 6 m, 9 m, and 12 m), with three replicate samples collected at each depth (for a total of 15 on each tide/ date/station combination). In 2014, samples were collected from three adjacent depths in three consecutive depth bins (upper-3 m, 4 m, and 5 m; mid-9 m, 10 m, and 11 m; lower-15 m, 16 m, and 17 m). A sea anchor was deployed to minimize vessel movement in relation to the sampled water body. Each complete set of samples was collected within the middle hour of either the flood or ebb tide. The physical structure of the water column was quantified immediately before and after each larval sample with a conductivity, temperature, and depth instrument (YSI Castaway).

Sample post-processing
Bivalve veligers were sorted and enumerated under a microscope. To identify mussel larvae, the hinge tooth morphology of a subsample of veligers (approximately 33 individuals per sample) was analyzed via scanning electron microscope (SEM). Soft tissue was dissolved with bleach and veligers were identified to family level via overall hinge teeth morphology (Le Pennec 1980). Each mytillid was subsequently classified to genus based on further details of hinge tooth morphology (Fuller and Lutz 1989), with Mytilus distinguished from the much less common Modiolus via established genus-specific relationships between valve length and the number of hinge teeth (Lutz & Hidu 1979). Finally, the fraction of each sample composed of M. edulis from SEM analysis was multiplied by the total veliger count to yield the total number of M. edulis at each depth/tide/station/date combination. The weighted mean depth (WMD) of M. edulis abundance across all sample depths per date/tide/station combination was calculated as the depth of the center of mass (Eq. 1; denoted Zcm in Raby et al. 1994;Hoffle et al. 2013), where D j is the midpoint of sample stratum j, WD j is the width of the individual sample stratum, and A j is the abundance of larvae.
WMD and values derived from it were the primary quantities used in most subsequent analyses because of their relative insensitivity to variation in sample distribution (both among our own sample years and among published studies in the meta-analysis). However, un-averaged depth-specific abundance data from our samples were compared to water column density data to assess the distribution of larvae relative to stratification and mixing.

Statistical analysis
The vertical distribution of larvae in the water column was described by fitting a spline function to depth-specific densities and then comparing the resulting curve to similarly analyzed water density data (downcast only). The relationship between WMD and tidal phase was assessed with paired ttests. Relative densities were log (x + 1) transformed prior to analysis to satisfy normality assumptions.

Comparison to passive particles
We compared the observed vertical distribution of bivalve larvae in a subset of our field samples with the predicted distribution of purely passive particles generated via a threedimensional circulation model. The vertical distributions of passive particles at the time of 2014 sample collections were estimated via a high-resolution coastal circulation model validated with acoustic Doppler current profiler data and coupled with an offline particle tracking model (Conlon et al. 2018). Vertical resolution was 20 sigma layers distributed in a depthdependent manner with greater resolution near the surface and bottom (Conlon et al. 2018). Horizontal resolution also varied with depth and ranged from 100 m near shore to 4 km at the open boundary (Conlon et al. 2018). An in depth spatio-temporal validation of the circulation model using current profile data from transects within the model domain (including both Western and Machias Bays) ensured that the modeled flow patterns were accurate (Conlon et al. 2018). Model skill for transects ranged from 0.99 to 0.60 (on a scale of 0-1, with 1 being a perfect match between observed and modeled) in the east-west current direction (average of 0.87 AE 0.04 SE) and 0.99-0.62 in the north-south direction (average of 0.83 AE 0.04), with the lower skill scores due to localized bathymetric errors; the model nevertheless predicted horizontal velocities very well across the remainder of those transects (Conlon et al. 2018). We were not able to directly validate the vertical velocities from the model using the current profile data because the vertical velocities were too small relative to the noise in the profiler data. However, the accuracy of the modeled horizontal velocities lends confidence to the modeled vertical velocities because the two should be correlated.
To simulate the movement of passive particles, we coupled the circulation model with a particle tracking model and released passive, neutrally buoyant particles at the surface in areas inshore from the coastal current 4 d before each sampling period (July 22, and August 5-6). The 4-d period gave the particles time to disperse throughout the water column and develop vertical and horizontal distributions that reflected purely passive processes, instead of requiring a priori assumptions about initial vertical distributions. We analyzed only the change in particle depth for particles passing through Western Bay during the model hour that was closest to the sampling time and calculated ΔWMD (the shift in WMD between ebb and flood tides) for the passive particles. Last, we used a linear regression to evaluate the extent to which sampled larval ΔWMD could be predicted by modeled ΔWMD.

Meta-analysis Approach and justification
We compiled and analyzed published shifts in the vertical distribution of bivalve larvae for effects of tidal phase. To assess possible passive effects, we used easily quantified proxies for tidal vertical advection. As tides interact with coastal topography, it is inevitable that some vertical advection occur, and the magnitude of the vertical velocity, w, is proportional to: where α is the topographic slope, U is the horizontal velocity, A is the amplitude of the tide, H denotes water depth (or height), and g denotes gravity. Tidal range and the angle of the slope (as derived from bathymetric data) are available or can be calculated for most sampled locations and can serve as useful proxies for variation in tidal vertical advection among locations. For the third component of the meta-analysis, we also included a third proxy, slope heterogeneity, as a measure of the extent to which local variation in the slope might complicate the effects of the idealized slope.

Proxy validation
To validate the use of slope and tidal range as proxies for tidal vertical advection, we used our coastal circulation model to test for relationships between predicted vertical velocities and tidal range and shelf slope at a random subsample of nodes, dates, and tides within our model domain. Our model domain includes greater than average variation in both tidal range (3-7 m average tides) and shelf slope, and overlaps with conditions at a large subset of the meta-analysis sites. However, the lowest tidal ranges in the meta-analysis fall outside of our model conditions. Data from a random subset (n = 1214) of 174,612 model nodes were included, with nodes spaced at least 6 km apart to avoid duplicate topographic effects on adjacent nodes. Each node was sampled on one randomly selected date and tide. Shelf slope was estimated for each node via the average of five parallel transects of 6 km length oriented across-shelf, centered on the node, and separated by 1 km so that the total area for slope estimation was a 4 × 6 km rectangle. Any rectangles that included dry areas were discarded (n = 156). We calculated the modeled tidal range (change in sea surface height) at each node and obtained vertical velocity estimates by averaging the absolute value of the modeled vertical velocity for the middle of the water column during that same tide. Effects of shelf slope and tidal range on modeled tidal vertical velocities were evaluated with a multiple regression.

Meta-analysis inclusion criteria and procedures
To compile a database of published bivalve vertical distribution studies, we searched the Web of Science and Google Scholar for the terms "bivalve vertical migration", "bivalve vertical movement," and "bivalve vertical distribution." Additional sources were acquired from references cited in papers obtained via these searches. We reviewed over 300 titles, 109 abstracts, and 67 papers, which were narrowed to 11 papers that met the following three criteria: (1) field samples were collected from at least two water depths during ebb and flood tides, (2) data were reported in a quantifiable form, and (3) sufficient information was reported to determine the coordinates for each study site. Data from these 11 studies from seven distinct regions within the Atlantic Ocean (the Irish Sea, the Gulf of St. Lawrence, the Scotian Shelf, the Gulf of Maine, the New Jersey coast, Chesapeake Bay, and the Carolina Outer Banks) were combined with our own original data to yield 12 distinct studies. Many of these studies sampled multiple taxa, and 11 within-taxon comparisons distributed among the 12 studies included sufficient replication to support paired analyses of tidal effects. The remaining taxa and studies contributed additional data for a combined analysis across all taxa. For a detailed outline of study inclusion and sample sizes for each analysis, see Supporting Information Supplement A.
We extracted WMD values (if available) directly from each paper using data digitization programs (GraphClick 3.0.3 [Arizona Software] and the open platform ImageJ; Schindelin et al. 2015). If data were not presented as WMD at different tidal or diel phases, we extracted raw data from tables or figures and calculated WMD. Data points that fell within 1 h of midtidal phase were categorized as ebb or flood tides samples. Although all sites are classified as having semidiurnal tides, we lack independent evidence that peak tidal flow occurred during the middle of the tidal phase and non-barotropic tidal mechanisms are beyond the scope of this analysis. If sampling occurred twice during 1 h of the mid-tide or if samples were collected within 1.5 h before and after the mid-tide, we averaged the two WMD values. Only daytime samples were included to avoid diel effects. Tidal range was extracted directly from figures or from tide tables (using NOAA harmonics and tidal data within U.S. waters or the Mr. Tides application [version 4.3.7; Hahn Software LLC] globally) to calculate the actual tidal range on the sample date. One study (Knights et al. 2006) did not supply the actual sample dates and so we substituted the average tidal range for that site. Maximum depth at each sampling site was recorded from the text where possible, and otherwise determined using bathymetric data (NOAA data available at https://www.ngdc.noaa.gov/mgg/bathymetry/hydro. html, last accessed 18 October 2017; Syvitski et al. 1987).
We generally calculated the local shelf slope around each sampling site by averaging the slope of five parallel bathymetric transects in a 4 × 6 km rectangle centered on the site. However, at a few nearshore sites, the complexity of the bathymetry or proximity to shore required modifications of this approach, as described in Supporting Information Supplement B. A third proxy, topographic heterogeneity for each site, was approximated using the average variance in the slope residuals among the five transects.

Meta-analysis statistics
All of the component parts of the meta-analysis tested changes in WMD in response to a variety of factors (i.e., tidal phase, vertical velocity proxies), but differences in the intent of each analysis required calculating WMD in two slightly different ways. To compare generally across studies and sites (i.e., the effects of tidal phase), we removed the effect of different water column depths among studies and sites (which varied from 1.5 m in Carriker 1951 to 66 m in Tremblay and Sinclair 1990b) by calculating proportional WMD (pWMD, Eq. 3).

pWMD = WMD Maximum Water Column Depth ð3Þ
By contrast, when we were interested in assessing the actual magnitude (in meters) change in depth of larval distributions (e.g., in response to vertical velocities), depth changes were calculated from WMD as previously defined (Eq. 1).
To test the effects of tidal phase on larval distribution depth across studies, between-group differences in ln(x + 1) transformed pWMD were first calculated and the results backtransformed to attain ΔpWMD (the shift in proportional vertical distribution). Three analyses were conducted with ΔpWMD values. We first tested the simple hypothesis that ΔpWMD was greater than zero (ebb deeper than flood) with a Wilcoxon signed rank test. ΔpWMD values of 0 were replaced with a value marginally greater or less than 0 (i.e., 0.000001) to satisfy analysis requirements.
Because this first step yielded a nonsignificant result, we next asked whether variation in physical conditions among study sites might affect statistical outcomes. In the 11 study and taxa combinations with sufficient sample sizes to permit within-study statistical analyses, the hypothesis that ΔpWMD was greater than zero was tested for each comparison with either a one-tailed t-test or Wicoxon signed rank test, depending on the outcome of a Shapiro-Wilk normality test. Only one study (Andrews 1983, for bivalves) did not pass the normality test (p = 0.0315), but t-test results were consistent with the nonparametric equivalent (Wicoxon signed rank test, V = 10, df = 6, p = 0.7656; t-test, t = −1.0297, df = 6, p = 0.8286) and were included for comparison purposes. To evaluate effects of vertical velocity proxies on variation in outcome among studies and taxa, the t-statistic from each comparison was regressed against both the mean tidal range (natural log transformed to linearize the data) and the shelf slope (cube root transformed), with linear regressions weighted by the degrees of freedom in the t-test.
Finally, to more fully assess the contribution of vertical velocities proxies to changes in WMD in the larger dataset (regardless of the level of replication in the individual study), flood WMD was subtracted from paired ebb WMD, yielding ΔWMD. Standard least squares regression was used to assess the influence of site tidal range (averaged across sampled tides and natural log transformed), bottom slope (cube root transformed), and slope heterogeneity (cube root transformed) on ΔWMD. Two outliers were removed to normalize the ΔWMD distribution (Tremblay and Sinclair 1990a, sample from 05 October 1985; Tremblay and Sinclair 1990b, sample from 01 October 1985; Shapiro-Wilk test before removal, W = 0.88, p < 0.0001; after removal, W = 0.99, p = 0.63). Supporting Information Supplement A outlines the data sources used in each analysis.

Larval vertical distributions in the Gulf of Maine
Mussel larvae in our samples were generally distributed below the pycnocline and changes in vertical distribution occurred beneath the mixed layer (Fig. 2). The pycnoclines were sharper on the 2012 and 2013 sample dates then on the 2014 sample dates because of lower surface salinity (Fig. 2). Treated as a group, the WMD of mussel larvae in our 2012-2014 samples was deeper on the ebb than the flood tide (one-way paired t-test; 2012-2013, t 4 = 8.31, p = 0.0006; 2014, t 3 = 4.40, p = 0.01, Bonferroni adjusted α = 0.017; Fig. 3). This tidal shift in the depth distribution was consistent across five sampling sites, 3 yr, and two bays ( Fig. 1 and Supporting Information Supplement A).
The four 2014 samples occurred during a time period when we were able to model the expected shift of passive particles circulating in Western bay. Sample particle tracks (Fig. 4) indicated substantial vertical movement, some of which occurred in synchrony with the tides, though bathymetric effects were also apparent. There was a very close correspondence between the vertical movement of passive particles and real mussel larvae, with the modeled shifts in Fig. 2. Continued on next page.
the WMD of neutral particles explained 99.1% of the variation in the observed shifts in the weighted mean depth of larvae (ΔWMD; Fig. 5, F 1,2 = 223.78, p = 0.004), suggesting that passive processes accounted for the vast majority of vertical movement observed in our samples.

Meta-analysis
To support our use of slope and tidal range as proxies for vertical advection, we first explored the relationships between tidal range, slope, and tidal vertical velocity within our coastal circulation model domain. A multiple regression with tidal range and slope as independent variables was highly significant and explained 58% of the variation in vertical velocity (Table 1). Both main effects and the interaction effect were highly significant (Table 1). Any process that alters the magnitude of the horizontal flow will also change the magnitude of the vertical flow. Variation in horizontal density gradients and wind are likely the main factors contributing to the remaining variation in the horizontal flow and hence the vertical velocities in this system (Conlon et al. 2018).
We initially asked whether, as one large group, the bivalve studies analyzed from the literature revealed tidal shifts in larval vertical distribution (ΔpWMD). Although the observed shift was in the expected direction, it was not significant (t-test; t 85 = 0.293, p = 0.77: normality assumption met; Shapiro-Wilk test, W = 0.355). We postulated that the absence of a tidal effect across the whole dataset was due to a mixture of effect sizes in the component studies, and further explored sources of this variation. Sampling for a subset of 11 taxa was sufficiently well replicated (> 2 paired ebb vs. flood measurements) to support within-study comparisons with a t-test (Table 2). Two of the 11 taxa/study combinations revealed significantly deeper pWMD on ebb than flood tides, while two others exhibited a marginally nonsignificant trend in that direction, and the remainder showed no clear effect (Table 2). Among this group, 77.6% of the variation in the t-statistic (regardless of significance level) was explained by Dashed lines show a spline (λ = 0.05) fit to water density data as an average of two casts conducted before and after larval sampling. Solid lines represent a spline (λ = 0.05) fit to M. edulis proportional abundance throughout the water column. Black circles represent larval data points. Note that the scale of the x-axes differs among panels to more accurately depict variation in the depth-specific relative density of mussel larvae among dates. mean local tidal range during the sampling periods (one of our proxies for tidal vertical advection; Fig. 6; F 1,9 = 31.14, p = 0.0003). The qualitative result did not change even if one point that disproportionately influences the regression (i.e., the highest value in Fig. 6) was removed as a potential outlier (F 1,8 = 10.44, p = 0.0120). By contrast, only 17.6% of the variation in the t-statistic was explained by local bottom slope and that relationship was not significant (F 1,8 = 1.71, p = 0.23). However, this analysis had less power than the tidal range analysis because one taxon/study combination  for bivalves) had to be excluded because of insufficient bathymetric data for slope calculation. In addition, the relationship with slope could have been complicated by variation in along vs. cross-isobath tidal flow among sites.
We then asked whether the effect size in the entire group of studies, regardless of level of replication within the individual study, could be explained by our tidal vertical advection proxies. A least squares regression model for all ΔWMD comparisons (regardless of within-study level of replication) explained 29.4% of the observed variation and revealed significant effects of tidal range, tidal range by bottom slope   interaction, bottom slope by topographic heterogeneity interaction, and the three-way interaction among tidal range, bottom slope, and topographic heterogeneity (Table 3). Of these effects, tidal range accounted for the vast majority of the variation in ΔWMD (Table 3, sums of squares 2-4 times greater than for other significant effects). A positive tidal range parameter estimate (1.74 AE 0.45 SE), indicated that the magnitude of the difference between ebb and flood WMD increased with tidal range. Variance inflation factors (Neter et al. 1989) indicated minimal multicollinearity among these independent variables, but prevented us from adding a water column depth effect (per Eq. 2) to the model.

Discussion
Vertical distribution shifts in M. edulis larvae Within our study region, tidal shifts in the vertical distribution of larvae appear to occur reliably and repeatedly (Fig. 2) and involve vertical movement below the pycnocline, as has been reported in other studies (Raby et al. 1994;Lloyd et al. 2012). However, distributions centered above the pycnocline have also been noted in some regions (Dobretsov and Miron 2001).
The modeled tidal vertical velocities in our study system ( Fig. 4 and Conlon et al. 2018) were substantially higher than commonly assumed (on the order of 10 −3 -10 −4 m s −1 compared to an assumed 10 −4 -10 −5 m s −1 ; Metaxas 2001; Lloyd et al. 2012) and comparable to the magnitude of larval bivalve swimming speeds (Chia et al. 1984). While the existence of vertical velocities of this magnitude in our study region will not be surprising to the readers with a background in physical oceanography, we have found no evidence in the larval vertical migration literature that vertical velocities are expected to vary systematically with tidal range or shelf slope. When our circulation model was used to generate passive expectations for passive particle depth associated with tidal changes, the modeled depths explained 99% of the variance in the empirical data from the same dates (Fig. 5). This strong relationship suffers from the inevitable limitations of a small sample size (i.e., only 4 points). Nevertheless, the magnitude of that association was strikingly high and inspired us to conduct the meta-analysis to further explore a passive physical explanation.

An argument for passive tidal vertical distribution shifts
Although numerous other authors (Carriker 1951(Carriker , 1961Kunkle 1957;Booth 1972;Maru et al. 1973;Gregg 2002;Knights et al. 2006;Peteiro and Shanks 2015) have reported the same tidal shifts in bivalve larval distributions that we demonstrated empirically (Fig. 3), we differ in supporting a primarily passive, physical interpretation, and suggest that evidence for an active behavioral component may currently be lacking. We argue that upward and downward vertical advection are inevitable byproducts of flooding and ebbing tides (respectively), that the magnitude of these vertical velocities is sufficient to explain much of the tidal shift in vertical distribution where tidal shifts have been observed, and that variation in the magnitude of tidal vertical advection explains why tidal vertical shifts have been demonstrated in some studies, but not others.
We develop our argument based on a number of lines of evidence that differ in the strength of inference supplied. Two lines of evidence are consistent with a passive mechanism, but do not rule out other processes. First, the reported tidal vertical shifts are virtually always in the direction of passive movement (up on flood tides and down on ebb, but not the reverse). Hence there is widespread evidence that larval distributions shift in the direction of tidal advection, but with no consensus on the magnitude or role of advection. Second, the pattern of larvae remaining beneath the pycnocline (Fig. 2) is consistent with a passive mechanism for tidal shifts. If larval distributions were routinely observed to shift across the pycnocline between tidal phases, it would be difficult to invoke passive processes. Vertical advection could not carry larvae through an observed pycnocline, because vertical advection in that depth range would result in a proportional vertical shift of the pycnocline.
Other lines of evidence provide more direct support for a passive mechanism and reduce the need to invoke a behavioral explanation. On the admittedly limited number of only four dates for which we were able to model the distributional shift of passive particles, passive movement driven by tidal vertical advection explained virtually all of the variance in mean larval depth between flood and ebb tides (Fig. 5). Although our study region experiences substantial tides, vertical velocities of approximately comparable magnitude should be characteristics of all regions with relatively large tidal ranges and steep shelf slopes. Published studies with data that support significant tidal shifts in the vertical distribution of bivalve larvae were all conducted in regions with at least moderate tidal ranges and slopes (Mid-Atlantic-Carriker 1951, 1961Kunkle 1957;Wood and Hargis 1971;Andrews 1983;Gregg 2002;Baker and Mann 2003: Sea of Japan-Maru et al. 1973: North Atlantic-Knights et al. 2006: Eastern Pacific-Peteiro and Shanks 2015. Vertical larval distribution shifts on the order of 1-2 m appear to be consistent with tidal vertical advection in most of these systems. Our meta-analysis supports extending this passive interpretation beyond our immediate study system to all bivalve larvae. First, when analyzed as a group, the aggregate set of bivalve studies did not indicate a consistent tidal vertical shift. When variation in tidal effect sizes among the different studies was subsequently explored with paired tests (for the subset of studies with adequate in-study replication to support such analyses), tidal range (a proxy for tidal vertical advection) explained 78% of the variance in the test statistic (Fig. 6). This pattern indicates that the aggregate group of studies can be decomposed into subsets exhibiting significant and nonsignificant effects (Table 2), with a proxy for tidal vertical advection largely determining whether a significant effect was detected (Fig. 6). If behavior played a major role in tidal shifts, we would expect to see significant effects in lower tidal range systems. In a multiple regression analysis of all bivalve studies (regardless of the level of replication in the individual study), tidal range and the interaction effects between tidal range and bottom slope (a second proxy for tidal vertical advection) as well as the three-way interaction between tidal range, bottom slope, and slope heterogeneity (which serves as a measure of local variance in slope, or the extent to which slope alone may be an imperfect predictor of tidal vertical advection) were all significant predictors of tidal shifts in vertical distribution ( Table 3). The relevance of tidal range and bottom slope as proxies for tidal vertical velocity is supported by first principles (see Eq. 2), but also demonstrated by significant relationships to vertical velocities in our empirical study region with predictions from the circulation model (Table 1).
We are hardly the first to question a behavioral interpretation for tidal larval vertical distribution shifts and argue for a passive explanation. As early as 1952, Korringa (writing about oyster larvae) was concerned about the absence of any apparent behavioral mechanism for tidal migration: "I have never understood how a certain group of larvae, swimming in a certain body of water, can distinguish whether it is ebb or flood" (Korringa 1952, p. 303). As Korringa recognized, there are very few candidate cues linked to tidal phase that would operate on a spatial scale that is accessible to larval sensory systems (Young 1995). Temperature and salinity cues were postulated in early years (e.g., Carriker 1961) but appear to be relevant only in limited circumstances because a larva has to encounter the cue (e.g., by moving into a thermocline or halocline) before it can respond. Thus, these cues are unlikely to trigger larger scale distributional shifts that occur throughout many meters of water column. A decade after Korringa (1952), the  Table 2) regressed on natural log transformed mean tidal range. Solid best fit line flanked by dotted 95% confidence intervals; R 2 = 0.78, F 1,9 = 31.14, p = 0.0003. discovery that bivalve larvae can respond to variation in pressure (Bayne 1963) offered a possible explanation. However, barokinesis does not link uniquely to tidal phase and has not really been embraced as the answer. It is common for bivalve larvae to temporarily stop swimming, and because they are denser then water, stationary larvae start to sink (Roughley 1933;Isham and Tierney 1953;Konstantinova 1966;Lough and Gonor 1971). An increase in pressure due to sinking is thought to be the cue to resume swimming to maintain a fixed depth (Raby et al. 1994;Young 1995). Consequently, even papers invoking behavioral mechanisms for tidal shifts in vertical distribution sometimes explicitly state that the cue is unknown (e.g., Lloyd et al. 2012). Nevertheless, stronger swimming crustacean larvae clearly demonstrate tidal shifts in vertical distribution (Queiroga 1998;Forward and Tankersley 2001;Criales et al. 2015), in spite of being subject to the same sensory scaling issues. Differences in behavioral vs. passive effects among taxa may partly explain why larvae of different taxa tend to co-localize in the field (Daigle et al. 2014a,b).
Although we did not explicitly analyze other taxa, we strongly suspect that the proposed passive mechanism will apply to tidal vertical distribution shifts in other taxa with weakly swimming larvae. For example, gastropod larvae swim at approximately the same speeds as bivalve larvae (Chia et al. 1984). Yet a recent paper interprets tidal shifts in the vertical distribution of gastropod larvae as prima facie evidence of the existence of a behavioral mechanism (Kunze et al. 2013).

Issues and limitations
If taxonomic information (Supporting Information Supplement A) is compared with the results from our paired analysis of well-replicated studies (Fig. 6), a potential pattern emerges. The two studies with the highest test statistics for tidal vertical shifts (values above 2.0 in Fig. 6) both involved mytilid mussels. However, in our more extensive multiple regression analysis (Table 3) Raby et al. 1994;Gregg 2002). Further, some Mytilus vertical shifts fell below the overall average, and one was even slightly negative (Knights et al. 2006). Overall, it seems unlikely that Mytilus alone exhibit a tidal signal in vertical distribution while other bivalve species do not. More likely, Mytilus exhibits a stronger signal because the species' ranges, and thus study locations, include higher latitudes with greater tidal ranges and/or steeper shelf slopes.
A number of published studies include observational comments to the effect that older larvae exhibit different vertical shifts than younger larvae (Carriker 1951;Kunkle 1957;Verwey 1966;Cragg 1980;Dobretsov and Miron 2001;Ishii et al. 2004). Unfortunately, neither our own empirical sampling nor the meta-analysis yielded sufficient data on specific size classes to evaluate these claims. Consequently, the empirical evidence for ontogenetic effects remains thin (Knights et al. 2006).
Several earlier studies reference tidal larval movement in or out of the sediments and/or benthic boundary layer, generally via passive sinking (Andrews 1983;Baker and Mann 2003;Knights et al. 2006). Our data and analysis were never intended to address larval distribution shifts at the very bottom of the water column. Accurately sampling larval distributions on that scale is very challenging and may require specially designed sampling technology (Knights et al. 2006). Thus, our interpretations apply only to the overall distribution of larvae through the bulk of the water column, but not in and around the benthic boundary layer.

Conclusions
Overall, our results question the need to invoke behavior to explain tidal shifts in the vertical distribution of bivalve and other weakly swimming larvae. The reported shifts in mean depth are typically fairly small, on the order of 1-2 m, and within the range that can be explained if tidal vertical velocities are on the order of 10 −3 -10 −4 m s −1 , as we believe they are in many steeply sloped, high tidal range systems. Nevertheless, the existence of passive processes that shift larval distributions upward on flood tides and downward on ebb tides does not negate the possible existence of a behavior operating in those same directions. However, detecting such a behavior (if it exists) will require much more carefully controlled studies that evaluate the magnitude of passive advection and subtract that transport from the net tidal distribution shift to estimate any remaining behavioral component. A simple tidal shift in distribution (Kunze et al. 2013) should not be interpreted as evidence of behavior in any taxa with weakly swimming larvae.