Methods Over 3 years, we measured annual oyster recruitment to standardized ceramic tiles on 9-16 intertidal reefs in coastal Virginia, USA. We paired these measurements with an 18-year time series of recruitment to natural substrate on 8 natural reference reefs and 44 restored reefs (0-16 years post-construction). Results Recruitment to tiles was highly correlated with recruitment to natural substrate, validating our methodology. Recruitment was positively spatially autocorrelated within 1 km and increased 9-14 × with moderate wind fetch. A one-meter increase in substrate elevation tripled recruitment. Recruitment was 4 × higher on natural reefs compared to restored reefs, regardless of elapsed time since restoration. Geospatial model predictions identified 6% (24 km 2 ) of intertidal areas as highly suitable for oyster recruitment, offering a refined target for restoration practitioners. Conclusions By integrating multi-year field studies, long-term monitoring, and habitat suitability modeling, our research identified environmental conditions favorable for oyster recruitment, offering insights that could enhance restoration planning and population resilience. Our findings provide actionable insights for optimizing oyster restoration by targeting areas with favorable wind fetch and elevation. These results offer valuable guidance for spatial planning in restoration and may inform strategies for other species where recruitment limits restoration success.

Restoration projects have developed globally to reverse extensive declines in habitat-forming foundation species and the ecosystem services they provide (Coleman & Williams 2002;Löf et al. 2019;Duarte et al. 2020). Long-term success in restoration relies on restoring the full life cycle of the targeted species to ensure self-sustaining populations (Hastings & Botsford 2006;Vanderklift et al. 2020;Temmink et al. 2021). A common bottleneck in habitat restoration is the failure to recruit sufficient offspring of the focal species. Recruitment is defined here as the first appearance of macroscopic individuals, reflecting settlement and post-settlement processes (Jenkins et al. 2009;Cruz & Harrison 2017;Toone et al. 2023). Recruitment failure can result factors such as nutrient pollution, turbidity (Cleland et al. 2013;Toone et al. 2023), habitat characteristics like topography and substrate availability (Whitman & Reidenbach 2012;Wang et al. 2020), restoration history (e.g., restored vs. naturally occurring substrate, elapsed time since restoration; Golet et al. 2008;Catterall 2018), and pressures such as predation or competition for space, food, or settlement area; Grabowski et al. 2005, Fodrie et al. 2014). Such recruitment bottlenecks have contributed to the mixed success of restoration efforts (Cruz & Harrison 2017;Esquivel-Muelbert et al. 2022;Larson et al. 2023), perpetuating the need for continuous recovery and incurring high restoration and monitoring costs (Bayraktarov et al. 2016). Addressing these challenges requires a comprehensive understanding of the environmental drivers of recruitment variation and effective methods for assessing outcomes, informing adaptive management, and learning from restoration successes and failures (Baggett et al. 2015;Eger et al. 2022). Achieving this goal often hinges on obtaining consistent and representative measurements of recruitment using standardized techniques that facilitate comparisons among multiple locations of differing environmental conditions.

Substantial restoration efforts have focused on re-establishing oysters globally following geographic declines during the twentieth century, especially the eastern oyster (Crassostrea virginica) along the U.S. Gulf and Atlantic coasts (1768 projects, 5199 ha restored, $299,999 per ha, on average; Hernández et al. 2018;Duarte et al. 2020). However, oyster restoration outcomes vary greatly (Schulte et al. 2009;Beck et al. 2011;Smith et al. 2022a;Smith & Castorani 2023). Differences in oyster reef restoration trajectories are often attributed to recruitment constraints, including dispersal limitations due to insufficient larval supply (recruitment limitation) or a lack of suitable settlement substrate (habitat limitation) (Lipcius et al. 2008;Schulte & Burke 2014). As for other species with dispersing propagule stages (e.g., plant seeds, invertebrate larvae), overcoming these constraints is critical to achieving oyster restoration goals (Atwood & Grizzle 2020;Vanderklift et al. 2020).

Because of stochastic dispersal and demographic processes, recruitment is inherently variable across time and space (Caffey 1985;Menge 1991;Jenkins et al. 2009). Spatial autocorrelation analyses are a valuable tool for quantifying this variability and identifying drivers (Koenig 1999;Chevalier et al. 2021;Ford et al. 2021). In the context of oysters, hydrological forces (Lenihan 1999;Knights et al. 2012;Hubbard & Reidenbach 2015) and elevation (Baillie & Grabowski 2019;Johnson et al. 2019) shape recruitment and reef development across multiple spatial scales, from fine-scale reef elevation changes (meters) to broader hydrodynamic patterns affecting larval dispersal (kilometers to tens of kilometers). Higher water flow may enhance larval delivery and extend the duration available for settlement and attachment (Kimbro et al. 2009;Byers et al. 2015;Hubbard & Reidenbach 2015), but it can also transport oyster larvae to unsuitable habitats or resuspend sediment, smothering settling larvae (Thomsen & McGlathery 2006;Whitman & Reidenbach 2012;Baillie & Grabowski 2019). Despite these dynamics, gaps remain in understanding how oyster recruitment varies across spatial scales, particularly in relation to hydrological forces and substrate elevation (Byers et al. 2015;Baillie & Grabowski 2019). Moreover, uncertainty remains about how environmental factors and restoration history together determine where oyster populations experience sufficient recruitment to sustain reefs and support successful restoration.

To address these gaps, we investigated the spatial patterns of intertidal oyster recruitment and how they are determined by geophysical conditions and reef restoration history (natural vs. restored and, if restored, the age of the restoration) by integrating a multi-year field study with a long-term oyster reef monitoring survey. First, we validated a standardized method for measuring recruitment by comparing recruitment densities on ceramic tiles with those on natural reef substrate (comprising live oysters and shell hash). Second, we used both datasets to: (1) investigate the spatial patterns of oyster recruitment and geophysical conditions (elevation, wind fetch, and water residence time, defined below) using spatial autocorrelation analyses; (2) quantify the effects of these environmental conditions on oyster recruitment; (3) determine how recruitment varies among natural reference reefs and restored reefs at different stages of recovery (as represented by the elapsed time since restoration); and (4) use model results to identify intertidal regions with suitable habitat for high oyster recruitment that can be targeted for future restoration efforts. Our study uncovers the spatial scales and patterns of variation in oyster recruitment and the accompanying geophysical conditions that facilitate this process, providing valuable insights for oyster restoration management and spatial planning for ecosystem restoration more broadly.

We studied populations of eastern oysters (C. virginica, hereafter 'oyster') within the Virginia Coast Reserve (VCR), a set of interconnected bays along the Atlantic coast of Virginia, USA, and site of the VCR Long-Term Ecological Research project. VCR salinities are euhaline (typically > 30 PSU) due to limited freshwater inputs. Low nutrient loading results in oligotrophy and high-water quality (McGlathery et al. 2001;McGlathery & Christian 2021). Water circulation is driven by semi-diurnal tides (1.2-1.5 m) that flow through deep inlets to connect bays with the ocean (Safak et al. 2015). The complex VCR landscape (~ 950 km 2 ) of barrier islands, marshes, and lagoons causes highly heterogeneous wind-wave energy (measured as wind fetch distance) and tidal water flushing (measured as water residence time) (see Geophysical measurements, below; Wiberg et al. 2015, Hogan & Reidenbach 2019, Besterman et al. 2021). Coastal development is minimal, with agriculture being the most significant land use, and contemporary wild oyster harvest is extremely limited (McGlathery et al. 2001, Virginia Marine Resources Commission 2022). Oysters are distributed throughout the VCR on intertidal fringing and patch reefs (and to a lesser extent, manmade structures), and our surveys confirm the absence of subtidal reefs.

Despite relatively pristine conditions and nowthriving oyster populations, the VCR has a legacy of oyster reef restoration following severe overharvesting and disease that decimated oysters in the 1900s (Smith et al. 2022a), mirroring declines in the nearby Chesapeake Bay (Schulte et al. 2009). Since 2003, The Nature Conservancy (TNC) and Virginia Marine Resource Commission have constructed over 24 ha of intertidal oyster reefs throughout the VCR, primarily by depositing natural reef substrate. A recent oyster habitat suitability model-parameterized using the current distribution of reefs and data on wind fetch, water residence time, and elevations surrounding the reefs-suggests that 52.4 km 2 of the VCR (or 12% of the 436.4 km 2 examined intertidal and subtidal area) may be suitable for oyster reefs (Hogan & Reidenbach 2019). Validation of this model with independent long-term oyster monitoring data showed 1.5 × greater adult oyster biomass in areas with predicted high-suitability habitat (Smith et al. 2022b). However, it is unknown whether such areas are also associated with high oyster recruitment, or why oyster reefs are found only in a small fraction of the predicted suitable area (Ross & Luckenbach 2009). Moreover, although mature (> 7 years) restored reefs match the density and biomass of adult oysters (> 25 mm shell height) on natural reference reefs that serve as restoration targets (Smith et al. 2022a), it is not known whether restored reefs of varying age (time since construction) recruit equivalent numbers of oyster larvae as reference reefs.

To investigate spatial patterns of oyster recruitment, we measured recruitment to standardized ceramic tiles and reef substrate at dozens of intertidal reefs throughout the VCR. We defined recruitment as the density of macroscopic oysters ≤ 25 mm shell height attached to the sampling unit (ceramic tile or natural reef substrate, see below), reflecting the joint influence of settlement and post-settlement processes (Jenkins et al. 2009).

First, in a standardized field study, we selected sites along the margins of the mainland, shallow lagoons, and barrier islands supporting either patch reefs or narrow bands of fringing reefs. These sites were selected to capture the spatial variability of reefs in the region, ensuring a representative range for the analysis. In this study, we measured recruitment to 100 cm 2 square unglazed, textured ceramic tiles at 9-16 reef sites in 2018, 2019, and 2021, with some sites revisited in multiple years. We deployed 3-10 tiles per reef site (n = 218 total tiles across all 3 years) with 3 tiles per site in 2018 due to initial study design, and an average of 9.6 tiles per site in 2019 and 2021 to increase replication (some tiles were lost due to foul weather). Deployments occurred in in early-tomid May, shortly before the typical start of the oyster spawning season in Virginia (Kennedy & Krantz 1982, Hubbard & Reidenbach 2015). At each site, we affixed tiles to PVC poles distributed haphazardly across the reef (1-2 tiles/pole, with poles separated by at least ~ 3 m), oriented vertically at the height of live oysters, with the rough, unglazed surface of the tiles facing outward (away from the poles) and the bottoms of the tiles at or slightly above the mud of natural reef substrate (Fig. 1a-b). Because tiles varied in their absolute elevation and across-reef location at each site, our data characterized site-level variation in recruitment to the reef, rather than recruitment to a particular microhabitat. After 3.5 months (in midto-late August), we retrieved all tiles and immediately photographed their outward-facing, unglazed surfaces in the laboratory. Prior to photography, we lightly scrubbed each tile with a toothbrush to remove algae, taking care not to remove oysters attached to the tiles. We manually counted attached oyster recruits in the photographs using image analysis software (ImageJ; Schneider et al. 2012).

Second, we analyzed data from long-term monitoring of oysters on restored and reference reefs. Data were collected by TNC from 2005 to 2022 across 8 reference reef sites and 44 restored reef sites (Lusk et al. 2022). During each sampling event, a 0.0625 m 2 quadrat excavation (15-30 cm depth) was collected by hand or a small trowel-sized hand rake from haphazard locations across the reef (1-3 replicates per site; n = 925 total excavations across all 18 years) and all live oyster recruits were counted (Fig. 1c-d). The date and frequency of sampling varied among sites and years (Lusk et al. 2022). a b

To test the validity of our standardized recruitment measurement approach, we compared recruitment density estimates per unit of area sampled between ceramic tiles (no. recruits per cm 2 of tile) and oyster shell from the reef substrate (no. recruits per cm 2 of bottom). To make appropriate comparisons, we selected data from TNC-surveyed reef sites located < 250 m from 8 locations where we measured recruitment to ceramic tiles in the same year (2018, 2019, or 2021). We used regression and correlation analyses to compare these estimates (see Statistical analyses).

To assess the effects of hydrodynamic conditions and elevation on oyster recruitment, we gathered existing data on wind fetch, water residence time, and substrate elevation. We used modeled wind fetch distances as a proxy for wind-wave energy and modeled water residence times as a proxy for tidal water flushing (Safak et al. 2015;Kremer & Reidenbach 2021). Fetch distances were estimated through wind measurements from the summers of 2014 and 2015 (30-m resolution; Kremer & Reidenbach 2021). Fetch was calculated as the maximum unobstructed distance over water in a constant direction, considering observed wind directions from a nearby meteorological station. Despite variability in fetch between 2014 and 2015 due to changes in wind direction, the fetch lengths at sampling sites were highly correlated (r 2 = 0.92). Thus, we assumed that the average fetch lengths between the two years represented typical site-specific fetch distances (Hogan & Reidenbach 2019). Water residence times were modeled using a validated, three-dimensional Lagrangian ocean circulation model (FVCOM;

Chen et al. 2006) with particle tracking (200-m resolution; Safak et al. 2015). We created wind fetch and water residence time raster layers in QGIS 3.20.3 (QGIS.org 2023) using model outputs and extracted data for each reef site location. We estimated substrate elevations (m relative to NAVD88) using data from a field-validated (Hogan & Reidenbach 2019) airborne-based LiDAR (light detecting and ranging) survey conducted in 2015 at low tide for each reef site location (< 1 m horizontal resolution with centimeter-scale vertical accuracy; Dewberry 2016). For reference, 0 m NAVD88 is approximately -0.04 to -0.15 m relative to Local Mean Sea Level (LMSL; Richardson 2013) and -0.15 to -0.81 m relative to Mean Low Water (MLW; Richardson 2013).

To determine how recruitment differed among reference reefs and restored reefs at different stages of recovery, we categorized each restored reef site in the TNC dataset as 'developing' (0-6 years since construction) or 'mature' (7-16 years since construction) because restored reefs in the VCR approximate the adult oyster density, biomass, and temporal stability of reference reefs in about 6 years (Smith et al. 2022a). We conducted separate analyses on the full TNC natural reef substrate data (52 reef sites in total) and a subset comprising 20 distinct restored reef sites, each paired with a nearby suitable reference. The categorical age of the restored reef site was used for further analysis and comparison with reference reef sites in multiple regression models (see Statistical analyses).

To validate a standardized method for measuring recruitment, we used correlation analysis and Deming regression to quantify the relationship between oyster recruitment density to ceramic tiles and natural reef substrate in R 4.1.1 (R Core Team 2023). We calculated the squared Pearson correlation (r 2 ) in R using the cor function to describe the strength of the relationship (Bossé et al. 2023). Deming regression is analogous to least-squares regression, but the method explicitly accounts for uncertainties in both x and y estimates by minimizing the sum of squared distances from the regression line in both directions (Fuller 1987). The relationship is significant if the 95% confidence interval of the Deming regression slope does not overlap zero (Farrer et al. 2021). This approach was appropriate because of the measurement error associated with both regressed variables. We fit the Deming regression in R with SimplyAgree 0.1.2 (Caldwell et al. 2022).

Next, addressing our first goal, we quantified how similarity in oyster recruitment density, wind fetch, water residence time, and elevation changed as a function of the distance between reef sites. We pooled oyster recruitment density estimates from ceramic tiles and the natural reef substrate, and averaged data for reef sites with repeated sampling over time (n = 130). We calculated correlations using Moran's I between all pairs of sites and created spline correlograms in R using ncf 1.3.2 (Bjørnstad 2022). A Moran's I value near -1 or 1 indicates strong spatial autocorrelation (negative for dissimilarity and positive for similarity), whereas a value close to 0 indicates weak or no spatial autocorrelation (Moran 1950). These correlations, estimating spatial dependence as a continuous function of distance (Chevalier et al. 2021), use maximum likelihood estimation with the denominator set to the sample size n for univariate data (Bjørnstad 2022). To determine how spatial correlation in recruitment may be induced by similarities in geophysical conditions, we repeated this procedure by fitting a negative-binomial generalized linear model (GLM; log link function) in R using glmmTMB 1.1.5 (Brooks et al. 2017) that predicted recruitment as a function of first-and second-order polynomials of wind fetch, water residence time, and elevation. Polynomials were necessary because preliminary models indicated nonlinearity. We then extracted model residuals and repeated the spline correlogram procedure described above to examine whether any of the spatial autocorrelation observed in the oyster recruitment data could be explained by the environmental predictors. To ensure the robustness of our findings, in addition to analyzing the pooled data we repeated the analysis separately for the ceramic tile data and the TNC natural reef substrate data (see Appendix S1: Figs. S1-2).

Addressing our second goal, to assess the effects of geophysical conditions on oyster recruitment, we fit a negative-binomial generalized linear mixed-effect model (GLMM) to predict the density of recruitment to ceramic tiles and natural reef substrate (no. per unit of area sampled) as a function of first-and second-order polynomials of wind fetch, water residence time, and elevation (random effects specified below; Zuur et al. 2009). We also estimated categorical effects of the month when the data were collected and the sampling method (ceramic tile or natural substrate). This analysis encompassed n = 1,143 observations of recruitment across all sites and years (218 measurements of recruitment to ceramic tiles and 925 measurements of recruitment to natural reef substrate). To ensure the robustness of our findings, in addition to analyzing the pooled data we repeated the analysis separately for the ceramic tile data and the TNC natural reef substrate data, dropping the irrelevant categorical predictors of month and/or sampling method as necessary (see Appendix S1: Fig. S3).

To address our third goal of determining how recruitment varied among natural reference reefs and restored reefs of differing categorical ages, we fit a similar GLMM as described above to a subset of the recruitment data for which we knew the restoration history of the reef site (n = 925 observations). These recruitment estimates were limited to TNC natural reef substrate data because the restoration history was not well documented for sites where ceramic tiles were used. We incorporated all the geophysical predictors and an additional categorical predictor describing reef restoration type and age (reference reef, mature restored reef, or developing restored reef). We calculated pairwise differences in oyster recruitment between the three reef restoration history categories using a Tukey's HSD post-hoc test with Bonferroni-corrected p-values (Underwood 1997).

For GLMMs (addressing second and third goals), we treated site as a random intercept to account for non-independence of data collected from the same reef and year as a random intercept to account for interannual variation in recruitment across all reefs (Zuur et al. 2009). We ensured that all model residuals met the assumptions of linearity, normality, homoscedasticity, and zero-inflation in R using DHARMa 0.4.6 (Hartig 2022). We square-root transformed water residence time for the model fitting due to gaps in the range of values and back-transformed values for data visualization. We estimated marginal means and 95% confidence intervals in R using ggeffects 1.1.5 (Lüdecke 2018).

Lastly, we addressed our fourth goal by developing a map that integrated the geophysical conditions (wind fetch, water residence time, and elevation) associated with above-average oyster recruitment based on model estimates. We created a prediction grid for each significant geophysical predictor identified in our models (p < 0.05) and re-fitted them with this new dataset. We calculated the mean oyster recruitment along the fitted line and identified the range of values for a given geophysical predictor that supported above-average recruitment. We integrated these ranges of geophysical conditions into multiple raster layers in QGIS, creating a new layer based on overlapping areas using the 'Intersection' tool to include only features (polygons) that overlapped among the input layers. We overlaid the new layer, which we refer to as the "suitable habitat for high or aboveaverage recruitment" layer, onto an existing oyster habitat suitability model for the VCR (Hogan & Reidenbach 2019). To determine the extent of intertidal area across the VCR, we calculated the maximum tidal range using water level data collected at a nearby NOAA station in Wachapreague, VA (Station ID: 8631044; https:// tides andcu rrents. noaa. gov) from 2005 to 2022. Building on the methods described earlier, where substrate elevations were estimated using an airborne LiDAR survey, the recorded maximum values for both high and low water levels (m relative to NAVD88) served to establish a range of intertidal elevations, representing the vertical difference between high tide and low tide. We repeated the procedure above using a 3-m resolution bathymetry model for the same range of elevations (Richardson et al. 2014). Using this additional data source ensured we had the maximum coverage of intertidal bottom. Subsequently, we used this available intertidal area (404 km 2 ) to assess the extent of overlap between the layers representing above-average oyster recruitment (this study) and adult oyster habitat suitability (Hogan & Reidenbach 2019). We also examined the overlap between mapped natural oyster reefs (Ross & Luckenbach 2009;Hogan & Reidenbach 2015) and areas with both above-average recruitment and suitable adult habitat to better understand how habitat influences oyster populations in the study region.

In validating our standardized oyster recruitment method, we found a strong correlation between recruitment to ceramic tiles and recruitment to the natural reef substrate from nearby (< 250 m) locations in the same year (squared Pearson correlation coefficient r 2 = 0.78, p = 0.002; Fig. 2). Applying Deming regression to account for uncertainties in both measurements, we found that the slope estimate (1.48) did not differ from one (t 1,7 = 1.2; p = 0.3; 95% confidence interval [CI] = 0.56-2.39) and the intercept estimate (-0.69) did not differ from zero (t 1,7 = 0.4; p = 0.4; 95% CI = -2.3-1.0).

Addressing our first goal-to investigate the spatial patterns of oyster recruitment and geophysical conditions-we found that oyster recruitment was highly variable across spatial scales, both on ceramic tiles (mean no. of oyster recruits per cm 2 of tile ± standard deviation [SD] = 2.1 ± 1.9) and on the natural reef substrate (mean no. of oyster recruits per cm 2 of bottom = 3.7 ± 5.2; Fig. 3, see Appendix S1: Fig. S4). Recruitment showed positive spatial autocorrelation at distances < 1 km (95% CI of maximum Moran's I = 0.06-0.4), but no significant trends at any farther distances (i.e., the 95% CI overlapped Moran's I = 0 for remaining spatial lags; Fig. 4a). Residuals from a GLM with wind fetch, water residence time, and elevation showed a reduction in the spatial autocorrelation at distances less than ~ 1 km (95% CI of Moran's I at distance 0 km = -0.08-0.2), suggesting that spatial similarities in recruitment at this scale were partially explained by geophysical conditions (Fig. 4b). In contrast, there was negative spatial autocorrelation of model residuals at distances between 22 and 26 km (the 95% CI of Moran's I < 0), likely an artifact due to the spatial clustering of sampling locations separated by such distances (see Discussion).

Addressing our second goal-to quantify the effects of geophysical conditions on oyster recruitment-we found that the fixed and random effects examined in our study explained 60% (conditional R 2 ) of the total variance in oyster recruitment, with wind fetch and elevation emerging as strong geophysical drivers. Recruitment had a unimodal relationship with wind fetch (χ 2 2,1123 = 34.4; p < 0.001; Fig. 5a), in which sites having moderate wind fetch (≈ 2000 m) had 9-14 × greater mean recruitment than sites with short (< 100 m) or long (4000 m) fetch distances. Recruitment was unrelated to water residence time (χ 2 1,1123 = 0.2; p = 0.6; Fig. 5b), indicating that oysters recruited across a broad range of tidal water flushing rates (residence time < 1 h to > 600 h). Recruitment was positively related to substrate elevation (χ 2 1,1123 = 32.3; p < 0.001; Fig. 5c), increasing 3 × across the approximately one-meter range in elevation (-1.0 m to 0.07 m NAVD88). Recruitment was influenced by the month of data collection Fig. 4 Spatial correlograms display spatial autocorrelation of oyster recruitment a and model residuals of oyster recruitment b after controlling for the effects of wind fetch, water residence time, and elevation. Larger positive values of Moran's I indicate greater spatial clustering of similar values, larger negative values indicate stronger spatial dispersion of dissimilar values, and values near zero indicate random spatial dispersion. The lines show the fitted spline correlogram, with shading indicating the 95% confidence intervals. The blue rectangles highlight areas where the confidence intervals do not overlap zero, denoting significant positive or negative spatial autocorrelation. Negative spatial autocorrelation at 22 km to 26 km in panel B was driven by five reefs at the southern margin of our study domain (see Aappendix S1: Fig. S7). Additional visualizations with all pairwise comparisons overlaid (see Appendix S1: Fig. S8) and recruitment density mapped by location (see Appendix S1: Fig. S9) complement Fig. 4 (χ 2 11,1123 = 37.1; p < 0.001; appendix s1: Fig. S5), but not the sampling method (tiles vs. natural substrate; χ 2 1,1123 = 0.03; p = 0.9; appendix s1: Fig. S6), consistent with the methodological validation carried out at a subset of paired sites, described above.

Addressing our third goal-to determine how oyster recruitment varies among natural reference reefs and restored reefs at different stages of recoverywe found that reference reefs had 4 × higher recruitment relative to mature (7-16 years old) and developing (0-6 years old) restored oyster reefs (mean no. recruits per cm 2 of natural reef substrate: 7.1 vs. 2.3-2.9; post-hoc pairwise comparisons: z = 3.5, p = 0.001 for mature restored reefs; z = 6.5, p < 0.001 for developing restored reefs; Fig. 6). On restored reefs, oyster recruitment on mature reefs was not significantly different from developing reefs (mean no. recruits per cm 2 : 2.9 vs. 2.3; z = 2.3, p = 0.06). In our secondary analysis with a smaller subset of spatially paired reefs, we found that reference reefs continued to have higher recruitment even when restored reefs were built close (< 2 km) to reference reefs (χ 2 2,549 = 29.5; p < 0.001).

Addressing our fourth goal-to identify potential restoration areas most suitable for (n = 1,143). The p-values indicate the respective main effect displayed from the mixed-effects model. Lines and shading show mean model predictions and 95% confidence intervals, respectively, for significant predictors (p < 0.05) after controlling for covariates oyster recruitment-we found that wind fetch distances ranging from 966 to 3413 m (mean fetch ± SD: 2160 ± 611) and elevations spanning -0.69 m to 0.07 m NAVD88 (mean elevation ± SD: -0.19 ± 0.18) were predicted to support above-average (i.e., 'high') oyster recruitment (> 2.2 mean no. oyster recruits per cm 2 tile or bottom). For the study region encompassing the total intertidal area (404 km 2 ), 24 km 2 (= 6%) fell within the range of suitable habitat for aboveaverage recruitment, as determined by wind fetch and elevation criteria (Fig. 7). When we compared the predicted areas supporting above-average oyster recruitment (totaling 24 km 2 ) with the predicted habitat suitable for adult oysters (51.7 km 2 = 12.8% of 404 km 2 encompassing only intertidal areas; Hogan & Reidenbach 2019), we found that there was much less area predicted to be suitable for both aboveaverage recruitment and adult oyster habitat (8.7 km 2 = 2.1%). We also examined the overlap between all mapped oyster reefs (natural and restored) and areas predicted to support both above-average recruitment and suitable adult habitat. We found that these existing reefs overlapped with suitable adult habitat areas (1.11 km 2 = 27.5%) slightly more than with above-average recruitment areas (0.99 km 2 = 24.5%).

We further examined the area where both above-average recruitment and suitable adult habitat overlapped existing reefs, which was more constrained (0.55 km 2 = 13.6%).

We integrated a multi-year standardized field study with long-term monitoring data to identify optimal intertidal geophysical conditions and restoration history for oyster recruitment. Our findings lead to three major conclusions with high relevance for oyster reef conservation and restoration: (1) Geophysical factors including fetch (a proxy for wind-wave energy)

Mapped oyster reefs Suitable habitat for high oyster recruitment and substrate elevation control larval recruitment for eastern oysters in our euhaline, intertidal system. Knowledge of these relationships should help optimize the spatial design of restored reefs by prioritizing building reefs in areas with fetch and elevation favorable to recruitment. Recruitment is the net result of larval supply, settlement, and post-settlement processes; understanding spatial variation in these and other demographic processes, such as survival and fecundity, will benefit oyster restoration planning (Puckett et al. 2018). ( 2) Natural reference reefs represent ideal restoration endpoints, as evidenced by the higher recruitment on reference reefs compared to restored reefs. However, restored reefs of any age (up to 16 years old)

were not yet recruiting oyster larvae equivalent to natural reference reefs. Continued study will determine whether restored reefs eventually catch up to natural ones. Both findings emphasize the importance of considering recruitment as an additional metric in evaluating restoration success.

(3) Areas ideal for both oyster larval recruitment and adult oysters should be prioritized for restoration, as areas ideal for adults are not necessarily good for recruitment, and vice versa. By incorporating spatial predictions of locations with above-average recruitment with an existing adult oyster habitat suitability model, we developed geospatial information helpful for identifying priority areas for restoration. Our particular results for coastal Virginia and general approach to informing oyster restoration planning are most applicable to recruitment-limited systems but also apply to habitat-limited systems in which additional recruitment enhances local population density. Therefore, our framework may prove useful for improving restoration science and planning for oyster populations in different regions and for other species for which recruitment limits population establishment, persistence, and size.

The positive spatial autocorrelation pattern observed in oyster recruitment at short distances (< 1 km) indicates the presence of spatial clustering of similar recruitment rates among nearby locations. After incorporating geophysical variables and examining model residuals, we found that this small-scale variability was partly attributable to the combined influence of wind-wave energy (fetch), tidal water flushing rate (water residence time), and elevation. This insight may imply that reefs in close proximity that share similar geophysical conditions tend to exhibit comparable recruitment because of 'environmental filtering' in which habitat selection and survival of oyster larvae and post-larvae drive recruitment. In theory, spatial autocorrelation in recruitment at short distances could also be controlled by similarity in oyster larval delivery and retention at scales < 1 km (Roughgarden et al. 1988). However, such distances are small relative to the scales of variation in tidally-driven water transport in our system (Safak et al. 2015); improved local understanding of oyster larval transport will help resolve this uncertainty (Lipcius et al. 2008;Theuerkauf et al. 2021). The spatial scales of similarity that we found for oyster recruitment echo those described for coral reefs, where positive spatial autocorrelation in benthic communities has been detected at scales of < 2 km (Ford et al. 2021). Thus, positive autocorrelation at such small scales may be a common spatial attribute of reef-building systems. Elsewhere, negative spatial autocorrelation can arise at similarly short distances (< 1 km), as in riparian plant communities, where competition operates at finer scales relative to favorable environmental conditions at a broader scale (Biswas et al. 2017). The detection of negative spatial autocorrelation in the model residuals of oyster recruitment (i.e., after controlling for geophysical conditions) suggests patchiness or a potential grouping structure in the system (Diniz-Filho et al. 2003, Kowe et al. 2019). Spatial separation (the geographical distance between locations) contributes to shaping recruitment dynamics where nearby sites exhibit similar recruitment patterns, whereas increased distance results in dissimilarity. This pattern-driven by the reefs at the southern margin of our study domain, situated a considerable distance away-likely arose from localized variations in recruitment responses to the environment. Notably, when these southern sites were removed from the analysis, the negative spatial autocorrelation in model residuals disappeared (see Appendix S1: Fig. S7). These nuanced responses may extend beyond the influence of the variables examined, emphasizing the need to consider additional factors that may differ for southerly (vs. northerly) sites such as the delivery and retention of oyster larvae (Lipcius et al. 2008;Theuerkauf et al. 2021); substrate suitability (e.g., sand vs. silt; Theuerkauf et al. 2017); predation pressure and the community of local predators (Tedford & Castorani 2022); or competition for space with other oysters and benthic organisms (e.g., mussels, barnacles; Boudreaux 2005;Diederich 2006).

Of the geophysical variables examined, wind fetch and elevation were the strongest drivers of oyster recruitment on the intertidal reefs that we studied. Wind fetch influences wave energy, and higher wave energy can resuspend sediments. Although we sampled relatively few reefs in areas with large windfetch distances that correspond with high wind-wave energy (9% of reefs had fetch ≥ 2500 m; Wiberg et al. 2015), most existing reefs with higher recruitment occurred in areas with moderate wind-wave energy (roughly 1500 to 2500 m). A moderate amount of wave energy may enhance sediment transport, which can create stable, suitable substrate for oyster settlement (O'Beirn et al. 1995). However, excessively high wave energy could make it more difficult for oyster larvae to settle and attach (Whitman & Reidenbach 2012), or cause sediment scouring that dislodges settled larvae or destabilizes the reef substrate, potentially contributing to the observed spatial pattern (Theuerkauf et al. 2017). We also observed diminished oyster recruitment in areas with lower windwave energy (short wind-fetch distances). Fine sediments, such as silt and clay, may accumulate quickly on reefs in such areas due to reduced water movement and greater sediment deposition, potentially smothering oyster recruits (Lenihan 1999;Housego & Rosman 2016;Baillie & Grabowski 2019). The hydrodynamic environment also modifies oyster feeding behavior and the quantity and quality of available food (Wilson-Ormond et al. 1997). Lastly, flow and mixing rates can modify densities of predators and spatial competitors on oyster reefs (Powers & Grabowski 2023). Further study will help determine the plausibility of these flow-related processes in contributing to the unimodal relationship between oyster recruitment and wind fetch in our system.

These potential flow-related mechanisms also align with our finding that recruitment is reduced at lower elevations, which correspond with intertidal reefs that are deeper, experience longer periods of inundation, and generally are subject to weaker hydrodynamic forcing (Lenihan 1999;Wiberg et al. 2019). Furthermore, if oyster predation plays a significant role in post-settlement mortality due to increased vulnerability on deeper reefs (Fodrie et al. 2014), this pattern could be explained by the fact that elevated areas provide advantages, such as reduced predation risk, in addition to reduced sedimentation.

Lastly, in contrast to the effects of wind fetch and elevation, the lack of an effect of water residence time that oyster recruitment was insensitive to a broad range of tidal water flushing rates (residence time < 1 h to > 600 h). This result was surprising, as shorter water residence times-often associated with faster water movement and increased flushing-can lead to higher nutrient levels and lower silt content (Kimbro et al. 2009;Wiberg et al. 2015). Likewise, as with many marine species that have low-motility larvae, tidal flushing rates can influence oyster larval dispersal, affect larval retention time for settlement, and impact the ability of settling larvae to attach to the substrate (Whitman & Reidenbach 2012;Schulte & Burke 2014). However, it is possible that in our system strong wind-driven vertical mixing and horizontal transport (large fetch distance) counteracted the effects of low tidal flushing (long water residence times) in some areas. In addition, hydrodynamic connectivity among oyster reefs-driven by tides, winds, waves-may have provided ample opportunities for oyster larvae to disperse and settle to the point where larval supply did not limit recruitment (Theuerkauf et al. 2021). If so, the importance of water residence time on recruitment might be diminished or overshadowed by sufficient larval supply.

Our finding that wind fetch and elevation exert control over oyster recruitment in coastal Virginia could apply to intertidal oyster populations elsewhere. Habitat restoration and management in the context of these geophysical factors may enhance oyster reef reestablishment efforts. For example, optimizing reef placement in areas with lower wave energy supports oyster recruitment (Pamlico Sound, North Carolina;Theuerkauf et al. 2017), while elevating reefs in sediment-prone areas might mitigate the impact of sedimentation on larvae, improving restoration outcomes (Ipswich River Basin, Massachusetts; Baillie & Grabowski 2019). Incorporating these site-specific geophysical considerations into restoration planning may help ensure more consistent recruitment success across diverse coastal systems.

Vol.: (0123456789)

Although mature restored reefs in our study system match the abundance and size of adult oysters (> 25 mm shell height) on natural reference reefs (Smith et al. 2022a), we found higher recruitment on reference reefs than on restored reefs regardless of their age. This discrepancy could be attributed to differences in larval supply, processes affecting settlement success such as substrate conditions, or post-settlement processes that influence oyster survival (e.g., competition, predation; Grabowski et al. 2005;Lipcius et al. 2008;Hubbard & Reidenbach 2015). Larval supply can vary over relatively small spatial scales (km or less) based on distance to the source of larvae and the hydrodynamic conditions that transport and retain larvae (Lipcius et al. 2008;Theuerkauf et al. 2021). In theory, reference reefs might have greater recruitment if they are in close proximity to strong sources of larvae, while restored reefs may have diminished recruitment if they rely on distant larval sources. However, many restored reefs in our study system were built close to reference reefs (< 2 km), and substantial differences in recruitment persisted when analyzing the smaller subset of spatially paired reefs (comprising 20 distinct restored reefs each paired with a reference reef). Moreover, larval transport usually varies over broader scales than the distances between paired restored and reference reefs (Haase et al. 2012).

If larval supply is similar among reefs differing in restoration history, substrate conditions on restored reefs may be less favorable for oysters to settle and attach than reference reefs, potentially due to factors such as substrate stability and availability. However, the restored reefs in our system are almost entirely built from natural substrates (mainly oyster shells, but also some whelk and clam shells), as opposed to 'novel' substrates such as concrete or rock. Still, it is possible that the rugosity of restored reefs in our system at scales of meters to hundreds of meters does not yet provide the structural complexity required for consistently high settlement, as reef habitat complexity is expected to increase with reef age (Bahr & Lanier 1984, Grabowski et al. 2005;Temmink et al. 2021). Additionally, we can consider the high post-settlement mortality of oyster recruits (Knights et al. 2012;Baillie & Grabowski 2019): it is plausible that settlement to restored reefs is similar to reference reefs, but mortality on restored reefs is high and occurs quickly after settlement (when oysters are < 1 mm in shell height). Post-settlement mortality could be due to increased sedimentation or sediment re-suspension on younger, flatter, and more open reefs (since older, porous and more complex reef structures trap more sediment; Lenihan 1999). Restored and reference reefs might also differ in the density of small predators or a lack of refuge areas to protect young oysters from predators (because older reference reefs might provide more interstitial spaces or vertical relief to shield small oysters from predation ;Grabowski 2004). This hypothesis aligns with our observation that Xanthid mud crabs, a common predator of small oysters (Kulp et al. 2011), were twice as abundant on older reference reefs compared to developing restored reefs (see Appendix S1: Fig. S10) with low oyster recruitment. Despite the higher predator abundance on reference reefs, the presence additional structural complexity could contribute to reduced predation and lower mortality among oyster recruits on these reefs. Continued long-term data collection will help determine whether restored reefs eventually match the recruitment found on reference reefs, if local conditions lead to a persistent recruitment deficit on restored reefs, or if there will be population-scale consequences of suppressed recruitment for restored reefs (e.g., reduced resilience to environdisturbances, disease, or overharvest).

The refined habitat suitability layer developed in this study, which is based on modeled estimates of how wind fetch and elevation support above-average recruitment, offers insights into target locations and restoration approaches to enhance oyster populations in our study system (Fig. 7a). In areas with overlap between the two suitability layers-indicating high predicted recruitment (this study) and high adult biomass (Hogan & Reidenbach 2019;Smith et al. 2022b)-suitable substrate may limit oyster reef coverage (Fig. 7b-c). These sites may be good candidates for creating new reefs through shell hash deployment, novel construction materials (e.g., alternative substrates), or simply allowing existing reefs to naturally expand over time. In places where predicted suitable habitat for adult oyster reefs is high but recruitment suitability is low (Fig. 7d-e), it is possible that such areas experience lower settlement rates or higher post-settlement mortality. But it is important to note that the predicted 'high' recruitment density in our study, represented by 2.2-5.6 oysters per cm 2 , which encompasses both the mean and upper quartile, is notably higher than the maximum recruitment densities reported in other studies, which typically range from 0.2 to 1.0 oysters per cm 2 (e.g., Knights et al. 2012;Whitman & Reidenbach 2012;Baillie & Grabowski 2019;Esquivel-Muelbert et al. 2022). This suggests that even our 'low' recruitment densities, represented by 0.1-2.1 cm 2 , which encompass the minimum up to the mean, may still represent relatively high recruitment compared to these benchmarks, highlighting the particularly favorable conditions in our study area. However, these results should be interpreted within the context of our specific environmental conditions and may differ in other systems, where baseline recruitment densities vary.

To address the challenge of restoring oyster reefs in regions or under geophysical conditions associated with low recruitment, one restoration approach could involve transplanting spat-on-shell cultch (young oysters attached to shell) to artificially enhance recruitment (Dinnel et al. 2009). Additionally, because our study system has areas suitable for above-average recruitment near areas with lower recruitment, opting for transplantation of wild spat-on-shell, as opposed to introducing hatchery-produced alternatives, would offer cost reductions while maintaining local genetic diversity. In recruitment-limited areas, this approach of transplanting early recruits or juveniles has been successful in restoring habitat-forming taxa with a propagule stage (i.e., seeds, spores, larvae; Vanderklift et al. 2020) such as Olympia oysters (Ostrea lurida; Dinnel et al. 2009), seagrass (Zostera marina; Orth et al. 2006), and kelp (Fredriksen et al. 2020).

Conversely, efforts to improve adult habitat suitability may be helpful to restoring areas predicted to support high recruitment but low habitat suitability for adult oysters. These locations should support robust oyster recruitment provided there is sufficient suitable substrate for settlement, attachment, and growth. There are several locations behind the southerly barrier islands of the VCR that match this situation (Fig. 7f-g). In these and other places, it may be revealing to evaluate the survival of new recruits and identify the factors limiting adult oyster establishment prior to initiating a substantial restoration plan. Although water residence time and wind fetch cannot be manipulated to improve habitat suitability for adult oysters, intertidal elevation can be modified by relatively small-scale restoration projects (Lenihan et al. 2001;Schulte et al. 2009). In areas that are too deep, building significantly elevated oyster reefs may accomplish multiple goals by increasing suitable substrate for settlement and recruitment (Fig. 5c) and adult oysters (Smith et al. 2022b). Analogous approaches have been taken to improve the restoration of submerged salt marshes through the placement of thin layers of sediment (VanZomeren & Piercy 2020). However, the sustainability of such restorations depends on the longterm stability of the new substrate. In our system, restorations behind geologically stable islands (e.g., Smith Island and Wreck Island; Fig. 7) may be more ideal for the creation of back-barrier reefs when compared with islands experiencing rapid erosion or sediment overwash (e.g., Cobb Island and Cedar Island; Fig. 7). In coral reef and kelp forest restoration, substrate stability is crucial for the successful recruitment of larvae and spores, as unstable substrates increase the risk of detachment and mortality, and ultimately failed restoration efforts (Burek et al. 2018;Ceccarelli et al. 2020). Decisions to restore oysters in areas predicted to support high recruitment but lacking suitable adult oyster habitat will likely involve trial and error given the dynamic geomorphology in many barrier-island systems like those in coastal Virginia (Robbins et al. 2022).

Lastly, the economic viability of restoring reefs in potentially transient, distant locations would depend on the outcome of recruit survival and local adaptation (Bible & Sanford 2016;Hämmerli & Reusch 2002). Limited gene flow between distant populations may reduce genetic diversity over time, potentially diminishing adaptability and population fitness (Bible & Sandford 2016). Prioritizing areas with both predicted above-average recruitment and existing reef locations or focusing on sites with predicted suitable oyster habitat despite the lack of observed recruitment could yield more effective resource allocation. Furthermore, a metapopulation approach to planning reef locations may be beneficial in challenging areas for isolated populations to thrive (Lipcius et al. 2008). Building reefs in suboptimal locations could enhance connectivity, facilitating gene flow and ultimately contributing to the long-term sustainability of oyster restoration efforts in these areas.

In summary, our findings demonstrate the value of integrating standardized field studies with long-term monitoring data to reveal how environmental conditions shape recruitment of an imperiled foundation species across a highly heterogeneous landscape. By combining geophysical data from models and remotely sensed data products, we offer a comprehensive perspective on regional restoration planning that incorporates multiple life history stages. Furthermore, the validation and refinement of habitat suitability models underscores the broader importance of improving spatial planning for restoration, which should ultimately enhance the success of population establishment, growth, and resilience (Puckett et al. 2018;Smith et al. 2022b). Our study advances knowledge for oyster restoration science and planning, and provides a framework that can be extended to other systems and species in which recruitment constrains restoration success.