Surveys were conducted between 2013 and 2019 across the U.S. Pacific as part of the NOAA Coral Reef Conservation Program’s National Coral Reef Monitoring Program (NCRMP) ( Figure 1A ; Table S2 ). Eight regions were surveyed: main Hawaiian Islands (MHI), Northwestern Hawaiian Islands (NWHI), northern Mariana Islands (NMI), southern Mariana Islands (SMI), Wake Island, Line Islands, Phoenix Islands, and American Samoa (SAMOA). Guam and CNMI were divided into the NMI and SMI based on the large variation in environmental conditions and heat stress exposure. The Line Islands, Phoenix Islands and Wake are part of the Pacific Remote Islands Marine National Monument but were separated into different regions for this study given their geographic distance and range of exposure to heat stress. Each region was surveyed once every three years, which resulted in two to three time points per region ( Table S2 ).

A total of 1,405 forereef sites across 34 islands, atolls, or reef systems, henceforth referred to as “islands”, were included in this study ( Table S2 ). Sites within each region, were selected using a stratified random sampling design; they were randomly selected within hard bottom forereef habitats stratified into three depth strata (shallow: 0–6 m, mid: > 6–18 m, and deep: > 18–30 m). Sample size within each stratum was proportional to the amount of hardbottom area and the variance of coral density in previous years (Smith et al., 2011) such that more sites were sampled in strata with higher variance and more hardbottom. We subdivided large islands (which consist of spatially varying oceanographic exposure, reef structure, local human population density and management area) into smaller units called “sectors” (Heenan et al., 2017). Sub-island sectors only exist for larger populated islands such as Tutuila, Guam, and the main Hawaiian Islands.

At each site, one 30-m transect line was deployed along the depth contour. Juvenile scleractinian colonies (1 – < 5 cm maximum diameter) were enumerated visually in situ as part of comprehensive demographic survey (Winston et al., 2020). Juveniles were counted within one to three 1 x 1 m segments along the transect (at 0–1 m, 5–6 m and 10–11 m). Three segments were surveyed at 93% of the sites, with the remaining sites having fewer segments due to high colony density, depth or logistical constraints. Juveniles were distinguished from tissue remnants by their growing edges and mostly circular nature. Segments were used to capture the heterogeneity in colony distribution, but were not considered spatially independent replicates, therefore colonies were pooled across segments to calculate total juvenile density (colonies m^-2) at each site. To assess benthic cover, digital images (3-10 MB and 0.7 to 1 m^-2 per image) were taken with a Canon PowerShot SD1200IS (2013–2014), Canon PowerShot S110 (2015–2017), and Canon PowerShot G9x (2018–2019). Thirty images were collected along the same transect as the juvenile surveys at 1 m increments 1 m above the benthos. Benthic cover data were extracted from images using Coral Point Count with extensions (Kohler and Gill, 2006) (2013–2014) and the web-based image annotation tool, CoralNet (Beijbom et al., 2015) (2015–2019). For each image, the organism or type of substrate was identified beneath each of ten randomly overlaid points. Points were pooled across all imagery to generate site-level percent cover of key benthic functional groups (300 points/site) (Lozada-Misa et al., 2017). Functional groups included hard coral, CCA, turf algae, macroalgae (including upright and encrusting macroalgae), and unconsolidated substrate (e.g. sand and rubble).

In this study, we conducted three sets of statistical analyses with total juvenile density as the response variable– a spatial analysis of regional patterns in juvenile density, a temporal analysis of changes in density across years, and a correlative analysis that identified relationships between juvenile density and ecological/environmental predictors. Due to the complex design across depth strata and unbalanced nature of our study design, we have employed the ‘survey’ package in R which allows us to generate unbiased estimates when using site-level data (Lumley, 2022). For the regional spatial analysis, 1,405 sites (only including depth strata that had at least 2 sites) were first inverse proportion weighted based on their selection probability within strata for each survey year (number of possible sites in a stratum/number of sites surveyed). The nested structure was defined as Survey Year, Region, Island, Sector, and Depth Bin using the ‘strata’ argument of the svydesign function. To test for statistical differences between regions, survey-weighted generalized linear models (GLMs) were fit using the svyglm function followed by post hoc tests (svytest), with juvenile abundance as a Poisson response variable, log(survey area) as an offset to account for variable survey areas, and region as a fixed effect. For the spatial analysis, ‘Survey Year’ was included in the nesting structure to account for variable sampling across years, but was not included as a fixed effect in the model. To visualize regional patterns in juvenile density, weighted means and standard errors were calculated at the regional level by averaging weighted site-level data across all years for a given region using the svymean function in the ‘survey’ package in R. We also calculated island-level weighted means to describe general patterns but did not conduct statistical analyses at the island-level. Year was not included as a fixed effect in the regional analysis because some regions and strata had inadequate temporal sampling. For the temporal analysis, only strata that had at least 2 sites and were surveyed in all years were included in the analysis (total of 1010 sites) to facilitate consistent trend detection. NWHI was not included in the temporal analysis because there were not enough strata that were sampled consistently across the survey years. The nesting structure was consistent with the regional analysis. To identify temporal differences in density within each region, we used the svyglm function to fit separate models for each region using juvenile abundance as a Poisson response variable, log (survey area) as an offset, and survey year continuous as a fixed effect. We then conducted post hoc tests (svytest) with a Benjamini-Hochberg multiple test correction to limit false discovery rate while preserving statistical power (Waite and Campbell, 2006). Temporal patterns were visualized using regionally weighted means and standard errors for each year.

To identify predictors of juvenile density, we summarized a range of ecological, habitat and environmental factors known to influence coral recruitment, post-settlement survival and coral communities ( Table 1 ). A total of 95 sites distributed across the 8 regions were dropped due to missing predictor data, resulting in a total of 1,310 sites. Site and sector-level percent cover area of scleractinian corals was included given its ability to affect juvenile abundance through mechanisms of larval supply and substrate availability for larval settlement. Island size is not an ideal predictor because regions such as the NWHI have very little emergent land, but often have considerable submerged reef area. Hence, we chose to represent island size/reef area as ‘sector coral cover area’. This metric is defined as mean cover (%) of all hard coral at the sector-scale multiplied by the hard bottom habitat area for a given sector. Habitat features such as depth, percent cover of CCA, macroalgae, and unconsolidated substrate were included due to their known influence on settlement preference and survival of early life stages. Biomass of herbivorous reef fishes was averaged to the stratum-level from stationary point count surveys in the same strata and years that juvenile surveys were conducted (for detailed methods see Heenan et al., 2017). NCRMP fish surveys were not co-located with juvenile coral surveys. Survey effort varied across strata due to the stratified random sampling design and contained at least 3 sites/stratum. Heat stress was included and accounted for using two methods. The first, heat stress severity (HS_sev), was quantified using the mean of NOAA’s Coral Reef Watch annual maximum degree heating weeks (DHW) (Liu et al., 2013) for years when DHW was greater than 0 °C-weeks in the ten years prior to the date of each survey. Years since a bleaching-level heat stress event (>4 °C-weeks; HS_ts) was also included. HS_ts at sites which had never experienced a heat stress event (245 sites), HS_ts was set to 32 years, the maximum of the climatological record, to retain them in the analysis. Gridded daily sea surface temperature (SST) data were derived from the NOAA Daily Global 5km Geo-Polar Blended Sea Surface Temperature Analysis (v1.0; ~ 0.05°). Range in SST is a commonly used metric to quantify short-term thermal variability on coral reefs (Sully et al., 2019). Here, we quantifed SST_var as the mean biweekly range in daily nighttime SST at each site in the ten years preceding the survey. Other satellite-derived environmental variables such as 8-day surface chlorophyll-a (Chla; Hollmann et al., 2013; ESA OC CCI v5.0; ~ 0.05°), a proxy of phytoplankton biomass, primary productivity, and turbidity, were included given their ability to affect juvenile growth and survival. Satellite-based Chl_a values were masked at its original resolution (~ 0.05°) prior to the analysis as a quality control measure. Chl_a rasters were first overlaid with a bathymetry raster at 15 arc seconds (Tozer et al., 2019). Chl_a pixels with an area greater than 5% characterized in depth of 30 m or shallower were then masked and substituted by using Chl_a values from the nearest unmasked pixel (Winston et al., 2022). These spatially and temporally resolved satellite data were statistically summarized over the ten years prior to every individual georeferenced and time-stamped survey record. The relative difference in wave energy between sites was calculated using wave power which is a combination of wave height and period derived from NOAA’s WaveWatch III global, full spectral wave model over the entire 1979–2012 climatology (Aston et al., 2019; Winston et al., 2022). Wave power was included due to its influence on larval settlement and survival. Gridded human density at 2.5 minutes (~ 5 km; 2000–2015; NASA GPW v4.11; doi.org/10.7927/H49C6VHW) was included as a proxy of human disturbance and its ability to affect coral communities. Latitude (absolute) was included as a proxy for factors that impact corals across our broad spatial gradient but were not explicitly measured here. Interaction terms were also included to examine the influence of each predictor on juvenile density under variable HS_sev. The interactions of SST_var x HS_sev and Mean Chl_a x HS_sev were not included in the full model because Mean SST and Mean Chl_a represented only a very narrow range of values at heat stress > 5°C-weeks, making interpretation difficult at high levels of heat stress.

We tested for potential collinearity between predictors using pairwise Pearson correlation coefficient (r), conservatively discarding variables when |r| > 0.65 (Dormann et al., 2013; Sully et al., 2019). Turf algae was strongly negatively correlated with CCA cover and was subsequently dropped. Likewise, latitude exhibited a strong correlation with Mean SST and was subsequently dropped. SST_var was also strongly correlated with HS_sev and subsequently dropped to allow us to test the role that heat stress severity has on juvenile density. The absence of collinearity was confirmed with variance inflation factors (VIF), which were lower than the a priori threshold of 3 for the subset of remaining predictors (Zuur et al., 2010). We used svyglm models to account for our sampling design using inverse probability weights as described above for the regional spatial analysis. To identify possible non-linear relationships between juvenile density and the predictors, we first visualized the data graphically, then used svyglm to model each parameter individually as a linear, 2^nd, and 3^rd order polynomial term. We used Wald’s tests and Akaike Information Criterion (AIC) values to identify orders that led to the best model performance. Depth was best modeled as a 2^nd order polynomial, and site-level coral cover was best modeled as a 3^rd order polynomial; all other variables were fit linearly. We scaled and centered all predictors to a mean of 0 and a standard deviation of 1 to facilitate meaningful comparisons of parameter estimates. The full model was fit with the predictors in Table 1 . To select the best predictive model, we used backwards model selection with Wald’s tests to sequentially drop non-significant variables at alpha = 0.05, then confirmed relative model performance using AIC. Validation of the final model revealed the presence of overdispersion (phi = 13.1), which was accounted for in the estimation of parameter standard errors. We observed no obvious structure in the standardized residuals. Goodness of fit was calculated using McFadden’s pseudo-R². We visualized predicted relationships by using the best fit model to predict juvenile coral abundance across the observed range of each variable, setting the survey area constant at 1 m² and all other predictors to their observed means. To visualize interactions between two continuous variables, one predictor was converted into a binned categorical term, and the predicted relationships between the response and other continuous predictor were plotted within each bin. All data were analyzed in R version 4.1.0 (R Core Team, 2019).

Regional mean (± SE) juvenile density averaged across years ranged from 1.51 ± 0.23 to 10.57 ± 0.70 colonies m^-2 and varied significantly across regions ( Figure 1B ). The highest regional densities were observed in the north Mariana Islands, Northwestern Hawaiian Islands, and south Mariana Islands. The lowest regional densities were observed at Wake Island (1.51 ± 0.23 colonies/m) and the Phoenix Islands (2.14 ± 0.52 colonies m^-2). Island-level weighted means (± SE) are reported in Figure S2 . Juvenile colony density varied significantly between survey years, but this change was not uniform across all regions. Post hoc tests revealed significant temporal changes in the Main Hawaiian, Line, and north Mariana Islands ( Figure 2 ). In the main Hawaiian Islands, juvenile densities increased significantly by 55% between 2016 and 2019. In the Line Islands, they increased by 48% between 2015 and 2018. In north Mariana Islands, juveniles also increased significantly—by 62% between 2014 and 2017. The other five regions experienced no significant changes in juvenile density through time, and no region experienced significant declines ( Figure 2 ).

Juvenile density was best predicted by a combination of (listed in order of parameter estimate magnitude) the interaction of HS_ts and HS_sev, depth, site-level coral cover, HS_ts alone, sector coral cover area, human density, unconsolidated substrate cover, macroalgal cover, herbivore biomass, the interaction of sector coral cover area and HS_sev, and the interaction of wave power and HS_sev sector coral cover area, and ( Figure 3 ; Table S3 ). This best-fit model accounted for 41% of the variability in juvenile density. Density showed a quadratic relationship with depth, with the model predicting a peak of 14 colonies m^-2 at 19 m depth, then a decrease in density down to depths of 27 m ( Figure 4A ). Site-level coral cover also demonstrated a non-linear relationship, with predicted juvenile density peaking at 12 colonies m^-2 at 15% coral cover then decreasing precipitously on reefs with greater than 20% coral cover ( Figure 4B ). Juveniles increased with the number of years since a 4 °C-week heat stress event ( Figure 4C ) and sector coral cover area ( Figure 4D ). While the increasing error with sector-scale coral cover is the result of one island in the NWHI with high coral area values —French Frigate Shoals, this predictor is still significant even after removing French Frigate Shoals sites (n = 30,results not shown). Density decreased significantly with increasing human density ( Figure 4E ). Juvenile corals also significantly declined with increasing percent cover of unconsolidated substrate ( Figure 4F ), percent cover of macroalgae ( Figure 4G ), and herbivore biomass ( Figure 4H ).

While HS_sev and wave power did not significantly predict juvenile density alone ( Figures 4H, I ), the interaction of HS_sev and HS_ts, and to a lesser degree the interaction of HS_sev and wave power, did significantly predict juvenile density ( Figure 3 ; Table S3 ). On reefs that experienced heat stress within three years prior to the survey, predicted juvenile density decreased from 12 to 3 colonies m^-2 with increasing severity of heat stress ( Figure 5A ). However, if the most recent heat stress event had occurred between 3 and 10 years prior to the survey, predicted density decreased only slightly with increasing heat stress ( Figure 5A ). If bleaching level heat stress occurred more than 10 years ago, predicted density increased from 7 to 13 colonies m^-2 at low heat stress ( Figure 5A ), but interpretation of predictions above 2°C-weeks is difficult given the high uncertainty. The interaction of sector coral cover area and HS_sev significantly predicted juvenile density ( Figure 3 ; Table S3 ). At low HS_sev (≤ 4°C-weeks), predicted juvenile density increased with increasing sector coral cover area ( Figure 5B ), but above 4 °C-weeks, juvenile density did not change with increasing sector-level cover ( Figure 5B ). At low HS_sev (≤ 4°C-weeks), predicted juvenile density increased slightly with increasing wave power ( Figure 5C ), but above 4 °C-weeks, juvenile density decreased from an estimated 9.5 colonies m^-2 to 4.8 colonies m^-2 ( Figure 5C ) as wave power increased.

As found in previous studies ( Table S1 ), juvenile density was highly spatially heterogeneous. There was a seven-fold difference in mean juvenile density across the regions (1.51 ± 0.23 to 10.57 ± 0.70 colonies m^-2), yet was within the range reported at other regional scales such as the Great Barrier Reef (6.9 ± 0.3 juveniles m^-2; Trapon et al., 2013) and Micronesia (~ 4 juveniles m^-2; Mumby et al., 2013). Our mean values are lower than those reported in an Indo-Pacific meta-analysis of juvenile density which reported a mean of 17.1 ± 2.08 colonies m^-2 across 76 sites between 1967 and 2012 (Edmunds et al., 2015). This difference could be attributed to a combination of factors. First, unlike the studies included in the meta-analysis, we excluded colonies less than 1 cm leading us to underestimate total abundance < 5cm. However, Trapon et al. (2013) noted that colonies < 1cm typically comprise a small portion of the total reported abundance due to difficulty detecting these colonies with the naked eye. Second, given the high spatial heterogeneity of juveniles, the discrepancy may be simply explained by the minimal spatial overlap between studies in Edmunds et al. (2015) and the current study. Lastly, given the link between juvenile density and heat stress ( Figure 5A ), the disparity in the temporal span of these two analyses may explain the differences in density.

In a recent global meta-analysis of coral recruitment, Price et al. (2019) and Edmunds (2023) found a steady decline in recruitment on Indo-Pacific tropical reefs (< 20° latitude) since 1974 and an increasing trend in coral recruitment on subtropical reefs (> 20° latitude). They suggest that the latitudinal shifts in recruitment are evidence that corals may be responding to rising sea surface temperature in the tropics by shifting their latitudinal range poleward. This recruitment study may help explain both the lower juvenile density in tropical regions such as the Phoenix Islands and American Samoa and the higher densities in the subtropical reefs of the NWHI, MHI and NMI. The especially low juvenile density observed at Wake Island (1.41 ± 0.20 colonies m^-2), does not follow the latitudinal patterns. At 19.3°N, Wake is situated at the boundary of the tropical and subtropical zones, but its low density may be explained by its high geographic isolation (> 1,500 km to closest shallow water reef) and its long history of human disturbance and military occupation. While our results demonstrate that high juvenile density can be found in certain regions, despite the recent mass bleaching, poleward expansion of corals is concerning given increasing heat stress and barriers to growth at high latitudes (Couch et al., 2017; Edmunds, 2023)

Depth plays a major role in structuring coral communities and the current study highlights the complex effect of depth. Our results suggests that Indo-Pacific studies should expand beyond the first 10m of depth that juveniles are typically surveyed at ( Table S1 ) to gain a more comprehensive understanding of juvenile dynamics. Previous studies have found a negative (Rogers et al., 1984; Carlon, 2001; Cardini et al., 2012), unimodal (Birkeland et al., 1981; Rogers et al., 1984), or no effect of depth (Edmunds, 2004). There are several possible mechanisms that may explain these non-linear effects of depth. First, shallow water corals exist in a more environmentally variable and stressful habitat. Corals are routinely exposed to higher heat stress, land-based inputs, and wave energy, all of which can be bottlenecks for survival (Edmunds, 2004; Fabricius, 2005; Reidenbach et al., 2009). The lower juvenile density deeper than 20 m is commonly observed on reefs (Birkeland et al., 1981; Rogers et al., 1984; Carlon, 2001; Cardini et al., 2012) and may also suggest that reduced light availability on deeper reefs is suboptimal for larvae settlement. In fact Maida et al. (1994) found that coral spat preferentially settled on the underside of settlement tiles along the edges. This suggests that larvae often prefer cryptic habitats, but require a certain level of light for settlement. Another hypothesis posed by Carlon (2001), is that deeper reefs have more competition for the cryptic habitats and higher abundance of filter feeding organisms that may impair coral settlement and growth. Thus intermediate depths provide the ideal habitat that reduces exposure to environmental stress, provides ideal light conditions and minimizes the bottlenecks to survival on deeper reefs.

Our study demonstrates that the relative importance of heat stress on juvenile corals depends on the severity and timing of heat stress. Heat stress can reduce fecundity by removing reproductively mature colonies and impairing gametogenesis (reviewed by Ritson-Williams et al., 2009; Edmunds, 2023), cause direct mortality of coral juveniles (Edmunds, 2004), and result in declines in juvenile density (Sheppard et al., 2017; Dajka et al., 2019; Koester et al., 2021; Huntington et al., 2022). Interestingly, we found that between 3 and 10 years following a heat stress event, the severity of heat stress had a much weaker effect on juvenile density. In the MHI, the 2019 surveys were conducted 3 years and 8 months after the 2015 severe bleaching event, suggesting that by July of 2019, the negative impacts of this event on juveniles may not be as apparent because new larvae had settled and survived past the coral spat stage. Across the Indo-Pacific, several studies have found evidence of increased juvenile density following mass coral bleaching. In Moorea, juvenile density rebounded to the pre-disturbance levels just four years after repeated disturbances (Adjeroud et al., 2018). In the Line Islands, Jarvis Island had extreme heat stress during the 2015/2016 event (30 °C-weeks) and while there were early signs of recovery in 2018, juvenile density remained 30% lower than 2015 values ( Figure S2 ). Alternatively, on reefs that experienced heat stress greater than 10 years prior to the survey, we found that juvenile density increased with heat stress. Although seemingly contradictory to the results discussed above, it is unlikely that small changes at low heat stress are meaningfully influencing juvenile density given the uncertainty around the prediction above 2 °C-weeks. While high recruitment following severe heat stress and bleaching is promising, our results suggest that potential stock-recruitment relationships disappear following elevated heat stress ( Figure 5C ). These results support findings by Hughes et al. (2019) who discovered that repeated mass bleaching events in the Great Barrier Reef have resulted in an unprecedented 89% decline in recruitment following the mortality of brood stock. Together, these results highlight the increasing threat of global climate to impair stock-recruitment relationships and further impede long-term recovery.

Across the Pacific, juvenile density varied non-linearly with site-scale coral cover, thus demonstrating the series of mechanisms driving coral demography. Juvenile density increased with site-level coral cover on reefs with < 15% coral cover. Underscoring this pattern, juvenile density also showed a simple positive association with sector-level coral cover area. Together these results suggest the importance of stock-recruitment relationships across spatial scales. Stock-recruitment relationships in corals has been identified elsewhere and indicates that larval supply and recruitment increases with the abundance of reproductively mature adult corals (Caley et al., 1996; Connell, 1997; Doropoulos et al., 2015; Hughes et al., 2019). Stock-recruitment relationships can be especially important in remote reefs that rely on self-recruitment and may explain why early recovery at Jarvis Island three years post-bleaching has been dominated by a handful of thermally tolerant taxa, while taxa locally extirpated during bleaching have yet to return (Huntington et al., 2022). Above 15% coral cover, juvenile density decreased precipitously with few juveniles predicted above 60% coral cover. These results support previous studies that found lower juvenile density on reefs with high coral cover (Hughes, 1985; Connell, 1997; Mumby, 1999) and suggests that adult corals monopolize the substrate needed for larval settlement. Furthermore, scleractinian corals can have strong allelopathic effects on early coral life stages, resulting in mortality and reduced growth rates of coral spat, demonstrating their ability to outcompete early life stages after settlement (Fearon and Cameron, 1997). Studying coral demography across this basin-wide scale provides the robust gradient in coral abundance that allows us to see the complex series of mechanisms driving juvenile abundance that are more difficult to identify in small-scale juvenile studies.

In the current study, we used human density as a proxy for human activity in coastal regions and found that juvenile corals in the U.S. Pacific decrease with local human density. These results are consist with a study in the Line Islands that have found lower coral recruitment on reefs near populated regions (Sandin et al., 2008). Human activity can negatively impact early coral life stages through a wide range of well-documented mechanisms. Coastal eutrophication from human activity increases stress to adult corals and enhances competition with algae that can lead to reduced fertilization and gamete formation, reduced larval settlement, and increased juvenile mortality (reviewed by Fabricius, 2005). Sedimentation can prevent settlement, reducing light availability causing corals to settle on exposed substrates thus increasing post settlement mortality due to smothering (reviewed by Fabricius, 2005). Future studies should address which aspects of human presence are negatively impacting juvenile density in this central and western Pacific region.

Percent cover of unconsolidated substrates such as sand and rubble are generally perceived as suboptimal habitats for early coral life stages. Our results support previous studies that have found a negative correlation between juvenile corals and unconsolidated substrate (Chong-Seng et al., 2014; Dajka et al., 2019). With the constant or periodic movement of the substrate in wave impacted areas, sand and rubble can impair post settlement survival by abrading or burying early life stages (Fox et al., 2003; Birrell et al., 2005; Kenyon et al., 2023). The role of unconsolidated substrates as a demographic bottleneck has the ability to impede recovery following disturbance. With the projected increase in the frequency and severity of heat stress and storm events (Cheal et al., 2017; Eakin et al., 2019) and the corresponding loss of live coral and increased bioerosion, reefs are projected to experience a loss of reef structure in the coming decades (Bozec et al., 2015). Continued deterioration of the reef structure could result in a feedback loop of impaired resilience to disturbance, whereby frequent disturbance prevents stabilization of the substrate, which increases mortality of early coral life stages, thus preventing coral populations from reaching reproductive maturity and re-stabilizing the substrate.

Our results confirm the deleterious, albeit weaker, role of macroalgae. Our study supports previous studies that have found a similar negative correlation between macroalgae cover and juvenile abundance (e.g. Edmunds and Carpenter, 2001; Mumby et al., 2013; Chong-Seng et al., 2014; Dajka et al., 2019). Macroalgae and turf algal mats are known to negatively impact early coral life stages by inhibiting coral settlement (Birrell et al., 2008; Vermeij et al., 2009; Doropoulos et al., 2014; Webster et al., 2015), impairing colony growth, and overgrowing coral juveniles thus increasing post settlement mortality (Box and Mumby, 2007). Our result is somewhat surprising given that unlike these previous studies, a majority of our sites had < 10% macroalgae cover, suggesting that even at low algal cover, macroalgae may be impacting juvenile corals in the Pacific. The low macroalgal cover also suggests that detection bias is less of a concern compared to regions that have high macroalgal cover.

Herbivores have a complex relationship with juvenile corals. While herbivorous fish are known to enhance coral settlement and post-settlement survival by controlling algae which outcompete coral spat and juveniles (e.g. Edmunds and Carpenter, 2001; Birrell et al., 2005; Box and Mumby, 2007), our study found that juvenile density declined with increasing herbivorous fish biomass. This negative relationship between abundance of juvenile corals and herbivores may be the result of corallivory. In previous studies, direct and indirect corallivory by parrotfishes has been associated with increased mortality and lower growth in early life stages of corals, highlighting complex tradeoffs between predation and algal competition (e.g. Penin et al., 2010; Doropoulos et al., 2016; Mumby et al., 2016). Unlike other Indo-Pacific reefs, the beneficial effects of herbivores may not have been as apparent given the low macroalgal cover across much of the U.S. Pacific. While a negative association between juvenile corals and herbivorous fish biomass may suggest that herbivores are detrimental to reef recovery, Mumby et al. (2016) demonstrated that the net impact of herbivorous fish during the first five years of life is positive for corals.

Wave power is a major environmental force that shapes coral reef communities, but relatively little is known about its role in early coral life stages. On reefs with low heat stress predicted juvenile density did not vary with increasing wave power, indicating that wave power may be buffering the effects of low heat stress. This hypothesis is supported by experimental studies which found that increased water flow has been associated with a lower bleaching response (Nakamura and van Woesik, 2001) and increased recovery from bleaching (Nakamura et al., 2003). Increased wave power reduces the boundary layer especially on small colonies, which enhances nutrient transfer and photosynthesis, thus improving survival during heat stress (Page et al., 2019). However, our study also indicates that above 4 °C-weeks, increasing wave power may reduce juvenile density, suggesting that the combination of physiological stress from heat exposure and shear stress compounds to impair early coral life stages. While a certain level of shear stress is needed to facilitate nutrient transfer, high shear stress can impede larval settlement (Reidenbach et al., 2009; Reidenbach et al., 2021) and result in elevated abrasion (Dollar, 1982). Although wave power was a fairly weak predictor, our results suggests that wave action may help buffer the impacts of low heat stress, but more exposed reefs may have a lower recovery potential following severe heat stress.

The abundance of early coral life stages can be highly temporally variable (Edmunds, 2023), and is further supported by the significant temporal changes in several regions ( Figure 2 ). In the MHI, we found that juveniles declined slightly between 2013 and 2016, but then almost doubled between 2016 and 2019 ( Figure 2 ). Similarly, juvenile density in the NMI and Line Islands increased 40 to 50% between surveys, respectively between 2014 and 2017. Over the same time periods, coral cover declined by 17 to 20% in all three of these regions (NCRMP unpubl. data), which has been attributed to a combination of the 2014-2017 mass bleaching in the MHI, NMI, and Jarvis Island (Line Islands) (Heron et al., 2016; Kramer et al., 2016; Vargas-Ángel et al., 2019) as well as crown-of-thorns starfish predation in the Line Islands (A. Pollock pers. comm.). Previous studies have suggested that juvenile corals more resilient to heat stress than adult colonies (reviewed by Edmunds, 2023), while others have reported declines in juvenile density immediately following severe prolonged bleaching (Sheppard et al., 2017; Dajka et al., 2019; Koester et al., 2021; Huntington et al., 2022). Our study highlights that despite the increasing heat stress, the negative effects of heat stress on juvenile abundance can be ephemeral with certain regions demonstrating promising recovery potential through increased juvenile density following mass bleaching. These results contradict much of the literature indicating that heat stress depresses gonad development, gamete viability and spawning for up to several years after the event (Ward et al., 2002; Mendes and Woodley, 2002; Hagedorn et al., 2016).