Study site and bleaching timeline

This study was conducted at Lizard Island, in the northern GBR, Australia (14°40′S, 145°28′E), across a 24-month sampling period, ranging from January 2016 to January 2018 (Supplementary Fig. 1). During this timeframe, the GBR experienced unprecedented back-to-back mass bleaching of scleractinian corals (i.e., February–April 2016, mean sea surface temperatures (SST) of 29.1, 29.1 and 27.8 °C, respectively; and from January–March 2017, mean SST of 28.8, 28.8 and 28.7 °C, respectively), as a result of prolonged elevated sea-surface temperatures12. Bleaching in 2016, severely impacted the northern 1000 km of the GBR, particularly the northern region between Port Douglas and the Torres Strait, while the 2017 bleaching event primarily impacted the reef further to the south, between Townsville and Cooktown12. We quantified changes in both coral cover and fish assemblages in response to consecutive mass bleaching across four sampling periods: (1) January–February 2016 (i.e., prior to mass bleaching), (2) April 2016 (i.e., during the peak of the first bleaching event), (3) October 2016 (i.e., 6-months after the peak of the first bleaching event) and in (4) January 2018 (i.e., 20-months after the peak of first bleaching event; ~10-months after the second bleaching event).

Sampling methodology

We employed a novel sampling methodology (Supplementary Notes 1) specifically designed to minimise diver effects, while yielding replicate, high-resolution counts of small, visually apparent coral-associated reef fishes (comprehensive methodological description provided in15). This method used a series of ‘photoquadrats’ of the same 1 m2 area across the entire 24-month sampling period, in which both coral cover and reef fishes were quantified simultaneously. Hence, each specific 1 m2 quadrat area was sampled four times. This method, therefore, provides a direct quantification of small-scale spatial overlap for both corals and fishes over time.

In total, we surveyed 19 locations across Lizard Island, predominately within the lagoon (Supplementary Fig. 3). At each site, divers swam a transect (range: 50–210 m; transect length dependent on individual reef length) along the reef ‘crest’ (at 0–4 m below chart datum; sites were chosen haphazardly) taking photographs of a 1 m2 quadrat (Camera: Nikon Coolpix AW130) at ~5 m intervals (range: 12–38 quadrats per indv. transect). Reef ‘crest’ habitats were chosen as they typically boast high coral cover37. Each replicate quadrat location consisted of three images: an undisturbed horizontal perspective photograph of the reef and coral-associated reef fishes (at a distance of ~2 m), taken within seconds of reaching the site and prior to the placement of the quadrat (i.e., reducing the so-called diver effect ; Supplementary Notes 1), ensuring all fishes 1.5 m above the substratum were included in the photograph; a second horizontal perspective photograph with the 1 m2 quadrat in place, using the identical camera placement as in the first image; and a planar perspective photograph (i.e., bird’s-eye view) of the 1 m2 quadrat in place over the substratum. Upmost care was taken when placing quadrats to not damage corals or other benthic organisms. Disturbance to the fish community by self-contained underwater breathing apparatus (SCUBA) divers was minimal and fish activity returned to normal soon after the photographs were taken. For subsequent sampling periods, transect starting points were identified using global positioning system (GPS) coordinates, while precise 1 m2 quadrat locations were identified using a second UW camera containing all previous images as a reference. All quadrat locations were photographed across four sampling periods, at 0, 2, 8 and 24-months. The duration of each sampling trip was ~2 weeks. All transects were conducted on SCUBA by two divers between 09:00 and 16:00 h (all photographs were taken by R.P.S. and S.B.T.). This study was observational in nature, no material was removed and all access was covered by permits granted by the GBR Marine Park Authority.

Image analysis

Per sampling period, a total of 132 (out of 451) 1 m2 photoquadrats were analysed. Since we were explicitly examining the response of coral-associated reef fishes to mass bleaching, only quadrats with a minimum live coral cover of 20% in the first sampling period were analysed. ‘Coral’ is used in the broad sense of the term and refers to taxa in the orders: Alcyonacea, Corallimorpharia, Helioporacea, Scleractinia; class: Hydrozoa (Millepora spp.)38. To quantify reef fishes, all selected photoquadrats were processed in Adobe Illustrator, by drawing an outline of the quadrat on the first photograph of the series (i.e., undisturbed), using the second photograph in the series as a reference. All visible reef fishes within the delineated 1 m2 (and ~1.5 m above the quadrat) were recorded to species level and categorised as either adult or recruit. Recruit was defined as the presence of juvenile colouration and/or <25% of adult maximum size, and therefore, this category included both newly settled reef fishes, as well as juvenile fishes a few months in age. Visually similar species Chromis viridis and Chromis atripectoralis were grouped into one species category C. viridis to avoid misidentification. Coral-associated damselfishes (Pomacentridae) were classified into two categories: obligate coral dwellers and facultative coral dwellers, in accordance to their documented live coral dependency, following3 (Supplementary Table 3).

Coral cover was quantified using the third photograph in the series (i.e., bird’s-eye view), using the software photoQuad39, which generated 40 randomly stratified points over each photoquadrat. For each point, the underlying benthic covering was recorded, i.e., coral to the lowest taxonomic level (generally genus or species), growth form, bleaching status, etc. Only 16 of 21,120 benthic points examined were categorised as unidentifiable, and were therefore excluded from analyses. The effect of random sampling, as well as variation in quadrat placement, was assessed by examining 40 random points from the first sampling period vs. the exact same 40 points from the third sampling period (October 2016, 6-months after peak bleaching; n = 15 randomly selected quadrats). The results showed just a 1.4% difference, and hence, our method appears to provide a good indication of benthic changes among temporal samples (analyses published in 15). For consistency, all images were processed by one person (S.B.T.).

Statistical analyses

We used GLMMs to test for differences in the proportional cover of total live corals, and Acropora spp., as well as the abundance of all fishes, coral-associated (facultative + obligate) damselfishes, facultative coral-dwelling damselfishes and obligate coral-dwelling damselfishes among the four sampling periods, i.e., before, during, 6-months after and 24-months after mass bleaching. Proportional coral cover data were examined using a GLMM with a binomial distribution and, where necessary, fitting an observation-level random effect to account for overdispersion. Fish abundance data were examined using a GLMM with a negative binomial distribution, to account for the non-normal and overdispersed nature of the count data. In all models, sampling period (before, during, 6-months after and 24-months after mass bleaching) was fitted as a fixed effect, while quadrat ID, nested within transect ID, were fitted as random effects, to account for the lack of spatial independence and the repeated measures sampling design. Model fits were evaluated using residual plots.

Data focussing only on the abundance of coral-associated damselfish recruits, and their association with coral cover, were also examined using GLMMs. Initially, the abundance of recruits in all coral-associated damselfish species and recruits in species with an obligate or facultative coral dependency were compared between summer sampling periods (January/February 2016 and January 2018), which aligns with summer recruitment pulses40. In this case, zero-inflated GLMMs with a negative binomial error distribution were used to account for the non-normal, overdispersed and zero-inflated nature of the data. In the two models, sampling period (January/February 2016 and January 2018) was fitted as a fixed effect, while quadrat ID nested within transect ID, was fitted as random effects to account for the lack of spatial independence and the repeated measures sampling design. Subsequently, relationships among the abundance of all coral-associated damselfish recruits and coral cover were explored separately for the January/February 2016 and January 2018 sampling trips. The relationships between recruits and coral cover were specifically explored separately between trips to assess if the nature of the relationship had changed. However due to this division, a Bonferroni correction was applied to subsequent models (α = 0.025). In these models, coral cover was considered as an explanatory variable in two ways: total coral cover and cover of corals preferred by damselfishes. The corals preferred by damselfishes was based on3 and included all Acropora spp., Echinopora lamellosa, Echinopora mammiformis, Porites cylindrica, Porites nigrescens, Pocillopora damicornis, Pocillopora verrucosa and Seriatopora hystrix. GLMMs were all based on a negative binomial error distribution to account for the non-normal and overdispersed nature of the data. Zero-inflated models were utilised to account for the large number of zeroes in the data. An extreme outlier was present in the January/February 2016 data and analysis was performed both with and without this data point. Where the outlier influenced model significance, this is noted in the results. Due to the nature of the Acropora spp. cover data, no formal analysis was conducted to examine the relationship between Acropora cover and recruit density. Furthermore, the relationship between facultative and obligate recruit density, and the cover of corals was visualised graphically, but again, the nature of the data prohibited formal statistical analysis. Statistical modelling was performed in the software R41, using the lme442 and glmmTMB43 packages.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.