The constraint of ignoring the subtidal water climatology in evaluating the changes of coralligenous reefs due to heating events

Predicting community-level responses to seawater warming is a pressing goal of global change ecologists. How far such predictions can be derived from a fine gradient of thermal environments needs to be explored, even if ignoring water climatology does not allow estimating subtidal marine heat waves. In this study insights about the influence of the thermal environment on the coralligenous community structure were gained by considering sites (Sardinia, Italy) at different temperature conditions. Heating events were measured (by loggers at 18 m, 23 m, 28 m, 33 m and 38 m deep) and proxies for their duration (the maximum duration of events warmer than the 90th percentile temperature), intensity (the median temperature) and variability (the number of daily ΔT larger than the mean daily ΔT, and the number of heating events larger in ΔT than the 90th percentile ΔT) were selected by GAM models. Reliable predictions of decrease in coralligenous richness of taxa/morphological groups, with relevant increment in turfs and encrusting coralline algae abundance at the expenses of bryozoans were made. Associations to the different types of heating descriptor have highlighted the aspect (intensity, duration or variability) of the heating events and the threshold for each of them responsible for the trajectories of change.

Scientific Reports | (2020) 10:17332 | https://doi.org/10.1038/s41598-020-74249-9 www.nature.com/scientificreports/ heating events, so to provide tools to draw the trajectories of change of this community due to future warming scenarios. The results may contribute to identifying heating descriptors relevant to any kind of subtidal habitat vulnerability, for which MHWs metrics cannot be applied for the lack of climatological multi-year data.

Results
Thermal regimes at the sites and heating descriptors. Important differences in the thermal regimes among sites depending on the depth were evidenced by graphical inspection (Fig. 2). AS was the site with the least differences in temperature among depths and the coolest maximum temperature at the shallower depth (24.36 °C at 18 m), while CP was the site with the largest variations in temperature (at the least down to 28 m). CC was the site with the largest differences in temperature among depths and TA was the site where water stratification was complete, as at the end of the summer temperature was homogeneous from the surface down 38 m deep ( Fig. 2 and Fig S1, supplementary 1). A more accurate description of such relevant heterogeneity in thermal regimes among depths and sites was entrusted to the four heating descriptors (D90, MED, NDD and F90, the explanatory variables) selected by the models (from a set of 26, Table 1). Therefore aspects of the heating events were quantified by each descriptor: (1) duration by D90, the maximum duration (in days) of events warmer than the 90th percentile temperature; (2) heat intensity by MED, the median temperature and (3) temperature variability by NDD, the number of daily ΔT larger than the mean daily ΔT, and by F90 the number of heating events larger in ΔT than the 90th percentile ΔT (Table 1 and Fig. 3). Consequently, AS is clearly ranked the coldest, while CC and TA are the hottest sites (with CP in between), for the consistency of MED patterns among depths. However, in terms of duration of heating events AS is the one with the longest D90 consistently across depths, while at CC D90 is quite variable and increases with depth; at TA and CP large variations among depths were estimated, and not following a depth gradient. Moreover, NDD and F90 have measured different aspects of temperature variability: AS had the highest NDD and the lowest F90, TA had increasing F90 with depth, but homogeneous and low NDD, and CC had a clear pattern only for F90 which decreased with depth.
Furthermore, based on the SST climatology of the four sites, May-October 2019 SST data (Fig. 2) were used to detect and characterize the MHWs during the study 15,48 : three MHWs were found in AS and CC, four in CP and two in TA. MHWs differed in intensity and duration: some strong waves have occurred at all the sites, but they were not directly traceable in the subtidal temperature (Table S1 and  Coralligenous structure-heating descriptors models. A total of 42 taxa/morphological groups composing the coralligenous reef community were identified at several taxonomic levels (18 macroalgae, 7 anthozoans, 7 bryozoans, 5 sponges, 2 polychaete, 2 tunicates and 1 hydrozoan, Table 2) and used to calculate the community richness and evenness.  www.nature.com/scientificreports/ The models have highlighted the significant influence of the site on the community (in terms of taxa richness and evenness) and on all the categories of taxa considered, while depth only affected a few of the response variables: encrusting coralline algae, red algae, turf algae, sponges and bryozoans (Table 3 and Fig. S3 supplementary  1). However, the most interesting result was about the influence of the selected explanatory variables (D90, MED, F90 and NDD), alone or in interaction with site and depth on the coralligenous community structure (richness and evenness) and every single conspicuous category (supplementary 2). Surprisingly, all the coralligenous response variables were influenced by some of the heating descriptors (Table 3). In fact, both intensity and variability of heating events influenced the number of taxa of the community, as richness was negatively influenced by the MED larger than 18.5 °C, while it had a non-linear response to NDD (Figs. 4, 5 and Table 3. Conversely, the community evenness was only affected by MED intensity depending on the site, except for TA (Table 3 and Fig. 4).
Particularly, the duration of heating events (D90, Table S2 in supplementary 1) has influenced all the categories of taxa (except for the green algae), though interactively with depth or site (Fig. 6). D90 had a positive influence on the abundance of turf algae at all depths (with 10 being the most common threshold) and encrusting coralline algae only at 18 m, while deeper the influence had an opposite direction, with the same threshold (10). Furthermore, the site also changed the type of D90 influence on some categories (red algae were negatively influenced at CC and CP, but positively at AS, while sponges were positively influenced at both AS and CP and bryozoans negatively at TA), although exceeding 10 D90 always seemed to make a difference.
The 18.5 °C MED intensity of heating events seemed to be an important threshold as encrusting coralline algae and bryozoans received a positive and negatively effect, respectively, for temperature larger than 18.5 °C. Furthermore, the type of effect on the green algae depended on the site but 18.5 °C remained the threshold (Table 3 and Fig. 4). Moreover, variability in heating events did have an influence on some categories: turf and red algae were negatively affected by NDD larger than 28 and F90 larger than 4.5, while the same descriptors had a positive influence on bryozoans and sponges again for NND larger than 28 and F90 larger than 4.5 ( Fig. 5 and Table 3).

Discussion
In a temperate location like Sardinia, the dynamics of subtidal temperature changes relevantly among sites and depths during the water stratification period. The whole study refers to the heating events that have occurred in the water column independently of the surface MHWs, whose lag of influence in time and space on the deep water temperature has likely been driven by several hydrodynamic mechanisms 57  www.nature.com/scientificreports/ not allow accurate predictions on the thermal subtidal environment 58,59 . Therefore, having detected the MHWs during the study period was useful to characterize the surface conditions, although their direct effects on the deeper thermal environment still remain to be deeply explored. The objective limitation of not having long time-series of subtidal temperature data has forced to find proxies of heating events, resembling the MHWs metrics 60,61 . All the heating descriptors here selected, D90, MED, NDD and F90, influenced several coralligenous response variables either of the community or the main categories, so that they all might be considered useful predictors for climate change investigations on the coralligenous reef. Accordingly to the correlations among temperature descriptors highlighted during data exploration, the predicting value of those selected can be extended to the correlated descriptors (see below data analysis). The possibility of exploiting the use of such descriptors of heating events for other subtidal efforts should be investigated, in addition to MHWs metrics for the surface 15,48 . Currently, the response of other subtidal community/species to the thermal environment is being examined based on satellite-derived SST, if the systems investigated are very shallow 62 . Alternatively, studies for deeper systems have relied on in situ loggers measures where, for very specific locations, temperature data have started to be collected long before the study was really designed and where researchers could formulate the hypotheses only a posteriori [60][61][62][63][64][65] . However, estimating subtidal temperature anomalies based on the subtidal climatology by long-time series of field temperature data has to be supported, although it will concern a limited number of localities: the recent promotion of research projects or temperature monitoring programs (i.e. TMED for the Mediterranean, http: //t-mednet) that have extensively implemented logger data acquisition will be fundamental to gain information in wide areas and will enable detecting the role of sea water warming on relevant ecological patterns (i.e. mass mortalities) by correlative analyses 66 . www.nature.com/scientificreports/ Overall, the arbitrary choices made to define the D90, MED, NDD and F90 descriptors were striking, as far as there have been many significant associations with the coralligenous response variables selected. Among Table 2. Presence (+) or absence (−) of coralligenous taxa/morphological groups at all the sites (AS, CP, TA and CC) and depth (18, 23, 28, 33 and 38 m). Taxa categories (Cat.) of GAM models are indicated by the letters B bryozoans, C encrusting coralline algae, G green algae, S sponges, R red algae and T turf algae. Taxa/morphological  groups  Cat  AS  CP  TA  CC   18  23  28  33  38  18  23  28  33  38  18  23  28  33  38  18  23  28  33  38 Articulated Coralline Algae   www.nature.com/scientificreports/ the whole portfolio of the coralligenous taxa/morphological groups initially considered (before pooling data) in the analyses, none could be considered due to the negligible sample size for most of them and results could be obtained only from models run on conspicuous groups of taxa (obtained by data pooling). This represents a relevant bias of the methodology, especially for this habitat since the coralligenous community is importantly structured by a multitude of taxa with a natural very low area of occupancy 67 and this approach does not allow deriving any prediction for uncospicuous taxa. This aspect involves even the cnidarian species, iconic taxa to the habitat, that have been repeatedly affected in the Mediterranean by relevant mortality events during summer-autumn heat waves (i.e. in 1985, 1999, 2003 and 2008) 36,37,46,68 and whose disappearance may cause shifts in the community composition, favoring filamentous algae 69 . Change in the structure of coralligenous reefs with loss of species (lower richness and diversity) is a wide described phenomenon 61 . Although for the study sites there is lack of historical complete quantitative data, information about gorgonians mortalities gained for TA site 36,37,40 and CC [FP, personal observation], induce inferring that even the structure of the communities currently under focus is the product of pressing climate alteration. Nevertheless, the approach used in this study can only lead formulating predictions on large categories of taxa. In general, the expectation about the coralligenous community change to heating events, based on the overall vulnerability to warming treatments 38,39,45,60 , concerns the reduction of sea fan and encrusting sponge species abundance and the consequent increase of algae such as turfs or fast growing weeds 61 . In this study, although with the limit of the taxonomic resolution, some unequivocal predictions of the effects of heating events on the coralligenous community could be made, as there was a consistence response across sites and depths of associations (descriptor-coralligenous): these include the high predictive power of the median temperature for the taxa richness, for the encrusting coralline algae and bryozoans. Based on these, increases in temperature intensity would drive to a lower number of taxa and, among the categories composing the community, to an increase and decrease of encrusting coralline algae and bryozoans, respectively. This prediction is in accordance with the recent investigations on the whole community 61 and on bryozoans 70 , although for these latter temperature tolerance seems species-dependent 44 and varies among areas 71 . Conversely, experimental evidence on the response of encrusting coralline algae to heat is badly needed 72 . Moreover, in the projection of longer lasting heating events, turf algae will increase the spread, as they were positively associated to D90, evidencing that the duration longer than 10 days is likely the key aspect of heating events responsible for the increase in turfs on the coralligenous reefs 61 .

Alcyonum coralloides
Measures of heating variability, NDD and F90, have provided some further important insights, as they were good predictors for the decrease in turf (NDD over 28 daily temperature shifts larger than the mean daily shift) and red algae (F90 over 4.5 heating events larger than the 90th percentile), and increase in sponges (F90 larger than 4.5) and bryozoans (NDD larger than 28). This result may be seen in apparent contrast with expectations, but the most likely interpretation is that large values of these two descriptors correspond to breaks to the trajectories of change driven by increases in intensity and duration of heating events, stressing the need of not neglecting heat variability in studies of community responses to sea water warming. At this regards it seems important to highlight that NDD and F90 do not estimate the frequency of heating events, an aspect of heating events that was not taken into consideration in this study.
Certainly, important thresholds for whatever type of effect (negative or positive) on the coralligenous community were found for all the temperature descriptors: 18.5 °C for the median temperature, 28 for NDD, 4.5 for F90 and 10 for D90. Although related to the specific study period that does not include the cold months, all these thresholds together provide evidence for formulating predictions on the coralligenous reef community on the basis of the heating events. In fact, independently of their absolute values, they seem edges highly consistent across the response variables of the community investigated, suggesting that specific temperature intensity, variability and duration of heating events affect the coralligenous community, both as a whole and its conspicuous taxa.
Overall, discriminating the effect and quantifying the specific magnitude of influence of the heating events, whether the duration, intensity or variability represents one of the main results of the approach used. The other Table 3. Summary of the significant GAMs results showing only the significant explanatory selected variables: site, depth and the temperature descriptors for the heating duration, temperature intensity and variability. N = 400. Dev = the deviance explained by the models. "X" = interaction. www.nature.com/scientificreports/ output regards the context-dependence of several predictions, either on the community or the categories of taxa, as the type of associations (negative or positive) changed among sites or depths. However, the categories of taxa are composed of numbered species whose contribution in abundance changes depending on the sitexdepth context and, because vulnerability to the heating may be species specific 44,45 and often depends on local adaptation 73 , the interpretation of unconsistent patterns could only be very oddly and speculative, until the relevancy of each single species is assessed. www.nature.com/scientificreports/ The present study uses the approach of relating different coralligenous communities to their thermal environment in order to contribute drawing the trajectories of change of the coralligenous community due to future warming scenarios. The innovative approach consists in providing some descriptors of heating events and the thresholds to which for this community would change. The same approach could be similarly applied to investigate the responses of any other subtidal community/species to climate change. The need of disentangling the effects due to the intensity, duration and variability of the heating events is stressed, as well as the importance of gaining spatially and temporally wide databases, both of the thermal environment and the biota. We hope the present contribution my also assist the implementing of persistent monitoring nets of subtidal habitats.   www.nature.com/scientificreports/ 166 days. The data obtained gave information on the variability in water stratification in the four study sites. Furthermore, 30-years satellite-derived SST data were obtained (AVHRR, https ://apdrc .soest .hawai i.edu) and used to detect the 2019 MHWs at each site and to explore the correlative pattern with subtidal water temperature.

Methods
Biota data collection. Samplings of the coralligenous community structure were done in May 2019 on vertical rocky walls using a non-destructive photographic method (underwater camera Lumix TZ30 with lighting achieved using two electronic strobes fitted with diffusers). At each site, two areas (about 100 m distant) were selected and at each depth (18 m, 23 m, 28 m, 33 m, 38 m) ten photographic samples of 0.2 m 2 of vertical surface were taken 74 . Organisms easily detectable from the photograph samples were identified to the lowest possible taxonomic level, while those not easily recognized were identified according to morphological groups 75 . In order to run reliable models, the cover of several unconspicuous taxa/groups was pooled in six categories (encrusting coralline algae, red algae, green algae, turf algae, sponges and bryozoans), based on the sample size and affinities ( Table 2). Red and green algae categories were mainly composed of Peyssonneliaceae and Udoteaceae, respectively, while turfs were algae smaller than 1 cm. The percent cover of each taxon/morphological group was estimated in each picture by ImageJ software (https ://image j.nih.gov/ij/). Richness of taxa/morphological groups and evenness index were calculated from the data acquired.
Temperature descriptors. A total of 26 temperature descriptors for all sitexdepth conditions were considered to estimating the duration of heating events, their intensity and temperature variability (Table 1). Overall, arbitrary thresholds were set choosing (1) temperature intensity for estimating S23-S27 (the number of days of daily temperature) and D23-D27 (the maximum heating duration) over 23-27 °C, respectively, (2) 4 °C as shift intensity for NHS and (3) 2 days of time for shifts in LTH. However, to estimate differences among sitexdepth conditions in terms of temperature anomalies, the daily temperature data were also explored to identify the number of days of maximum duration of a heating event (D90 and D95), the number of days of high temperature (S90 and S95) and the number of heating events (F90 and F95) respect to the 90th and 95th percentile, calculated on the temperature data collected during 2019. The temperature variability was measured by other several descriptors: NHS, the number of times that heating shifts occurred in two days in a row were larger than 4 °C; LTH, the largest heating event (ΔT, in °C) occurred in two days in a row; NDD, the number of daily shifts larger than the daily mean shift; LDH the largest daily heating (in °C, Table 1).

Data analysis.
Data exploration was carried out following 76 . Outliers were inspected with Cleveland dotplots and normality with histograms and Q-Q plots. Collinearity between continuous explanatory variables was inspected with pair-plots, and variance inflation factors (VIFs) were calculated. Several significant correlations were found among covariates: (1)  Since data exploration indicated non-linear relationships between each response variable and the explanatory variables, generalized additive models (GAMs) were run to correlate separately the different response variables with the explanatory variables, allowing detection of the eventual effects of Site and Depth on the continuous variables. GAMs are nonparametric extensions of linear regression models that allow the evaluation of highly non-linear relationships between explanatory and response variables thanks to the use of smooth functions 77 . To avoid over-dispersion, a negative binomial distribution was applied to all the response variables, with the exception of S (Poisson distribution) and evenness (gaussian distribution). The choice of the best fitting explanatory variables used in the final model was undertaken using AIC (Akaike information criterion), following a forward selection approach 76 . Models validation were run calculating and plotting the Pearson residuals against (i) the fitted values, (ii) each explanatory variable in the model, (iii) each explanatory variable not in the model 76  www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.