Eﬀects of Ocean Climate on Spatiotemporal Variation in Sea Urchin Settlement and Recruitment Relationships between larval settlement and biotic and abiotic conditions Effects of maternal and larval nutrition on growth and form of larvae. A Site Profile of the Kachemak Bay Research Reserve, a Unit of the The role of fluctuating food supply on recruitment, survival and population dynamics in sea among eggs shifts to cooperation along a sperm "Spiral asters" and cytoplasmic rotation in sea urchin eggs: induction in Strongylocentrotus purpuratus eggs by elevated Lagrangian descriptions of marine larval Supplementary Information - Appendices for “Geographically Opposing Responses of Sea Urchin Recruitment to Changes in Ocean Climate”

1 Sea urchins are voracious herbivores that influence the ecological structure and function 2 of nearshore ecosystems throughout the world. Like many species that produce planktonic 3 larvae, their recruitment is thought to be particularly sensitive to climatic fluctuations in 4 temperature that directly or indirectly affect adult reproduction and larval transport and 5 survival. Yet how climate alters sea urchin populations in space and time by modifying larval recruitment and year-class strength on the time-scales that regulate populations remains 7 understudied. Using an unprecedented, spatially replicated weekly-biweekly dataset spanning 8 27 years and 1100 km of coastline, we characterized seasonal, interannual, and spatial patterns 9 of larval settlement of the purple sea urchin ( Strongylocentrotus purpuratus ). We show that 10 large spatial differences in temporal patterns of larval settlement were associated with 11 different responses to fluctuations in ocean temperature and climate. Importantly, we found a 12 strong correlation between larval settlement and regional year class strength suggesting that 13 such temporal and spatial variation in settlement plays an important role in controlling 14 population dynamics. These results provide strong evidence over extensive temporal and 15 spatial domains that climatic fluctuations shape broad-scale patterns of larval settlement and 16 subsequent population structure of an important marine herbivore known to control the 17 productivity and function of marine ecosystems.


19
Large scale climate oscillations (e.g., El Niño/Southern Oscillation, North Atlantic 20 Oscillation) lead to changes in ocean temperature, biogeochemistry and the severity and 21 frequency of disruptive events that affect ocean circulation, upwelling and primary productivity 22 (Cai et al. 2014;Mantua et al. 1997). Such shifts impose wide-reaching ecological impacts, in 23 part by altering animal recruitment and food web structure in space and time (Sydeman et al. 24 2015). Hence, understanding how climate variability alters the recruitment of marine species is 25 particularly important for effective conservation and management of the ocean's resources. 26 Yet for recruitment variability, climatic fluctuations give rise to shifts in numerous factors that 27 shape both adult reproduction and larval supply, including primary productivity, temperature, 28 and advection and transport. Thus, significant challenges remain in achieving such 29 understanding for benthic species with planktonic larvae due to the substantial effort needed 30 to characterize spatial and temporal variation in larval settlement and the numerous sensitive 31 vital rates that contribute to it. 32 33 For benthic species like sea urchins, understanding causes and consequences of 34 recruitment variability has both ecological and economic implications. Sea urchin grazing can 35 alter the structure of some of the world's most diverse and productive marine ecosystems, 36 including coral reefs (Edmunds & Carpenter 2001), seagrass meadows (reviewed by Valentine & 37 Heck Jr 1999) and kelp forests (reviewed by Filbee-Dexter & Scheibling 2014). In addition, sea 38 urchins form the basis of important nearshore fisheries in many regions of the world (e.g. 39 Andrew et al. 2003;Kato & Schroeter 1985). As a result, climate-driven changes in sea urchin 40 populations have the potential to profoundly affect the functioning of marine ecosystems and 41 the economic value of the fisheries that they support. Much of the research on controls of sea 42 urchin population dynamics has focused on the roles of predation and disease in controlling 43 adult abundance and their cascading influence on community structure (e.g. Burt et al. 2018;44 Estes & Duggins 1995;Filbee-Dexter & Scheibling 2014;Lafferty 2004). Yet short-term 45 empirical studies (months to a few years) suggest that environmentally regulated larval supply 46 is likely an important driver of adult urchin dynamics (Hernández et al. 2010;Ling et al. 2009). 47 Despite the widespread recognition of the importance of recruitment variation in controlling 48 population fluctuations in many marine species (Shelton & Mangel 2011), relatively few studies 49 have examined the biotic and abiotic processes controlling the supply of sea urchin larvae in 50 nature (but see Hernández et al. 2010;Ling et al. 2009), and the degree to which they affect the 51 abundance and dynamics of older life stages over time scales that impact population dynamics. 52 53 Fluctuations in climate can affect spatial and temporal patterns of larval supply by 54 influencing the production of larvae by benthic adults, transport of larvae to adjacent habitats, 55 and the survival of larvae in the plankton. Increases in ocean temperature can impact larval 56 production by: (1) increasing adult mortality via the spread of water-borne pathogens 57 (reviewed by Feehan & Scheibling 2014), and (2) reducing adult fecundity and inhibiting 58 gametogenesis by altering food quantity and quality (Basch & Tegner 2007;Cochran & 59 Engelmann 1975;Foster et al. 2015;Okamoto 2014). Because sea urchins produce feeding 60 larvae that spend weeks to months in the plankton, increases in ocean temperature can also 61 affect larval development, growth and survival, either directly, or indirectly by altering the 62 by incorporating a Bayesian prior in the form of a Beta distribution with hyperparameters as the 129 observed ratio and total subsample count (the Jeffrey's prior for ratios from binomial counts). 130 Third, initial examination of the data indicated that there were strong temporal components 131 (interannual and seasonal trends, substantial and multiplicative among-sample variability 132 within each period, and serial autocorrelation) and spatial components that needed to be 133 simultaneously accounted for. 134 To account for the above described nuances of the data we estimated biweekly 135 settlement density ( ",$ ) as the sum of a site-specific mean (β &,$ ), log annual trend 136 (Annual_Trend 2," ), log seasonal trend (Seas_Trend ",$ ), and a spatially correlated lognormal 137 process error (ε ",$ ); 138 Eq. 1 ln ",$ = β &,$ + ln Annual_Trend ",$ + ln Seasonal_Trend ",$ + ε ",$ (Trend Equation) 139 We estimated trends on the log-scale so that process errors, and seasonal and interannual 140 trends were multiplicative and strictly positive on the original scale. We used a Poisson 141 observation likelihood to link the biweekly mean estimate from the trend equation ( ",$ -see 142 Tables 1 and 2) to observed counts of larval settlers (N) for each brush (b) within each site 143 within each year . 144 The annual trends were assumed to be correlated in both space and time, where spatial 145 and temporal covariances were assumed to be independent. The annual spatial covariance was 146 assumed to be unstructured (i.e. no formal distance decay structure) because we only had 147 seven sites. The annual temporal covariance was determined by a Matérn 3/2 kernel because it 148 provided sufficient flexibility to capture interannual trends. The seasonal trend (log-scale) 149 within each site was estimated using a seasonal periodic temporal kernel within the model to 150 capture the cyclical nature of seasonality. The process error (ε ",$ ) was modeled with an 151 unstructured spatial covariance to account for correlations among sites in their deviations from 152 model expectations. We estimated the model posteriors using Stan (Carpenter et al. 2016) with 153 three 1000 iteration chains after a 1000 iteration burn-in. Stan model code is provided in the 154 supplement. For model equations and symbology see Tables 1 and 2.  155 156

Relationship between benthic year-class strength and larval settlement. 157
If year class strength of sea urchins is limited by larval supply, then we should see 158 increases (decreases) in the abundance of juveniles in years following anomalously high (low) 159 larval settlement. We therefore tested whether larval settlement in the Santa Barbara Channel 160 (mean of Anacapa, Stearns Wharf, Ellwood Pier and Gaviota Pier) was predictive of subsequent 161 juvenile recruitment (hereafter referred to as "year-class strength") on natural reefs in the 162 region. We calculated densities of juveniles (<2cm test diameter) from the Channel Islands Kelp 163 Forest Monitoring (KFM) Program (Kushner et al. 2013), including only the XXX sites with time 164 series extending from 1990 through 2016. We examined the relationship between larval 165 settlement and year-class strength using a generalized linear mixed effects model (GLMM) with 166 a negative binomial likelihood (a Poisson likelihood indicated overdispersion) and survey site as 167 a random effect. We tested the hypothesis of no correlation using a likelihood ratio test 168 comparing the model with versus without settlement as a covariate. Models were estimated 169 using glmmTMB (Magnusson et al. 2017) in R (see Electronic Supplement for further details). 170

Relationships between larval settlement and biotic and abiotic conditions 172
We estimated the strength of relationships between biweekly larval settlement and 173 various hypothesized physical and biological drivers using an integrated regression model that 174 included a regularized regression within a Bayesian time-series modelling framework (Fig. 2). 175 Specifically, the model considered estimated log biweekly density (ln ",$ ) as a function of the 176 centered and standardized covariates (denoted by the vector , ) while directly estimating and 177 accounting for seasonal trends (Seasonal_Trend ",$ ), and separable spatially and temporally 178 correlated, multivariate lognormal process error (ε ",$ ): 179 Eq. 2 ln ",$ = &,$ + , ′ + ln Seasonal_Trend ",$ + ε ",$ (Regression Equation) 180 We used a Poisson observation likelihood to link the latent trends ",$ to the observations (see 181 Tables 1, 2) and we assumed the temporal correlation in the process error followed a first order 182 autoregressive model. 183

184
Because these analyses were correlative, we focused on unbiased parameter estimation 185 rather than hypothesis testing or model comparison per se. We estimated two separate 186 models: one included all locally measured (i.e. for a given site) environmental variables and the 187 other included a common set of composite indices of oceanographic climate ("global 188 covariates" -ENSO Index, Pacific Decadal Oscillation, North Pacific Gyre). Because inclusion of 189 numerous covariates can cause overfitting that leads to bias and uncertainty in the explanatory 190 power of the covariates, we erred on the side of sparsity and assigned the vector of regression 191 coefficients ( ) a regularized horseshoe prior (i.e. the Finnish Horseshoe). Sparsity in this case 192 assumed that only a few of the covariates were meaningful without a priori knowledge of which 193 covariates were relevant and which were not. A sparse regression encoded this assumption to 194 allow the data to inform which covariates were relevant and how they were correlated with the 195 response variable. The Finnish horseshoe prior is a Bayesian version of the lasso (Carvalho et al. 196 2010;Piironen & Vehtari 2017) that produced a data-driven reduction in the influence of 197 weaker covariates by regularization of those coefficients given the data. Oscillation (NPGO -Di Lorenzo et al. 2008). The MEI provides a metric of the intensity of El 206 Niño Southern Oscillation (ENSO) fluctuations, which persist for 6-18 months and explains much 207 of the oceanographic variability in the tropics (Di Lorenzo et al. 2013;Wolter & Timlin 1993;208 Wolter & Timlin 1998). In contrast, the PDO and NPGO exhibit longer-period fluctuations, 209 explain much of the oceanographic variability in higher latitudes, and the low frequency 210 variability of these metrics are shaped through ENSO forcing of the Aleutian Low and the North 211 Pacific Oscillation, respectively (Di Lorenzo et al. 2013). 212 213 (ii) Coastal upwelling index (monthly, 1997-2016, all sites): The Bakun index (Bakun 214 1973) provides an index of large-scale coastal upwelling and specifically describes the volume of 215 water that is transported offshore from Ekman transport (http://www.pfel.noaa.gov/ -sites = 216 33N-119W and 39N-125W. This index has been used as a large-scale proxy for processes that 217 may favor larval retention and advection from shore in addition to its value as a predictor of 218 coastal productivity (Menge et al. 2011;Menge & Menge 2013;Shkedy & Roughgarden 1997;219 Wing et al. 2003). We emphasize that this index is not a location specific metric nor does it 220 consider fine-scale nearshore hydrodynamic processes that are also likely to affect larval 221 retention (Fisher et al. 2014;Morgan et al. 2009;Morgan et al. 2016;Morgan et al. 2018 (monthly, 1997-2016, all sites): Satellite imagery of sea 228 surface chlorophyll a provides a spatially and temporally resolved estimate of phytoplankton 229 biomass (mg m -3 ) that is not available from in situ sampling. We used version 3.1 of the OC-CCI 230 merged ocean color time series (Sathyendranath et al. 2018) that combines SeaWIFs, MERIS, 231 MODIS and VIIRS to provide the temporal and spatial coverage required for this study. For each 232 larval settlement collection site, we aggregated data into a 30-day moving average (30 days 233 prior to brush collection) within a buffer of 10 km from shore that stretched 150 km alongshore 234 (the average Lagrangian estimate for dispersal distances for species with a planktonic larval 235 duration (PLD) of 30-days (but see Shanks 2009;Siegel et al. 2003)). For the site at Anacapa 236 Island, we included any point within a 150 km radius and within 10 km of any coastline. We 237 used 30-day moving averages for chlorophyll and temperature (below) because larvae were 238 exposed to these conditions for at least 30 days prior to settlement, and because averaging 239 over 30 days minimized the effects of serial autocorrelation in the data. 240 241 (iv) Sea surface temperature (monthly, 1997-2016, all sites): We used the 30-day moving 242 average of sea surface temperature, derived from Pathfinder AVHRR (Reynolds et al. 2007) 243 (advanced very high resolution radiometer) that was optimally interpolated at daily and 0.25 244 degree latitude/longitude resolution. We spatially aggregated data in the same way as sea 245 surface chlorophyll. 246

247
(v) Fall kelp canopy biomass (annual, 1996-2015  settlement corresponded with a nearly three-fold increase in the mean density of juvenile 318 urchins two summers later ( Figure 5). Over the time series, average settlement densities in the 319 Santa Barbara Channel varied by more than three orders of magnitude among years. This 320 interannual variation in settlement (averaged across the sites in the Santa Barbara Channel) 321 corresponded to, on average, more than a three-fold increase in the density of benthic 322 juveniles. 323

Relationships between larval settlement and biotic and abiotic conditions 325
Larval settlement of purple urchins showed strong correlations with local sea surface 326 temperature (SST, Figure 6 a-g, Figure 7 a) as well as the major climate indices (e.g. ENSO - Fig 6  327 h-n, Figure 7f). For temperature and ENSO, the sign of the correlation at Fort Bragg was 328 opposite of that at sites in the Santa Barbara Channel and San Diego. The correlation between 329 settlement and temperature varied from strongly negative at sites in San Diego and Santa Fairweather 1989). In the case of sea urchins, larval supply has been linked to dramatic 384 changes in the state of benthic communities (Estes & Duggins 1995;Hernández et al. 2010;Ling 385 et al. 2009). Yet empirically demonstrating such links has historically proved challenging 386 because recruitment can be attenuated by high mortality of newly settled larvae (Connell 1985;387 Rowley 1989), increasing the need for high frequency data of larval settlement spanning 388 multiple years. Data such as ours that meet these criteria are rare. 389

390
The potential factors affecting larval supply that we examined other than temperature 391 and climate showed no meaningful correlations with larval settlement. Of particular note was 392 the lack of a correlation with ocean chlorophyll revealed by our analyses. This finding differed 393 from that observed for the tropical sea urchin heterogeneity not only influences larval transport, but also when and where larvae encounter 407 productive food environments. Consequently, if larval food limitation plays a role in shaping 408 the dynamics of urchin settlement at our sites, then it likely arises out of more complex 409 dynamics than our spatially aggregated composite metrics were able to detect. For example, 410 processes that shape stratification, front formation and spatially isolated phytoplankton blooms 411 in combination with larval behavior may play a more important role than predicted by 412 simplifications of the higher or lower spatially averaged phytoplankton production or greater or 413 lesser coastal upwelling. Integrating quantitative measurements of these processes was 414 beyond the scope of our study, yet such approaches are needed to develop a mechanistic 415 understanding of processes that control settlement dynamics across large gradients in 416 oceanographic settings. 417 418 Although climate related changes can alter larval supply via their effects on the 419 production of larvae by adults, we found no relationship between larval settlement and 420 regional adult abundance or adult food (kelp). While fisheries research relies heavily on stock-421 recruit dynamics, there is continual debate about whether adult dynamics actually control 422 recruitment patterns (Gilbert 1997;Szuwalski et al. 2015). The lack of a positive correlation 423 between adults and larval settlement that we found for purple sea urchins adds to this debate 424 by showing that high abundances of adult urchins and kelp averaged across large spatial scales 425 did not translate into high larval supply. We note that spatial structure in the dynamics of kelp 426 and adult urchins may have affected this outcome. For example, modest but uniform densities 427 of adult sea urchins in food rich kelp forests can provide very different outcomes for the 428 production of embryos compared to mosaics of dense adults in barren patches interspersed 429 with forested patches having few adults (Okamoto 2014;Okamoto 2016). Indeed,previous 430 work has demonstrated that assessing population productivity in spatially structured 431 populations requires a spatially structured analysis to account for processes of population 432 regulation and competition (Chesson 1996(Chesson , 1998Chesson et al. 2005;Thorson et al. 2015). Niño caused anomalously strong poleward flow and decreased local retention in coastal areas 455 off southern California (Lynn & Bograd 2002;Mitarai et al. 2009). Yet El Niño events also create 456 stochastic small-scale patterns of water motion in the form of fronts, eddies, and bores (Pineda 457 et al. 2018). Spatial variability in such patterns could explain the opposing climatic responses in 458 larval settlement between Fort Bragg and our sites in southern California. El Niño events in 459 southern California are associated lower stratification and reduced internal waves that promote 460 the onshore transport of larvae (Pineda 1994;Pineda et al. 2018;Shanks 1983). In contrast, El   Larval Settlement