Kay, Fisher, & Hoekstra, 2012) and genotyped 791 recruits and 1,729 adult A. clarkii at 1,340 146
Single Nucleotide Polymorphisms (SNPs) using the bioinformatics pipeline dDocent (Puritz, 147 Hollenbeck, & Gold, 2014) (see Supporting Information). We then used Colony2 to identify 148 parent-offspring matches (Wang, 2012) . Colony2 uses maximum likelihood to jointly estimate 149 parentage and sibship relationships from SNP data. We ran Colony2 for each year by including 150 only recruits sampled in that year. We cumulatively pooled adults across years so that once an 151 adult was sampled, it was included as a possible parent in the analysis of that and all subsequent 152 years (Table S1 ). We retained parentage assignments found in the maximum likelihood where Γ denotes the gamma function. The generalized Gaussian is defined by a scale parameter k 165 that sets the dispersion of the distribution and a shape parameter θ that affects its kurtosis (how 166 much of the density is contained in the tail). At a value of θ = 1, the kernel is Laplacian and 167 leptokurtic, exhibiting exponential decline with increasing distance and a fat tail. In contrast, θ = 168 2 is a Gaussian kernel, while θ = 3 is a platykurtic Ribbens function with a broad central peak 169 and a thin tail. 170 We calculated mean dispersal distance, as well as the distance at which 50% (median) 171 and 90% of larvae are expected to recruit for each kernel. Then, we calculated a distribution of 172 mean and median distances from 1,000 parameter sets of k and θ drawn with replacement from 173 each kernel 2D likelihood surface, weighting parameter sets by the difference between the log-174 likelihood value and the maximum log-likelihood value (parameter sets closer to the maximum 175 log-likelihood are thus more likely to be drawn). Finally, we integrated the kernels across the 176 width of our study region (30 km) to estimate the proportion of recruits retained in the 177 metapopulation. 178
To examine variation in kernel shape, we quantified the excess kurtosis (the density in 179 the tails relative to a Gaussian kernel) using the following equation (Clark et al., 1998) and June). In contrast, recruits greater than 4.5 cm but less than 6 cm were likely to be 6 to 11 193 months old and to have dispersed during the Southwest Monsoon. 194 We tested whether the observed interannual and seasonal variation in kernel fits was 195 greater than expected from a null model that simulated no differences among years or seasons. 196
We randomly shuffled the years in which recruits were observed, creating a seven-year time 197 series that broke apart any relationship between year and dispersal kernel. We fit kernels to each 198 of the seven years, calculated the standard deviation of k across years and calculated the 199 coefficient of variation for θ, mean dispersal distance, and median across years. We then 200 repeated the reshuffling 1,000 times. We compared these null distributions of temporal variation 201 to the empirical estimates. We used the same approach for seasonal dispersal variation. 202
Finally, for both the annual and seasonal kernel fits, we explored the sensitivity of our 203 results to the proportion of the population sampled. We down-sampled the all-years data (2012-204 2018 for both assigned and unassigned recruits) at proportions ranging from 0.01-1.0 and fit a 205 kernel to each of n = 100 samples. We then plotted the relationship between the proportion of 206 data sampled and either k or θ. 207 208
Next, we tested whether average dispersal direction changed among years. We fit 210 generalized linear models (GLMs) with binomial errors where the response variable indicated 211 whether a parentage match was observed (1) or not (0) for each pairwise site-to-site connection 212 (dispersal route) each year (D' Aloia et al., 2013) . Predictor variables included the distance 213 between sites, the sampling effort at the source (adult) and destination (recruit) sites, the source 214 and destination site size (number of anemones), the sampling year, and the direction from source 215 to destination (North = +1, same site = 0, South = -1). When larvae recruit to their natal (source) 216 site, we refer to this as self-recruitment. We calculated sampling effort as the proportion of 217 anemones sampled at a site in each year (see Supporting Information). To test directional 218 switching of dispersal among years, we included an interaction term between year and 219 directionality. We used Akaike information criterion (AIC) to sequentially remove terms from 220 the full model. We used the same approach to explore variation in direction between seasons, 221 using season as a predictor variable rather than year. 222 223
Connectivity between pairs of sites (dispersal routes) underlies metapopulation dynamics 225 (Hanski, 1994; Hastings & Botsford, 2006; Roughgarden & Iwasa, 1986 ). To examine temporal 226 variation in each dispersal route, we discretized the annual kernels according to the known 227 widths of and distances between sites, creating a connectivity matrix for each year where each 228 cell (a single dispersal route) represented the probability of dispersal from a source site (column) 229 to a destination site (row). We then calculated the coefficient of variation for each dispersal route 230 across the seven years. The parentage analysis identified 71 parent-offspring matches from 791 recruits and 241 1,729 adults for an overall assignment rate of 9% (Table S2 ). We also identified parents for 11 of 242 132 offspring that likely dispersed during the Northeast Monsoon and for 35 of 428 offspring 243 from the Southwest Monsoon (Table S3 ). Of the assigned offspring, 28 were assigned to both 244 parents. When both parents were identified, they had always been sampled on the same 245 anemone. Overall, the frequency of observed dispersal events generally decreased as distance 246 increased (Fig. S1 ). Across all years, 27 individuals dispersed northward, 23 dispersed 247 southward, and 21 recruited to their natal site (Fig. 2) . 248 249
The best fit parameter sets describing the scale (k) and shape (θ) of the dispersal kernels 251 differed among years (Fig. 3 , Table S2, Table S4 ), though the bivariate likelihood surfaces 252 overlapped for some years ( Table S2 ). Mean dispersal distance point estimates 259 differed almost ten-fold across years (CV = 1.03), and these differences were greater than 260 expected from the null model (p = 0.001, Fig. S3C ). 261
Median dispersal distances ranged from a minimum of 5 km in 2013 to 12 km in 2012 262 (CV = 0.34), corresponding to a two-fold difference among years (Fig. 4B ). In every year, 50% 263 of larvae recruited within 15 km or less of their natal origin (Fig. 4B , Table S2 ) but the variation 264 among years was greater than expected from the null model (p = 0.041, Fig. S3D ). The distance 265 at which 90% of larvae recruited ranged from a minimum of 15 km to a maximum of 157 km, 266 while the proportion of recruits retained in the 30 km study site was between 0.33-0.49 across 267 years (Table S2) . 268 Driven by differences in kernel shape, we observed variation in the excess kurtosis 269 among kernels ( (Table  284 S3). 285
Fitting kernels to down-sampled datasets demonstrated that k and θ were relatively 286 insensitive to samples larger than 50% of our dataset, but that k and θ were biased high for 287 smaller datasets (Fig. S5, Fig. S6 ). 288 289
The most parsimonious regression model suggested that the probability of dispersal 291 between any two patches declined with distance, was higher for dispersal from north to south, 292 increased with the size of destination and source habitat patches, and varied among years (Table  293 S5). There was some evidence for including recruit sampling effort (DAIC = 1.8) as well as adult 294 sampling effort (DAIC = 2). The interaction between year and direction was not included in any 295 of the models with DAIC < 5, suggesting no evidence for change in the predominant direction of 296 dispersal among years (Table S5) . 297
The best model for the seasonal kernels was similar, but there was some evidence for 298 variation in dispersal direction among seasons (Table S6) Table S6B ). 303 304 Correlation among dispersal routes 305
Self-recruitment to sites was low, from a minimum of 0.23% to a maximum of 8.6%. 306
However, dispersal probabilities along individual routes had high interannual coefficients of 307 variation that were generally highest for self-recruitment routes (range = 0.45-1.0, mean = 0.53, 308 We also found strong correlations between dispersal routes. Coefficients ranged from -313 0.94 to 1, with a mean correlation of 0.21 ± 0.67 (SD) (Fig. 6B ). Clustering identified three main 314 groups (Fig. S7 ). The first cluster included routes between sites that were both geographically 315 close (including self-recruitment) and distant (-0.28 ± 0.61). Within this cluster, routes between 316 close sites were positively correlated, and distant routes were negatively correlated. The second 317 cluster included routes between sites that were somewhat close to each other and were more 318 weakly correlated (0.23 ± 0.53). The third group was of routes between sites that were 319 intermediately distant and mostly positively correlated (0.36 ± 0.66). 320 Levin, 2006) . Here, we directly quantified 325 dispersal patterns across seven years in a coral reef fish and found substantial variation in 326 dispersal kernel scale and shape across years. We also observed seasonal variation in dispersal 327 that was of similar magnitude to interannual variation. Though there was temporal variation in 328 dispersal distance, we found some evidence of variation in direction between monsoon seasons 329 but did not find strong evidence for variation in the dispersal direction among years. Self-330 recruitment at the scale of individual sites and of the metapopulation varied strongly through 331 time, and dispersal probabilities between sites were often correlated. 332
Marine larval dispersal has been called the "black box" of marine ecology (Robert K. While we studied a marine system, seasonal variation in dispersal is also likely in 391 terrestrial systems. In deciduous forests, the seasonal loss of foliage can increase the potential for 392 long-distance dispersal by allowing vertical eddy motions to transport seeds higher in the canopy 393 where winds are stronger (R. Nathan & Katul, 2005) . Broadly, seasonal variation in reproduction 394 should increase variability in connectivity patterns, not only because environments can differ 395 between seasons (e.g., monsoons or storm frequency), but also because seasonal reproduction 396 will sample a smaller range of environmental conditions and therefore be more sensitive to 397 chaotic or stochastic environmental variation (Siegel et Although these results are the most complete time series of dispersal kernels of which we 434 are aware and cover more than a generation in our focal species, seven years remains limited. 435
Future work assessing longer time series, ideally on decadal scales that may illustrate the 436 importance of long-term climatic oscillations, is an important though logistically challenging 437 next step. There is also a need for studies of dispersal variability across diverse environments and 438 taxa to assess the generality of these results. Additionally, we described dispersal using dispersal 439 kernels because they provide a convenient framework for considering dispersal across spatial 