Sensitivity of a cold‐water coral reef to interannual variability in regional oceanography

We assessed the effects of regional oceanographic shifts on the macrofaunal biodiversity and biogeography of cold‐water coral reefs (CWCRs). CWCRs are often hotspots of biodiversity and ecosystem services and are in the frontline of exposure to multiple human pressures and climate change. Almost nothing is known about how large‐scale atmospheric variability affects the structure of CWCRs’ communities over ecological timescales, and this hinders their efficient conservation. This knowledge gap is especially evident for species‐rich macrofauna, a key component for ecosystem functioning.


| INTRODUC TI ON
Large-scale fluctuations in the atmosphere can affect the water mass characteristics and shape hydrography, biodiversity and food webs that are critical components of the marine environmental status (Mittermeier et al., 2011;Myrberg et al., 2019). In the North Atlantic, the dominant mode of atmospheric variability is the North Atlantic Oscillation (NAO), which can be described through the North Atlantic Oscillation Index (NAOI, hereafter) using an empirical orthogonal function of sea-level pressure anomalies over the Atlantic sector (Hurrell et al., 2003).
The NAO is a major driving force for basin-scale interannual variability in circulation and water mass characteristics. The NAO influences the strength and extent of the subpolar gyre (Lohmann et al., 2009) and upper ocean temperature/salinity values in the subpolar North Atlantic (Hátún et al., 2005;Johnson et al., 2013). A negative NAOI is associated with westward subpolar gyre retraction enabling the intrusion of relatively warm and salty subtropical waters into the eastern North Atlantic (Hátún et al., 2005;Johnson et al., 2013) and subsequently onto the western European shelf (Johnson et al., 2020 and references therein). The NAO can mediate the timing and biomass of primary production (Basterretxea et al., 2018), shifts in zooplankton biogeography and biodiversity (Beaugrand et al., 2002) and cause regime shifts in food webs (Molinero et al., 2008). However, studies on effects of the NAO on benthic marine biodiversity hotspots remain scarce.
Cold-water coral reefs (CWCRs hereafter) are structurally complex ecosystems characterized as hotspots of biodiversity, biomass, elemental cycling and socioeconomic values (Armstrong et al., 2019;De Clippele et al., 2021;Henry & Roberts, 2017). CWCRs' spatial distribution is mainly controlled by the availability of hard substrates, seawater temperature, chemistry and hydrodynamics which facilitate the supply of food (Roberts et al., 2006;Hebbeln et al., 2019). As a response to their environmental requirements, CWCRs are often found in the (upper) bathyal zone and continental shelves.
While CWCRs located at relatively shallow depths on the continental shelf might act as refugia for the effects of ocean acidification (Turley et al., 2007), they are still threatened by warming, physical damage and microplastics (Bindoff et al., 2019;La Beur et al., 2019).
As human activities expand, adaptive spatial management measures are required. This, in turn, needs a better understanding of variability in biodiversity and ecosystem functioning over space/time as well as the development of integrated monitoring programmes (Cormier et al., 2019).
Macrofauna are a key constituent in ecosystem functioning globally: they are extremely species-rich (Grassle & Maciolek, 1992) and shape benthic elemental cycling (Janas et al., 2019). Macrofaunal biodiversity and functional traits are usually correlated with a system's capacity to buffer environmental stressors and offer a useful proxy of marine ecosystem health (Gunderson, 2000;HELCOM, 2018;Hooper et al., 2005). Studies have examined the effects of size and type of biogenic habitat (Bourque & Demopoulos, 2018;Jensen & Frederiksen, 1992;O'Hara et al., 2008) as well as bathymetric and hydrographic gradients in shaping macrofaunal biodiversity (Henry et al., 2013;. Notably, these studies focussed on a single time point, and work examining macrofaunal communities in CWCRs over time is lacking. This is a major knowledge gap that hinders our understanding about the response of macrofauna to temporal environmental changes. Addressing this knowledge gap is crucial as studies over geological timescales have found compelling links between the decline of CWCRs' biodiversity and climatic/ oceanic perturbations (Douarin et al., 2014).
Our work aims to advance our understanding on the effects of spatial gradients in seafloor characteristics and current speed as well as the effects of interannual variability in regional oceanography on CWCRs' macrofauna. We focus on the Mingulay Reef Complex, a hotspot of marine biodiversity on the western European shelf.

Mingulay Reef Complex is included in the East Mingulay Marine
Protected Area (Figure 1). Increasing our understanding on the functioning of CWCRs enables efficient implementation of marine policies and the protection of these extremely fragile ecosystems for the future.

| Study area
The Mingulay Reef Complex (MRC, hereafter) is one of the beststudied CWCRs in the world and is located on the northeast Atlantic continental shelf, in western Scotland (e.g. Henry et al., 2013;Roberts et al., 2005). The MRC contains reef mounds that are 13-60 m wide, 16-108 m long and between 2 to 34 m tall, which are formed by the cosmopolitan scleractinian coral Lophelia pertusa (De Clippele et al., 2017). X-ray computed tomography of seabed cores revealed cycles of different ecosystem phases (reef initiation, framework expansion, framework collapse, coral rubble) over the mid-Holocene to late Holocene (Douarin et al., 2014 (Johnson et al., 2013;Porter et al., 2018). Over multi-annual time scales, the bottom temperatures and salinities at MRC show a statistically significant inverse relationship with the NAOI . Specifically, in years with weak atmospheric forcing and a negative NAOI, warmer and more saline bottom conditions are found at MRC, while cooler and fresher bottom values are seen during a positive NAOI . In contrast, the fresher Scottish Coastal Current has a minor influence at the MRC, particularly at the seabed (Inall et al., 2009). Considering data resolution and previous findings at MRC (NAOI-Hurrell et al., 2003;temperature, salinity-Findlay et al., 2013;bathymetry, topography, current speed, macrohabitat-Roberts et al., 2009;Henry et al., 2013), we classified a priori the bathymetry, topography and current speed as "spatial variables" (i.e. mainly changing across space), and the bottom potential temperature (hereafter mentioned as "bottom temperature"), bottom salinity and the NAOI as "temporal variables" (i.e. variables mainly changing across time).

| Benthos
Benthic samples (total n = 79) were collected using a modified Van Veen Grab (sampling 0.1 m 2 per deployment) in 2003,2005,2009,2010,2011 (Figure 1; Table S1). Samples were sieved at 1 mm, macrobenthos collected and stored in 4% seawater formalin and  Table S1 (a) (b) transferred to 70% industrial methylated spirit. Specimens were identified to the lowest possible taxonomic level using best available taxonomic keys for North Atlantic marine invertebrates and guidance from expert taxonomists.
Accounting for the existence of both colonial and solitary taxa, we have recorded species' presence/absence (Henry et al., 2013).
We acknowledge that the use of presence-absence (and not abundance/biomass) data may have some effect on the outcomes of our analyses but these effects are not expected to be major (Buchner et al., 2019). Furthermore, we feel that presence-absence data was the optimal approach as for many species (e.g. encrusting bryozoans) measurements on abundance and/or biomass would not be accurate.
Following previous biodiversity studies at MRC (Henry et al., 2013), we have also excluded sponges due to a lack of taxonomic resolution across this group. Presence/absence data were assembled in a species assembly matrix using PRIMER (Clarke & Gorley, 2015).

| Spatial variables: bathymetry, topography and current speed
Topographic terrain variables were calculated from bathymetric data (depth). Specifically, the following metrics were derived from the bathymetry data: slope, ruggedness (equivalent to rugosity), northness and eastness (via aspect), standardized broad-and fine-scale bathymetric position index (inner radius of 1 and outer radius of 5 and 3 cells, respectively). We used these bathymetry-related variables to have a good coverage of the seabed terrain parameters that change across space and to quantify their role in macrofaunal biodiversity patterns. Previous studies at the MRC (Henry et al., 2010(Henry et al., , 2013 and elsewhere (e.g. Savini et al., 2014;Sundahl et al., 2020)  Macrohabitat was also treated as a spatial variable. Based on previous studies at MRC (Henry et al., 2010(Henry et al., , 2013, macrohabitat was classified in six types, representing six levels of structural complexity. These types were as follows: "muddy sand" (i.e. the least structurally complex macrohabitat), "rubble," "rock," "live coral," "dead framework" and "mixed live and dead framework" (i.e. the most structurally complex habitat).

| Data analysis
To avoid statistical bias and Type I error inflation, all explanatory variables were checked for correlation (R package "corrplot" version 0.84; Wei & Simko, 2017) and variables with significant correlation scores (Pearson's r > .7 or <−.7) were discarded ( Figure S1). The species data were Hellinger-transformed, giving lower weights to rare species and maintaining linear correlations between species and gradients of environmental variables (Legendre & Gallagher, 2001).
The variance of macrofaunal communities explained by spatial and temporal environmental variables was examined using redundancy analysis (RDA). This multivariate statistical method combines regression and principal component analysis and can model response variables as a function of multiple explanatory variables (Borcard et al., 2018). The RDA model was built with the species assembly matrix as the response variable and all the non-correlated environmental variables as the explanatory variables. A forward stepwise selection was applied to retain the statistically significant explanatory variables (R package "vegan" version 2.5-6; Oksanen et al., 2019). In order to quantify the variance of macrofaunal communities explained by (a) spatial, (b) temporal and (c) interactions between spatial and temporal variables, a series of RDA and partial RDA (pRDA) models were constructed (Borcard et al., 1992).
The significance of fitted environmental vectors was determined by calculating the goodness of fit defined as the squared correlation coefficient r 2 using permutation tests with the command envfit() of package "vegan" (Oksanen et al., 2019). Species under the influence of each specific environmental variable were identified through visual inspection of the total RDA and pRDAs biplots, with species selection based on RDA scores (R package "vegan" and "RcolorBrewer" version 1.1-2; Neuwirth, 2014). All RDA and pRDA analyses were conducted in R. R code is provided in the Supporting Information.  (Burk, 2018). When the Kruskal-Wallis test was used, the pairwise comparisons were carried out through the Dunn test (Dinno, 2017). Accounting for the multiple comparisons, the p values were adjusted using the Bonferroni correction (Armstrong, 2014).

| Drivers of macrobenthic community structure
The forward stepwise selection of environmental variables retained four spatial and four temporal variables that had a statistically significant contribution (p < .05) in explaining macrofauna variance. In total, 26.7% of the variance was explained: 13.7% by the temporal variables, 7.2% by the spatial variables and 5.7% by the interactions between spatial and temporal variables. The spatial and temporal variables with a statistically significant contribution in explaining macrofaunal variance can be seen in Table 1.
Species richness (SR) did not show statistically significant changes over the interannual changes in bottom temperature (Figure 2; Table S4).

| Macrofaunal groups under the influence of environmental gradients across space
Species linked to higher fine-scale bathymetric position index values were mainly sessile suspension feeders, while those associated with structurally complex macrohabitats included a rich mixture of functional traits with predators (e.g. Porania (Porania) pulvillus, Eunice norvegica), mobile deposit feeders/grazers (e.g. Galathea strigosa) and sessile suspension-and filter-feeding species (e.g. the bivalve Hiatella arctica) (Figure 3a; see also Figure S2a and Table S5). Species found in highly rugged areas were mainly erect suspension feeders with larger more robust morphologies (e.g. the anthozoan Parazoanthus anguicomus), while about one-third of species found in less rugged areas were encrusting bryozoans (e.g. Schizomavella (Schizomavella) linearis). Species linked to high average current speed areas were mainly suspension-and filter-feeding species (the erect bryozoan Buskea dichotoma) and omnivores, while those in low average current speed areas were mainly deposit feeders (e.g. the isopod Uromunna petiti). In areas with higher structural complexity, there were more predators than suspension-and deposit-feeding species; on the contrary, in areas with lower structural complexity there were more suspension-feeding species than predators and deposit-feeding species (Figure 3a; see also Figure S2a and Table S5).

| Macrofaunal groups under the influence of environmental gradients across time
A single assemblage seemed to be associated with higher bottom temperatures (13 species) and another one associated with lower bottom temperatures (seven species) (Figure 3b; see also Figure S2b and  (Figure 3b; See also Figure S2b and Table S6).
Three species were associated with higher bottom salinity values (e.g. the isopod Janira maculosa) while 29 species (molluscs, polychaetes, bryozoans and echinoderms) (Figure 3b; see also Figure S2b and Table S6) were associated with slightly fresher bottom waters.
Seven species were associated with the higher NAOI (same year) (e.g. the suspension-feeding polychaete Serpula vermicularis) and the higher NAOI (previous year) values (e.g. the erect anthozoan Parazoanthus anguicomus). In contrast, a much larger number of species (n = 22) were associated with the lower values of the NAOI (same year). These were mainly encrusting suspension feeders (e.g. the bryozoan Amphiblestrum flemingii), small suspension feeders (e.g. the bivalve Modiolula phaseolina) and deposit feeders (e.g. the ophiuroid Amphipholis squamata) (Figure 3b; see also Figure S2b and Table S6).

| Macrofaunal functional traits
Almost all the groups of functional traits showed statistically significant changes over the interannual changes in bottom temperature, bottom salinity and NAOI (Figure 4; Table S7). The frequency (%) of the suspension/filter feeders and the sessile erect organisms had their lowest values during highest bottom temperatures (Figure 4).
This lowest frequency of suspension/filter feeders and sessile erect organisms was statistically different (p < .05) from the frequencies F I G U R E 4 Changes in (a) feeding types, (b) mobility/orientation types and (c) body size across years at the Mingulay Reef Complex. The outcomes of the pairwise comparisons are shown in Table S7 recorded under lower temperatures (Figure 4; Table S7). On the contrary, the frequency of predators had its maximum values under the highest bottom temperatures ( Figure 4); this peak in the distribution of predators was statistically different from the frequency of predators recorded under lower bottom temperatures (Figure 4; Table S7).

| D ISCUSS I ON
This study has unravelled the effects of both spatial and temporal environmental gradients on macrofaunal communities in a coldwater coral reef. Interannual changes in regional oceanography explained almost double the variance explained by spatial gradients, the latter which were solely considered in previous studies.
( Bourque & Demopoulos, 2018;Buhl-Mortensen et al., 1995;Cordes et al., 2008;Henry et al., 2013). Interestingly, there were no significant changes in SR over bottom temperature changes but there were changes in community composition, Δ + and functional traits and marginal changes for d. This highlights the necessity for thorough examinations of macrofaunal communities to identify changes in ecosystem structure and functioning.
Due to the relatively short time series examined here, it is challenging to draw a definite conclusion on the effects of temperature, salinity, and NAOI interannual variability on CWCR macrofauna. Models predict an increase of up to 2.0°C for bottom temperatures at Mingulay by 2100 (Findlay et al., 2013) while temperatures between 200-1,000 m in the North Atlantic are predicted to rise, on average, by 1.8°C over the same period (Puerta et al., 2020;Sweetman et al., 2017). Considering that changes in biodiversity and functional traits were observed over a range of 0.32°C, we highlight the need for the establishment of monitoring programmes at MRC and at CWCRs in general.

| On the role of spatial gradients
Many species associated with topographic highs or high average current speed were sessile suspension feeders suggesting that these conditions are highly favourable, possibly due to enhanced food supply .
Macrohabitat was also an important driver of species distribution. Areas dominated by live and dead coral framework hosted a richer mixture of functional traits than areas with muddy sand. This is likely related to the fact that higher structural complexity increases the number of ecological niches (Tews et al., 2004). The presence of predators (e.g. Eunice polychaetes) is also possibly related to the fact that they feed on anthozoans and tunicates which are abundant in coral framework . Regarding ruggedness, larger sized and erect suspension feeders were found in highly rugose areas while encrusting delicate forms proliferated in smoother areas. Large/erect/robust body forms likely proliferate in rugose areas under turbid waters due to enhanced food capture rates. In contrast, low-lying and delicate forms are preferably found in less energetic environments (Henry et al., 2013).

| On the role of temporal gradients in water mass characteristics and regional oceanography
A recent study showed a statistically significant inverse relationship between bottom conditions in the eastern North Atlantic and the strength of the NAO . The westward retraction of the subpolar gyre during a weak NAO (Lozier & Stewart, 2008), enables the northward extension of warm and salty subtropical waters (Hátún et al., 2005(Hátún et al., , 2017Johnson et al., 2013). A similar signal was observed at the MRC, with warmer and more saline bottom conditions again observed during a weak NAO . conditions enhance species' reproductive success (Byrne, 2011;Holopainen et al., 2016;Lewis et al., 2002;Rupp, 1973). Second, during a weak NAO/westward retraction of the subpolar gyre, the northward expansion of the warm and salty North Atlantic Water facilitates the extension of species from lower to higher latitudes.
Most species associated with warmer waters at MRC are mobile polychaetes and arthropods. Their higher mobility and pelagic larvae (Iannotta et al., 2007) could support their dispersal.
NAOI (same year) had a stronger contribution than the NAOI (previous year) in explaining biological variability, suggesting a relatively rapid response of macrofauna to NAOI. This mirrors findings from the intertidal/upper subtidal zone in the North Sea (Birchenough et al., 2018 and references therein), the Mediterranean bathyal (Cartes et al., 2009) and the abyssal northeast Atlantic (Smith et al., 2009) where the winter NAO caused shifts in the structure and functioning of benthic communities in subsequent seasons and years.
A high-resolution ocean model shows that the NAO has a positive correlation with bottom kinetic energy at the MRC . A high NAOI/increased kinetic energy likely enhances food supply supporting the proliferation of some suspension/filter feeders such as P. patelliformis and S. vermicularis (Table S6). In contrast, deposit feeders and encrusting organisms (Table S6) were associated with a low NAOI/decreased kinetic energy which likely enhances the food capture by these low-lying species.
It is important to note that a large part of biological variability remains unaccounted for. About 80% of the polychaetes species are distributed over a relatively wide range of bottom temperatures and salinities (OBIS, 2019). Thus, environmental gradients experienced at the MRC are likely within their niche envelope and will not explain any variability. Unexplained variability could also be linked to limited dispersal. Indeed, 40% of the species not associated with environmental gradients were found at only one station. Limited dispersal obscures any relationship between biological and environmental sampling scales (Bullock et al., 2006;Selmoni et al., 2020).
Third, the temporal resolution of environmental sampling may be somewhat coarse and biased towards species with perennial life spans.

| Implications for biogenic ecosystem functioning
The changes seen at the MRC macrofauna over time, including common species (Henry et al., 2013;, indicate that the macrofauna may be susceptible to interannual changes in water mass characteristics. When the MRC was subjected to higher bottom temperatures, there was a significant decrease in suspension/filter feeders and a significant increase in predators (Figure 4; Table S7). The mechanisms driving these shifts are not known, and longer time series are needed to draw a firm conclusion. It is noteworthy, though, that ~70% of the suspension/filter-feeding species that were not recorded during higher temperatures have to date also not been recorded in (sub)tropical areas. This indicates their preference for higher latitudes and colder waters (Hayward & Ryland, 2017;Jensen & Frederiksen, 1992;Klitgaard, 1995). In addition, predators recorded during higher bottom temperatures at MRC have been found in (sub)tropical areas indicating their affiliation with a wide range of bottom temperatures (Froglia & Speranza, 1993;Simon et al., 2014).
Suspension and filter feeders are usually primary consumers which channel organic matter from the lower to the higher trophic levels. Studies in tropical reefs and CWCRs have shown their key role in food web functioning through the high assimilation of dissolved organic carbon and its conversion to biomass (de Goeij et al., 2017;Rix et al., 2016). Their role is also pronounced in areas where strong bottom currents prevail, for example CWCRs, as they trap food sources which otherwise would have been transported away (Gili & Coma, 1998;Roberts et al., 2009  and thus, the establishment of regular monitoring programmes is necessary for reef conservation. Our findings contribute also to the conservation of CWCRs on a regional scale. Specifically, biological and environmental (meta)data from MRC are archived in online repositories (e.g. PANGAEA) serving the upscaling of findings from the local to the regional level, which is necessary for the achievement of Good Environmental Status (Kazanidis et al., 2020). In addition, our findings contribute to the identification and management of climate change refugia, a mitigation strategy for the conservation of species susceptible to climate change (Morelli et al., 2016). The shoaling of the aragonite saturation horizon by 2100 (Orr et al., 2005) is expected to cause major reductions in habitat suitability for North Atlantic CWCRs as they will be bathed in more acidic waters . However, CWCRs on the continental shelf, such as the MRC, are expected to be less affected by the shoaling of the aragonite saturation horizon, and thus, they may serve as climate change refugia Turley et al., 2007). Data continuity, accessibility and affordability are important obstacles in the identification of fine-scale climate change refugia (Morelli et al., 2016).

| Implications for CWCR conservation
Advancing knowledge on the temporal dynamics of continental shelves' CWCRs and sharing of gathered (meta)data serve the establishment of monitoring programmes and CWCR conservation for the future.

ACK N OWLED G EM ENTS
We would like to thank the captains and crews of the expeditions involved in the collection of samples used in the present study. We would like also to thank the following taxonomists for providing valuable guidance in the identification of the speci- his comments regarding oceanographic settings in the Scottish continental shelf and finally the captains and crews of the vessels used for the collection of benthic samples.

CO N FLI C T O F I NTE R E S T
The author declare there is no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13363.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the data used in the present study have already been deposited