Decades of biomass loss in the shallow rocky subtidal vegetation of the south-eastern Bay of Biscay

This study seeks to assess changes over time in the structure of subtidal macroalgal assemblages across depth in the south-eastern Bay of Biscay. The results reveal a large-scale decline in total macroalgal biomass between 1982 and 2014. However, the temporal pattern of shift differs from one depth to another: total biomass decreased at depths from 3 to 10 m, but increased at depths of 2 and 11 m. The strong decrease in biomass detected in the 3–10-m depth range is a consequence of a sharp net decline in large macroalgae biomass which was not offset by increased biomass of small species, mainly corresponding to turf-forming algae. The dominant canopy-forming Gelidium corneum in 1982 had practically disappeared by the end of the study period and its biomass loss was far from being offset by the small increase detected in the fucoid Gongolaria baccata. By contrast, at depths of 2 and 11 m, the most notable result is an increase in large species, mainly Halopithys incurva and Codium decorticatum at 2 m and G. baccata at 11 m; however, at both depth levels, a new canopy was far from being developed. These findings evidence that biomass and habitat provision, two pivotal roles of canopy-forming species in ecosystem functioning, have been altered. Further research into potential changes in primary productivity and biodiversity linked to the shift detected in assemblage structure needs to be conducted in order to get information for conservation and management decisions associated with the loss of habitat-forming macroalgae.


Introduction
Canopy-forming macroalgae usually form extensive stands in rocky benthic subtidal communities in most temperate regions (Steneck et al. 2002;Smale et al. 2013;Strain et al. 2014). These large, often perennial macroalgae play a very important role in marine ecosystems as they create structurally complex assemblages analogous to forests on land (Ballesteros et al. 2009;Reed and Foster 2012;Gianni et al. 2013), i.e. they are foundation species. In this regard, these macroalgae increase three-dimensional complexity by providing biogenic habitats and protection for a great variety of marine organisms (Steneck et al. 2002;Wernberg et al. 2011;Smale et al. 2013). They also act as ecosystem engineers, since they substantially modify the local environment by changing light conditions (Wernberg et al. 2005), water flow (Rosman et al. 2007) and sedimentation rates (Eckman et al. 1989) in ways that favour the settlement of other organisms. Canopydominated assemblages constitute some of the most diverse and productive ecosystems anywhere in the world, contribute significantly to nearshore primary productivity and also enhance secondary productivity (Mann 2000;Steneck et al. 2002;Tait and Schiel 2011;Smale et al. 2013). In addition, these foundation species supply many valuable ecosystem services such as reducing coastal erosion, CO 2 sinking, nutrient cycling and water quality control (Airoldi and Beck 2007;Smale et al. 2013;Wernberg et al. 2016).
However, in the last 30 years, drastic declines in these foundation species have been observed as a consequence of multiple anthropogenic pressures comprising harvesting, pollution, sedimentation, invasive species, overgrazing due to a decline in grazer predators, fishing nets, recreation and ocean warming (Steneck et al. 2002;Serisawa et al. 2004;Connell Communicated by A. F. Bernardino et al. 2008;Perkol-Finkel and Airoldi 2010;Smale et al. 2013;Mineur et al. 2015;Krumhansl et al. 2016;Wernberg et al. 2016). Whatever is driving the decline, a shift to less structurally complex communities dominated by turf-forming, filamentous or ephemeral macroalgae has been widely documented Mangialajo et al. 2008;Perkol-Finkel and Airoldi 2010;Tait and Schiel 2011). To date, there has been little evidence of foundation species recovering from disturbances (Dayton et al. 1992, Martínez andCárdenas 2003), mainly because once less structurally complex communities are established they may inhibit recolonisation by canopy species, resulting in alternative stable states (Strain et al. 2014;Wernberg et al. 2016). Loss of large perennial macroalgae in favour of less structured communities may have consequences for the whole ecosystem, and in this connection, numerous studies have linked the retreat of these ecosystem engineers to a decrease in species richness and abundance (Graham 2004;Norderhaug et al. 2007;Schiel and Lilley 2007;Wikström and Kautsky 2007), homogenisation with neighbouring habitats (Mangialajo et al. 2008) and loss of productivity (Tait and Schiel 2011;Crowe et al. 2013).
In the particular case of the south-eastern Bay of Biscay, local retreats in subtidal perennial canopy-forming Gelidium corneum, Laminaria ochroleuca and Gongolaria baccata (formerly Cystoseira baccata) have been reported in the past few decades (Díez et al. 2012;Borja et al. 2013Borja et al. , 2018Muguerza et al. 2017Muguerza et al. , 2020. In addition to canopy decline, there has been an increase over the same period in richness and abundance of warm-affinity species (mainly ephemeral forms with simple morphology), coralline algae and crustose species have become abundant, and non-indigenous species have expanded (Díez et al. 2012;Muguerza et al. 2017). These studies point out that higher temperatures are probably the main driver of the changes observed. Indeed, a warming of 0.26 ± 0.03°C every 10 years was detected in the Bay of Biscay for 1982-2014 (Costoya et al. 2015). However, other local factors such as nutrient availability, solar radiation, sunlight hours and wave height have been suggested as potential co-acting factors of change in combination with warming (Díez et al. 2012;Borja et al. 2013Borja et al. , 2018Muguerza et al. 2017).
Most of the aforementioned research papers provide information on species abundance in terms of cover, but little information about loss of biomass due to canopy decline is available for the south-eastern Bay of Biscay. Only Borja et al. (2013) provide data on changes in G. corneum biomass for the eastern Basque coast, where a drastic reduction in its standing stock of about 7800 t across 30 km has been documented for 1993-2012. However, no data on the variability over time of the biomass of other species and areas are available. This paper seeks to assess changes over time (three sampling surveys: 1982, 2007 and 2014) in the structure of subtidal macroalgal assemblages in terms of composition and taxon biomass in the westernmost part of the Basque coast. These assemblages can be considered representative of the south-eastern Bay of Biscay since they are distributed along large stretches of coastline in this region (Gorostiaga et al. 1998;Díez et al. 2003).

Study area
The study area lies at the eastern end of the Cantabrian Sea (Northern Spain), on the south-eastern Bay of Biscay. It is over 192 km in length and is open to strong waves coming mostly from the NW, with an average height of 1.9 m (Díez et al. 2003;González et al. 2004;Galparsoro et al. 2010). The south-eastern Bay of Biscay is exposed to highly exposed to the prevailing NW swells, with high, mostly erosional, energy. It also features extensive vertical cliffs and abrasion platforms interspersed with sandy beaches. In shallow waters in the study area, the rocky bottom is almost continuous, but it becomes sandy as depth increases (Chust et al. 2011). The flora belongs to the warm temperate NE Atlantic Region according to the biogeographical scheme proposed by van den Hoek and Breeman (1990).

Field sampling and processing
The biomass of taxa composing in macroalgal assemblages was studied at six different depths (2, 3, 6, 9, 10 and 11 m) along seven transects in three sampling surveys (1982, 2007 and 2014). Each transect was set perpendicular to the coastline following a north-south orientation with a starting point at 2 m below extremely low tides (Fig. 1). Within each transect, a surface area of 2000 cm 2 was delimited at each depth using quadrats of 40 × 50 cm placed systematically at the midpoint of the cross section of the transect. This means that there were six quadrats per transect across each depth profile and sampling survey. In some transects, it was not possible to sample at certain depths due to a lack of appropriate substrate (continuous bedrock with slight to moderate slopes of <30°). Each surface was destructively sampled, with all macroalgal species within the quadrat being collected except the mostly calcareous crustose layer, which was not sampled. Once in the laboratory, samples were kept frozen in labelled plastic bags. For analysis, samples were thawed and the macroalgae were separated and identified. Algal taxonomy was updated following AlgaeBase (Guiry and Guiry 2021). The dry weight (DW) in grammes (100-110°C, 24 h) was then obtained for each species.

Statistical analysis
Multivariate analyses were performed in order to explore the spatio-temporal variability of the structure of assemblages in terms of composition and taxon biomass. Prior to analysis, the biomass values for each taxon were square root transformed to reduce the influence of the dominant ones. The similarity between pairs of samples was calculated using the Bray-Curtis index. The hypothesis that time and depth have no influence on the structure of macroalgal assemblages was tested by means of PERMANOVA (permutational multivariate analysis of variance, see Anderson et al. 2008) with an a priori chosen significance level of α = 0.05. The design of the experiment was as follows: time (Year; set with three levels: 1982, 2007 and 2014) and depth (Depth; random with six levels: 2, 3, 6, 9, 10 and 11). Post hoc pairwise comparisons were performed using Gosset's t-statistic to investigate the significant terms of the PERMANOVA (Anderson et al. 2008). In order to graphically visualise this spatio-temporal variation, a non-metric multidimensional scaling (nMDS) was conducted. Given that PERMANOVA tests the null hypothesis that centroids and/or dispersion of the groups defined by the factors of the experimental design is equivalent, we performed a permutational test for homogeneity of multivariate dispersions (PERMDISP) to check for differences in dispersion between the levels within the Year factor. By applying classification analysis (CLUSTER), samples were segregated into different groups according to their similarities. This analysis was followed by a similarity percentage (SIMPER) analysis to calculate the contribution of each taxon (%) to the dissimilarity between the clustering groups. All statistical analyses and the aforementioned routines were performed using the PERMANOVA+ for PRIMER6 software package (Clarke and Gorley 2006).

Results
The survey carried out in 1982 identified a total of 55 taxa, 94 species were documented in 2007 and 65 in 2014. Across all three surveys, the most widely represented phylum was Rhodophyta, with a total of 92 taxa, followed by Ochrophyta with 13 and Chlorophyta with 12 (Table 4).
PERMANOVA results (Table 1) show the partitioning of sample variation in the multivariate space based on the Bray-Curtis similarity in response to the Year and Depth factors. Differences between depths in the pattern of temporal change were detected, given that the Year x Depth interaction was significant (p = 0.0028). Pairwise comparisons reveal that at a depth of 2 m, significant changes in the structure of assemblages start to occur in 2007. At depths of 3 and 6 m, significant differences are detected from 1982 onwards, whilst at depths from 9 to 11 m, the structure of assemblages changed between 1982 and 2007 but has not changed significantly since then (Table 1). Non-metric multidimensional scaling (nMDS) (Fig. 2) shows the PERMANOVA results in graphic form. A gradual shift in the structure of assemblages is found across the sampling surveys (1982, 2007 and 2014 In the dendrogram resulting from classification analysis (CLUSTER) of samples, two main groups, A and B, are distinguished at a similarity level of 14% (Fig. 3). Group A (average similarity between samples 42.9%) consists mainly of samples from 1982 but also contains some samples from 2007 and one from 2014. The samples from 2007 are all from depths of 3 and 6 m except for two which come from 2 and 9 m (Fig. 3). Group B (average similarity between samples 24.8%) comprises the rest of the samples from 2007 and 2014 and two samples from 1982 taken at a depth of 2 m (Fig. 3). The similarity percentage (SIMPER) analysis points to Gelidium corneum as the main species responsible for the separation of these two groups, with a contribution of 26.9%. The mean biomass of this rhodophyte is higher in group A. Apart from G. corneum, other species which make noteworthy contributions (>2%) to group A are Pterosiphonia complanata, Plocamium cartilagineum, Asparagopsis armata and Dictyopteris polypodioides. Of these species, P. complanata, P. cartilagineum and D. polypodioides are more abundant in group A, whilst the biomass of A. armata is similar in both groups. In addition to the lower presence of G. corneum, group B also shows a mosaic distribution of many species, dominated by Gongolaria baccata, Codium decorticatum, Halopithys incurva and Corallina spp. It is also noteworthy that this group shows a greater abundance of morphologically simple forms such as Aphanocladia stichidiosa, Lychaete pellucida, Aglaothamniom pseudobyssoides or Microcladia glandulosa, among others (Table 3).
Group A is divided into two subgroups (A1 and A2) at a similarity level of 24% (Fig. 3). Subgroup A1 consists of   (Table 3). By contrast, subgroup A2 comprises samples from 1982, mainly from deeper waters. In this case, the abundance of G. corneum is notably lower and the species that shows the greatest abundance is P. complanata (Table 3). Other species that show considerable biomass values are A. armata, D. polypodioides, Calliblepharis ciliata and Heterosiphonia plumosa, all of which are more abundant in subgroup A2 than in subgroup A1. The latter subgroup (samples from depths of 3 and 6 m) is divided, in turn, into two groups (A1.1 and A1.2) at a similarity level of 59% (Fig. 3) (Table 3). Group B is divided into two subgroups (B1 and B2) at a similarity level of 17% (Fig. 3). Subgroup B1 is represented mainly by deeper samples from 2007 and 2014, with G. baccata as the dominant species (Table 3). Other species with noteworthy biomass values include Phyllophora crispa, and Corallina spp. Subgroup B2 mainly comprises shallower Table 3 Summary of the SIMPER test indicating the average biomass (Av.Bio; g DW ·2000 cm −2 ) of each taxon and its contribution (C (%)) to the differentiation of the subgroups identified in the classification analysis. The biomass values expressed in this table are untransformed.
The trend over time across depths of the mean total biomass and the biomass of those species that SIMPER indicates contribute most (>7%) to the groups detected in the CLUSTER which is shown in Fig. 4. The highest total biomass figures for 1982 are found at depths of 3 and 6 m with 194.56 and 154.47 g DW · 2000 cm −2 , respectively (Fig. 4). Intermediate values are found at depths of 9 and 10 m (77.36 and 71.51 g DW · 2000 cm −2 ) and the lowest values at depths of 2 and 11 m (46.18 and 30.18 g DW ·2000 cm −2 ). These biomass records in 1982 at depths of 3, 6, 9 and 10 m are consistent with the high abundance of G. corneum (189.36, 149.51, 63.74 and 43.67 g DW ·2000 cm −2 , respectively). However, at 2 m, H. incurva is the most abundant species (18.65 g DW ·2000 cm −2 ), whilst at 11 m, the biggest contributor to biomass is P. complanata (7.31 g DW · 2000 cm −2 ). In 2007 and 2014, the total biomass does not exceed 93 g DW · 2000 cm −2 at any depth and it is more homogeneously distributed than in 1982.

Discussion
This research reveals a large-scale decline in total macroalgal biomass between 1982 and 2014. However, the temporal pattern of the shift differs from one depth to another: total biomass decreased at depths from 3 to 10 m but increased at 2 and 11 m. The strong decrease in biomass detected in the 3-10-m depth range is a consequence of the sharp net decline in large macroalgal biomass, which was not offset by the increased biomass of small species, mainly turf-forming forms. The dominant large macroalga in 1982, the canopy-forming G. corneum, had practically disappeared by the end of the study, and its biomass loss was far from offset by the small increase detected in the fucoid Gongolaria baccata. By contrast, at depths of 2 and 11 m, the most notable result is the increase in large species, mainly Halopithys incurva and the annual chlorophyte C. decorticatum at 2 m, and G. baccata at 11 m; however, at both depth levels, a new canopy is far from being developed.
The decline of G. corneum detected in this study is in line with that recorded along the Basque coast in the last 30 years (Díez et al. 2012;Muguerza et al. 2017). Increased water temperature, irradiance and wave height are suggested as being among the main drivers underlying this regression (Díez Fig. 4 Mean biomass (g DW · 2000 cm −2 ) of the whole community (Total) and the five species contributing most (>7%) to differences detected between surveys at different depths. At depths of 3 and 6 m (*), the biomass ranges from 0 to 200 g dry weight · 2000 cm −2   At the end of the study period, scattered individuals and small patches of the large macroalga G. baccata were detected. However, the ability of this fucoid to colonise shallow rocky reefs in coastal stretches exposed to strong waves is rather limited (Díez et al. 2003). Accordingly, the increase in wave energy detected since the early 1990s on the Basque coast (Borja et al. 2013) is suggested as the main factor of change explaining the decline of G. baccata in some pristine locations along this coast (Muguerza et al. 2017). This macroalga is a warm-temperate species whose distribution is expected to expand northwards in the context of ongoing climate change (Hiscock et al. 2004). It has a high capacity for acclimation to increased temperature and irradiance levels (Miguel-Vijandi et al. 2010), which may be the reason for the lengthening of its growth period and the increase in biomass detected in semi-exposed coastal stretches of the southern Bay of Biscay since 2007 (Méndez-Sandín and Fernández 2016). Likewise, the physiological traits of G. baccata may explain the increase detected in the present study at a depth of 11 m, where the bottom friction exerted by waves is low enough for G. baccata to cope with it, whilst the increases in water temperatures and irradiance registered in the study area (Quintano et al. 2019) may have favoured its development.
Concurrently with the strong net decrease in canopyforming macroalgae at depths of 3-10 m, turf-forming species comprising articulated coralline macroalgae (Corallina spp., J. rubens) and morphologically simple and ephemeral algal forms (Aphanocladia stichidiosa, Lychaete pellucida, Aglaothamniom pseudobyssoides, M. glandulosa) have expanded. By contrast, at a depth of 2 m, where the warmtemperate affinity species H. incurva and C. decorticatum increased, and at the depth of 11 m where G. baccata increased, no expansion of turf-forming species was detected. Large macroalgae competitively exclude some species by monopolising resources, particularly space (Maggi et al. 2009), and numerous studies indicate that opportunistic species readily colonise space made available by canopy loss (Benedetti-Cecchi et al. 2001;Bulleri et al. 2002;Airoldi et al. 2008). The net loss of biomass linked to the shift detected in assemblage structure may be long-lasting (Tait and Schiel 2011;Crowe et al. 2013), since once turfing vegetation is well established, it may inhibit canopy-forming macroalgae from recruiting . In this regard, a shift to stable turfing assemblages in the vicinity of the study area has already been documented (Díez et al. 2014).
Canopy-forming macroalgae act as ecosystem engineer species and directly or indirectly modulate the availability of resources (Steneck et al. 2002). One consequence of canopy loss is a decline in the richness and abundance of the associated flora and fauna (Steneck et al. 2002;Graham 2004;Norderhaug et al. 2007;Schiel and Lilley 2007;Wikström and Kautsky 2007). In the study reported here, this functional role provided by G. corneum in 1982 was not replaced in 2014 by the expansion of H. incurva and C. decorticatum at a depth of 2 m, or by G. baccata at 11 m, since at both these levels, a new canopy was far from being developed. The complex habitat formed by G. corneum is essential for the functioning of the ecosystem since it preserves understory and epiphytic assemblages of smaller macroalgae, as well as sessile and vagile invertebrates (Borja et al. 2004;Bustamante et al. 2017). This habitat has become extinct in most of the study area, which may have consequences for the whole benthic ecosystem .
Branched canopies of species of the genus Cystoseira (G. baccata was formerly Cystoseira baccata) increase coastal primary production (Ballesteros et al. 2009), preserve biodiversity (Bianchelli et al. 2016), offer nursery areas for juvenile fish (Cheminée et al. 2013) and provide a home for outstanding species richness and density in coastal fish assemblage (Orlando-Bonaca and Lipej 2005). Therefore, the potential development of deepwater forests of G. baccata may take over the ecosystem functional role previously played by G. corneum. However, at present, this functional replacement has not yet occurred, since G. baccata stands are poorly developed. In addition, the newly established turf-forming macroalgae typically consist of species with less ecological and functional value than those replaced (Crowe et al. 2013). A decrease in the richness and abundance of associated organisms may therefore be expected. Furthermore, recent research in the study area has shown a sharp decrease in invertebrate taxonomic and functional density and diversity after canopy loss (Bustamante et al. 2014. Although species richness increased only slightly by the end of the study period for this research, significant increases (mainly in ephemeral forms) have been detected along the Basque coast (Muguerza et al. 2017(Muguerza et al. , 2020. The latter finding has been related to the intermediate disturbance hypothesis, which predicts maximal diversity at intermediate levels of disturbance (Connell 1978).
Foundation species are very important for marine food webs as they facilitate the capture and export of carbon (Dayton 1985;Krumhansl and Scheibling 2012;Smale et al. 2013). Moreover, much of their biomass is not consumed directly by herbivores, so canopy species are a major source of nutrition for other nearshore ecosystems (Duggins and Eckman 1997, Mann 2000, Steneck et al. 2002. Although no specific productivity measurements are made in this study, there may be a high risk of future primary productivity being significantly lower if the biogenic habitat previously provided by G. corneum is not recovered. Thus, previous experimental research on the impact of canopy loss on ecosystem functioning (Tait and Schiel 2011) has reported long-term reductions in primary productivity of macroalgal assemblages following canopy removal.
In addition to changes in ecosystem functions, canopy loss and subsequent decrease in macroalgal biomass in the study area could also impair the ecosystem services that canopy species provide. As mentioned above, these ecosystem engineers provide refuge for numerous other species, including many that are economically important for humans (Graham 2004;Smale et al. 2013). In this regard, the fact that the species that have replaced G. corneum cannot perform the same functional role means that the exploitation of some species in the study area might decline. G. corneum used to be the main raw material for agar extraction along the Atlantic shores of Spain, Portugal and Morocco (McHugh 1991). In the particular case of the study area, this resource has not been exploited since 1999 due to its decline (Borja et al. 2013).
In summary, the findings reported here evidence that two pivotal roles of canopy-forming species in coastal ecosystem functioning (biomass and habitat provision) have been altered in shallow rocky bottoms in the south-eastern Bay of Biscay. It seems that at present there are no native canopy-forming species that can adapt to the environmental Table 5 List of all the taxa recorded, showing their mean biomass (g DW ·2000 cm −2 ) for the six depths and the three sampling surveys.  conditions in the 3-10-m depth range where G. corneum previously thrived. Further research on potential changes in primary productivity and biodiversity linked to the shift detected in assemblage structure needs to be conducted to obtain information of elementary importance for conservation and management decisions associated with the loss of habitat-forming macroalgae.