Effects of Natural and Anthropogenic Stressors on Fucalean Brown Seaweeds Across Different Spatial Scales in the Mediterranean Sea

Algal habitat-forming forests composed of fucalean brown seaweeds (Cystoseira, Ericaria, and Gongolaria) have severely declined along the Mediterranean coasts, endangering the maintenance of essential ecosystem services. Numerous factors determine the loss of these assemblages and operate at different spatial scales, which must be identified to plan conservation and restoration actions. To explore the critical stressors (natural and anthropogenic) that may cause habitat degradation, we investigated (a) the patterns of variability of fucalean forests in percentage cover (abundance) at three spatial scales (location, forest, transect) by visual estimates and or photographic sampling to identify relevant spatial scales of variation, (b) the correlation between semi-quantitative anthropogenic stressors, individually or cumulatively (MA-LUSI index), including natural stressors (confinement, sea urchin grazing), and percentage cover of functional groups (perennial, semi-perennial) at forest spatial scale. The results showed that impacts from mariculture and urbanization seem to be the main stressors affecting habitat-forming species. In particular, while mariculture, urbanization, and cumulative anthropogenic stress negatively correlated with the percentage cover of perennial fucalean species, the same stressors were positively correlated with the percentage cover of the semi-perennial Cystoseira compressa and C. compressa subsp. pustulata. Our results indicate that human impacts can determine spatial patterns in these fragmented and heterogeneous marine habitats, thus stressing the need of carefully considering scale-dependent ecological processes to support conservation and restoration.


INTRODUCTION
Predicting the effects of abiotic and biotic stressors on marine vegetation changes has been a central concern of marine science in the last decades world-wide (Lotze et al., 2006;Airoldi and Beck, 2007;Coleman et al., 2008), and in the Mediterranean sea (Viaroli et al., 2008;Macreadie et al., 2014;Papathanasiou et al., 2015;Tsioli et al., 2019). Among the different algal groups, special attention has been devoted to canopy-forming fucalean brown algae (Strain et al., 2014;Coleman and Wernberg, 2017), characterized by slow growth and limited propagule dispersal and, therefore, unable to respond rapidly to anthropogenic and climatic changes (Buonomo et al., 2017;Bermejo et al., 2018). In turn, since the fucalean brown algae form some of the most productive, diverse, and valuable marine habitats (Benedetti-Cecchi et al., 2001), providing essential ecosystem services (Cheminée et al., 2013;Mineur et al., 2015;Bianchelli and Danovaro, 2020), they are recognized as a subtype of a natural habitat type (Reefs, code 1170) in need of conservation in Europe (Directive 92/43/EEC). Data from laboratory experiments with artificial seagrass showed that leaf area index (LAI), which combines the effect of leaf length and shoot density and therefore can be used as a metric of meadow abundance, was positively related to the delivery of ecosystem services such as wave attenuation (Paul et al., 2012), highlighting the need for a better understanding of the scale based abundance patterns (Townsend et al., 2018) under marine pristine and degraded ecosystems (Benedetti-Cecchi et al., 2001).
Algal forests composed by the genera Cystoseira, Ericaria, and Gongolaria (Fucales), with the exception of C. compressa, likely represent the most endangered habitat in the Mediterranean Sea (Barcelona Convention-Annex II; United Nations Environment Programme/Mediterranean Action Plan-UNEP/MAP; Verlaque et al., 2019) and have undergone a major decline in the last decades (Thibaut et al., 2005;Blanfuné et al., 2016;Rindi et al., 2020). Although the processes driving fucalean species diversity, abundance, and coexistence along environmental gradients are still unclear, the cumulative impacts of local anthropogenic stressors such as coastal development, habitat destruction, pollution, and fisheries along with climate change are considered the main causes Gianni et al., 2013;Strain et al., 2014;Mineur et al., 2015;Bianchi et al., 2018;Fabbrizzi et al., 2020). In addition, human overexploitation of top predators of sea urchins (e.g., the sparid fishes Diplodus sargus and D. vulgaris (Guidetti, 2004), has reduced the control over this component, thus leading to sea urchin overgrazing on algal forests, creating barrens (Ling et al., 2015). In addition, disruption of climatic patterns with more frequent and stronger storm events, as well as heat waves, which are becoming increasingly frequent during summer in the whole Mediterranean Sea (Darmaraki et al., 2019), are expected to have severe effects at local scales on fragmented populations of fucalean species (Verdura et al., 2021).
The knowledge required to disentangle the key drivers and abiotic variables that influence the distribution and abundance of fucalean seaweed species is still incomplete, mainly when the link between species decline and environmental conditions is considered (de Caralt et al., 2020). Therefore, understanding the reasons of the persistence and loss of benthic forests is crucial to evaluate the role of ecological and evolutionary processes and to manage direct and indirect anthropogenic stressors.
The determination of the most appropriate spatial scales for investigating species/environment relationships, directly related to the scale of ecological processes, is another crucial issue (O'Neill et al., 1989). For example, to explain differences in fucalean seaweed abundance and distribution, an understanding of the underlying ecological processes is necessary (Rindi and Guiry, 2004;Orfanidis et al., 2008), and therefore, quantification of their spatial patterns is needed (Wiens et al., 1993;Underwood, 1996;Benedetti-Cecchi et al., 2001;Fraschetti et al., 2005). Despite this interest in scale-specific patterns, effective implementation of multiscale approaches in theoretical and empirical research is still limited (Benedetti-Cecchi et al., 2001;Fraschetti et al., 2005;Mancuso et al., 2018).
To explore the stressors (natural and anthropogenic) estimated semi-quantitatively that may cause fucalean habitat degradation, we studied the patterns of variability of 11 species and subspecific taxa in abundance (expressed as percentage cover) to identify relevant scale of spatial variation. Indeed, five species or subspecies of the genus Cystoseira [C. compressa (Esper) Gerloff and Nizamuddin, C. compressa subsp. pustulata (Ercegovic) Verlaque, C. corniculata (Turner) Zanardini, C. foeniculacea (Linnaeus) Greville, C. foeniculacea f. tenuiramosa (Ercegovic) A. Gómez Garreta, M.C. Barceló, M.A. Ribera, and J. Rull Lluch], four of the genus Ericaria [E. amentacea (C. Agardh) Molinari and Guiry, E. barbatula (Kützing) Molinari, and Guiry, E. crinita (Duby) Molinari, and Guiry, and E. mediterranea (Sauvageau) Molinari and Guiry], and two of the genus Gongolaria [G. barbata (Stackhouse) Kuntze, G. elegans (Sauvageau) Molinari, and Guiry] were studied. The study was performed at three spatial scales (ranging from a few meters to 10 s of kilometers) in five countries and eight different locations across the Mediterranean Sea (2 in Greece; 1 in Albania; 2 in Italy; 2 in Spain; 1 in Turkey). In terms of functional morphology, C. compressa and C. compressa subsp. pustulata are semi-perennial or pseudo-perennial (i.e., algae in which most of the thallus is lost every year and the species persists in unfavorable seasons in the form of a small holdfast), whereas the others are perennial (i.e., algae in which a more or less large part of the thallus persists continuously for many years). We also explored the potential correlation between stressors individually or cumulatively (MA-LUSI index), including natural stressors (confinement, sea urchin grazing), and the functional groups (perennial, semi-perennial) percentage cover at forest spatial scale.

Sampling Locations
Sampling was conducted selecting randomly from a number of locations of the Mediterranean Sea where the presence and longterm persistence of the fucalean genera Cystoseira, Ericaria, and Gongolaria was documented by previous data or observations (Figure 1). A description of the locations is provided in Supplementary Material.

Sampling Design
Cystoseira, Ericaria, and Gongolaria forests were sampled following a random nested design with a hierarchy of three spatial scales: transect, forest, location (Figure 2). An exception was the locations I_GNI and T_DIDIM, where the sampling was implemented at two spatial scales, i.e. transect and forest, due to the limited extent of the only forest present. Sampling was undertaken during the warm time of the year (corresponding in the Mediterranean to late spring and summer), from May to October 2019 (and in July 2020 for T_DIDIM). In this period, sea surface temperature ranges between 20 and 25 • C in most of the Mediterranean. This is the most suitable time of the year to survey fucalean algae and evaluate their percentage cover, because both perennial and semi-perennial species occur in the field in their fully developed habit (consisting of a crust-or disk-like holdfast, a stipe and many branched fronds). During this period, in each location, one or more forests were selected, representing a spatial scale that equals the site's spatial scale. Several transects parallel to the shore were selected randomly within each forest, depending on the forest's extent (Figure 2). A georeferenced 10m graduated transect line was placed in a randomly selected spot within a forest down to 1-1.5 m depth. Three to twelve (n = 3-12) metallic or PVC quadrats (20 × 20, 25 × 25, 30 × 30, or 50 × 50 cm) were placed randomly in parallel to each transect at 0.5 m distance from the line. The fucalean species within the quadrats were recorded, and their percentage cover was estimated by visual census or photographs that were subsequently analyzed by means of imaging software (PhotoQuad; Trygonis and Sini, 2012). While the two forests sampled in each location were at least 1,000 m apart, the distance between transects was higher than 10 and less than 900 m. The numbers of locations, FIGURE 2 | Scheme illustrating the random nested sampling design on a hierarchy of three spatial scales (transect, meadow, and location) used to study fucalean forests percentage cover and frequency.

Anthropogenic and Biological Stress Index
MA-LUSI is a cumulative index of anthropogenic stress (GIG, 2013; Papathanasiou and Orfanidis, 2018) specific for coastal water benthic macrophytes and inspired by the LUSI index (Flo et al., 2011(Flo et al., , 2019. To calculate the quantitative index, information of key indirect and direct stressors in a 3-km buffer zone around the sampling forests is needed. The ESRI GIS software was used to create the buffer zones on the Corine Land Cover database 2018 1 and to assess the extension of indirect stressors such as urban, agricultural (irrigated land), and industrial pressures in terms of percentage of land cover accounted for by the respective activities. The semi-quantitative information on stressors is classified and assigned in a score (Supplementary Table 1  inputs, Harbor-related impacts (Supplementary Tables 2, 3).
All the scores were summed together and multiplied by correction factors related to hydrology and coastline confinement (Supplementary Table 4) to obtain a numerical value. Herbivory on fucalean forests was estimated by the number of sea urchins found close to the sampled transects.

Statistical Analysis
Frequency indicates the number of times a species was present within a given number of sampled quadrates. It was measured by noting the presence of a species in randomly sampled locations, which are distributed as widely as possible throughout the study area. Since frequency is highly influenced by the size of the quadrats used, we calculated it at a transect scale for which the chosen size was 10 m across the Mediterranean locations sampled. Frequency (%) = (Number of sampling transects in which the species occurs)/(Total number of sampling transects employed for the study) * 100. The percent (%) area of the quadrat covered by a species was used as a measure of abundance as follows: The mean percentage cover = (Total percentage cover of a species in transects in which the species occurs)/(Total number of transects in which the species occur) * 100. The mean percentage cover of each species at each transect was calculated as the mean percentage cover of the species at quadrates sampled in the transect.
To test whether the number of quadrats sampled in each location was sufficient to represent the existing fucalean species variance, the T-S species accumulation curves of Ugland et al. (2003) were used. The species accumulation curve describes the accumulation rates of new species over the sampled area and depends on species identity. Differences in patterns of distribution across spatial scales were tested using Permutational Multivariate Analysis of Variance (PERMANOVA) of nontransformed percentage cover data based on Bray-Curtis dissimilarities, using location (6 levels as random factor), forests (3-7 levels, nested in location, random) and transect (2-12 levels, nested in forests, random) as random variables, n = 3-12 for each transect. Due to the unbalanced design, results were interpreted with a more conservative significance level of α = 0.01 (Underwood, 1996). The locations I_GNI and T_DIDIM weren't included in these analyses since each location hosted a single forest. Multivariate analyses (nMDS, SIMPROF, and SIMPER) were plotted to visualize and explain patterns of dissimilarities at the scale of locations, forests, transects, and spatial relationships among the fucalean species. For all the above analyses, Primer version 7 was used with the add-on package PERMANOVA+. Relationships between percentage cover and anthropogenic (urbanization, agriculture, mariculture, sewage discharge, harbor-related pollution) and natural (confinement, sea urchin grazing) stressors of the MA-LUSI were explored with Redundancy Analysis (RDA) using the CANOCO 5 software. Collinearity of natural and anthropogenic stressors was tested by Spearman rank correlations analysis using R 3.5.0 environment (R Core Team, 2020). All plots were designed using the "ggplot2" package (Wickham, 2016).

Species Frequency and Percentage Cover
At the transect scale, Cystoseira compressa was by far the most frequent species (57.8%) and, therefore, the species with broader distribution in the locations studied (Supplementary Table 5).

Spatial Variability of Species Percentage Cover Across Scales
PERMANOVA of fucalean species percentage cover showed statistically significant differences at all spatial scales (p < 0.001; Table 2). Figure 3 shows the mean percentage cover values of the fucalean species at different spatial scales. Based on FIGURE 3 | Mean percentage cover of eleven fucalean brown seaweed species at the three different spatial scales (transect, meadow, and location) studied. In comparison, the perennial fucalean are indicated by colour columns, the semi-perennials by black and white. the components of variance (Table 2B), the highest variability in percentage cover (51%) was observed at the scale of location, followed by a high variance component associated with residuals (25%). Lower components of variance were observed at forest (15%) and transect (9%) scales. The test of homogeneity of dispersion (PERMDISP) also confirmed that there was a significant difference in within location variance [F (5 , 680) = 62.283; P(perm): 0.001; Supplementary Table 6].
Although locations were originally chosen randomly and so treated as a random factor in the analysis, it was nevertheless of interest to describe and discuss differences among locations. There were differences in species composition of the fucalean assemblages between the locations sampled (Supplementary Table 7). While C. barbatula and C. corniculata were found to be the species with the highest percentage cover in the Greek forests, they were absent in all other locations. In both Italian locations, C. compressa had the highest percentage cover, while C. crinita dominated in Sazan, Albania. The open coasts of Eastern Macedonia and Menorca were the two locations with the highest number (4) of coexisting species.
nMDS analyses on the Bray-Curtis similarity index indicated differences among locations, forests and transects (Figure 4). SIMPROF analysis identified statistically distinct subsets (groups) of species at forest and transect scales. At forest scale seven groups were identified (A-F; Figure 5A). SIMPER analysis results showed species contribution to the dissimilarity between significant SIMPROF groupings based on species-level percentage cover. The species principally responsible for the forests groups were C. corniculata (Group A), G. barbata and C. foeniculacea f. tenuiramosa (B), E. barbatula and C. compressa (C), C. compressa (D), E. amentacea and E. crinita (E), E. crinita (F). At transect scale eight groups were identified (A-H; Figure 5B). The species principally responsible for the transects groups were C. corniculata (Group A), E. mediterranea (B), G. barbata and C. foeniculacea f. tenuiramosa (C), G. barbata and C. compressa subsp. pustulata (D), E. barbatula (E), C. compressa (F), E. amentacea and E. crinita (G), and E. crinita (H).

Spatial Relationship Between the Fucalean Species
nMDS analysis based on fucalean percentage cover at the scale of 10 m transect revealed low species relationship in different transects across the Mediterranean Sea and therefore being heterogeneous, except E. crinita and E. amentacea (Figure 6).

Anthropogenic and Natural Stressors Metric Relationships
The collinearity between the anthropogenic (urbanization, agriculture, mariculture, sewage discharge, harbor-related pollution) and natural (confinement, sea urchin grazing) stressors was weak (ρ < 0.7, Supplementary Table 10). The anthropogenic stressors identified in the studied forests in a decreasing rank were changes in coastline caused by the construction of harbors, by agriculture and by mariculture (Supplementary Table 9). According to MA-LUSI, three forests from the Conero Riviera in Italy (Grotta Azzurra, Passetto, and Passetto-Scalaccia) are featured by the highest anthropogenic stress (MA-LUSI > 4.5). The lowest anthropogenic stress was estimated for the forests of the Albanian coasts (0) and for two Open Eastern Macedonia coasts, Greece (1). The primary natural stress identified for the studied forests was sea urchins grazing, with its highest value estimated at the Grame forest in Albania (3) and the lowest in the Italian, Spanish and Turkish forests (0). Confinement values ranged from 0.75 (convex coastline) to 1.25 (concave coastline). The RDA full model results showed that 72.5% of the response data variance were explained in the first two axes, while 93.5% in the first three axes. As shown in the ordination graph ( Figure 7A) the percentage cover of perennial fucalean species was negatively correlated mainly with mariculture, urbanization and harbors, as well as to MA-LUSI index. The same activities and stress showed a positive correlation with the percentage cover of C. compressa. Mariculture, urbanization, MA-LUSI and harbors also explained statistically (p < 0.05) the fucalean percentage cover variability (%, Table 3).
The RDA forward selection model results showed that only two of the explanatory variables added significantly to the explanatory power of the analysis ( Figure 7B and Table 3). These were mariculture (73.6% to the explanatory power, p = 0.001) and urbanization (that added 11.1%, p = 0.01). The two ordination axes accounted for 61.4% of the total variation in the response data. A very clear pattern emerged in this analysis, with the percentage cover of semi-perennial fucalean species being positively correlated to both explanatory variables, while the perennial fucalean species showed a negative correlation.
Cystoseira compressa showed a different pattern in relation to anthropogenic stress compared to perennial species. This is illustrated in Figure 8, where the polynomial regressions between perennial and semi-perennial fucalean species and MA-LUSI are shown. While the perennial fucalean percentage cover was negatively affected by the increase of MA-LUSI (R 2 = 0.61), the percentage cover of C. compressa and C. compressa var. pustulata increased (R 2 = 0.73). In both cases, changes in percentage cover wasn't linear, with thresholds at values of MA-LUSI greater than about 1.5 and 2.5, respectively. FIGURE 5 | Dendrogram from group-average clustering of the fucalean species percentage cover at forest (A), transect (B) spatial scales. Continuous lines indicate the groups which were significantly differentiated by SIMPROF tests (at the 5% level). Within each of these groups, the null hypothesis that all pairs of species have the same association to each other cannot be rejected, the subgroup structure identified by cluster analysis thus having no statistical support (dashed lines).

DISCUSSION
Even though this is not a manipulative experiment able to disentangle cause effects relationships, this large scale approach allowed to better identify the potential effects of anthropogenic and natural stressors on fucalean species percentage cover in the Mediterranean Sea. Sampling was carried out at eight locations, characterized by extensive rocky shores and more or less accessible from the coast. The forests were selected (see T-S species accumulation curves; Supplementary Figure 1; Ugland et al., 2003) along Mediterranean coasts, from Catalonian coasts in the West to the Aegean Sea in the East (approximately between 4 and 24 • E), and to cover a less extensive latitudinal gradient in the central Mediterranean (approximately between 44 and 40 • N). This approach enabled us to assess distinct spatial patterns in fucalean assemblages and correlate the percentage cover of Cystoseira, Ericaria, and Gongolaria species with key stressors considered the most important for the decline of fucalean forests in the Mediterranean (e.g., Fabbrizzi et al., 2020), either singly or in combination.

Relationship Between Species Percentage Cover and Anthropogenic and Natural Stress
Spatial-scale-based mapping allowed to understand better the effects of anthropogenic and natural stressors on fucalean species percentage cover in the Mediterranean Sea. By examining the effects of spatial scales and main stressors singly or cumulatively on the percentage cover of Cystoseira, Ericaria, and Gongolaria, we could explore by correlative means the stressors that may cause species degradation.
Indeed, the results of this study indicated significant differences in the percentage cover of fucalean species at all spatial scales, with the highest variance detected at the scale of location, i.e., 10 s of kilometers (Figure 4 and Tables 2A,B). Variance at such a broad spatial scale is likely to be related to environmental factors such as local geomorphology, or chronic pollution, e.g., eutrophication (Benedetti-Cecchi et al., 2001;Sales and Ballesteros, 2009;Cefalì et al., 2016). This explanation is also supported by the multivariate analyses results (nMDS, SIMPROF, and SIMPER), where the sampled forests and transects, but not the locations, were sub-grouped based on different species percentage cover (Figures 5A,B). The criteria for this grouping were ecological, e.g., the group of open high hydrodynamic coasts inhabited by C. corniculata or of anthropogenic stressed coasts inhabited by C. compressa, as well as geographical, e.g., groups of a broadly distributed species, G. barbata, and groups of C. foeniculacea f. tenuiramosa or C. compressa subsp. pustulata, with a more restricted distribution. Similar patterns have also been observed in studies of other benthic macrophytes, where the functional metrics like species percentage cover reduced the spatial complexity and showed similarities of habitats with similar ecological conditions (Orfanidis et al., , 2010. In Posidonia oceanica meadows, the highest variations of the standing crop (g dry biomass m −2 ), a comparable metric with percentage cover used in this study, have also been observed at the largest spatial scale (10's kilometers apart). Such a result might reflect differences in the features of habitats in different localities, such as wave exposure, substrate type (rocky vs. pebble), sediment characteristics, and grazing pressure (Balestri et al., 2003). Therefore, the location spatial scalespecific management, also adopted by the European Water Framework Directive (2000/60/EC), i.e., the spatial scale on which the most significant variation at a habitat exists, is relevant to a wide range of present and future conservation and restoration actions.
Fucalean species exhibiting different functional traits, i.e., perennial vs. semi-perennial species, responded differently to anthropogenic and natural stressors. While the existence of mariculture, harbors, and the cumulative anthropogenic stress (MA-LUSI index) negatively correlated with the percentage cover of perennial species, the same stressors seem to produce favorable conditions, at least up to a certain intensity, i.e., MA-LUSI values between 2 and 5, for the growth of C. compressa and C. compressa subsp. pustulata. This result agrees with the well-known pattern of replacement of perennial fucalean species by the more tolerant, relatively fast-growing C. compressa under stressing conditions (Panayotidis et al., 2004;Airoldi and Beck, 2007;Devescovi and Iveša, 2007;Falace et al., 2010;Giakoumi et al., 2012;Kletou et al., 2018). Competitive release is likely to be a critical determinant of fucalean diversity and abundance when certain anthropogenic stressors limit the growth of perennial species, allowing the semiperennials to expand and dominate in the community (Segre et al., 2016). However, although the specific mechanisms behind these changes have not been fully understood yet, habitat destruction and decrease in water quality are likely to play a major role in the decline of perennial species (Tsiamis et al., 2013;Thibaut et al., 2015;Iveša et al., 2016;Rindi et al., 2020). These processes are indicators of nutrient enhancement, water turbidity, and high eutrophication levels, which are invoked in several studies as the main causes for the regression of fucalean species Mancuso et al., 2018) and seagrasses  in the Mediterranean Sea. The input of nutrients and changes in water transparency are considered among the processes affecting the growth of macrophyte communities (De Jonge et al., 2002;Viaroli et al., 2008). In a recent paper (Fabbrizzi et al., 2020), geomorphological features were recognized among the most relevant drivers predicting presence of fucalean seaweeds, followed by anthropogenic variables such as distance from ports and urbanization.

Species Frequency and Percentage Cover
At the transect scale, Cystoseira compressa was by far the most frequent species. Variation in distribution may be caused by several factors like growth pattern, amount and dispersal of zygotes, and grazing (Falace et al., 2005;Mangialajo et al., 2012). In terms of growth pattern, fucalean seaweeds in the Mediterranean generally undergo a morphological shift from the period of main growth in late winter-early summer to dormancy in late summer-autumn, when many species shed a large part of their fronds (Orfanidis et al., 2017). Indeed, fronds of C. compressa undergo senescence and get detached in summer, after the alga has released the gametes (Sauvageau, 1912). In late summer and autumn this species usually consists of only a small, perennial holdfast and a few short flattened branches (the so-called rosetta form, Cormaci et al., 2012), which will issue new fronds in the subsequent winter. Zygotes of Cystoseira compressa tend to adhere to parental receptacles, remaining entrapped in a layer of mucilage formed on the surface of the alga (Sauvageau, 1912;personal observation). We suggest that this may be a strategy of C. compressa to resist strong hydrodynamic and unsuitable habitat conditions or to expand distribution by fragments of fronds floating far away from parental populations as reported for Sargassum muticum (Yendo) Fensholt (Deysher and Norton, 1981).
Knowing the abundance patterns of different species can provide insight into how a community or ecosystem functions and how the processes link the local abundance of a species and its regional distribution (Brown, 1984). The present study indicated a large variation locally in the mean percentage cover of fucalean seaweeds, confirming that assemblages formed by these algae are fragmented and heterogeneous across the Mediterranean Sea. Differences in species composition among the coasts have been suggested as a key biological feature of Mediterranean biogeography (Coll et al., 2010;Sales et al., 2012). However, since fucalean seaweeds unequivocally dominated communities characterized by good water quality with low anthropogenic stress, we argue that a severe local decline may be caused mainly by habitat destruction, decrease in water quality, and overgrazing by herbivores (Tsiamis et al., 2013;Thibaut et al., 2015;Iveša et al., 2016;de Caralt et al., 2020). This leads to replacement by relatively fast-growing species (e.g., C. compressa) or a shift to less-structured assemblages formed by morphologically simple algae (i.e., turf-forming, or other ephemeral, opportunistic seaweeds), mussels, or barren grounds (Airoldi and Beck, 2007;  3 | Simple term effect of explanatory variables of RDA full and forward selection models between the percentage cover of perennial fucalean seaweeds and semi-perennial C. compressa and C. compressa var. pustulata found in the sampled locations, and the main anthropogenic (urbanization, agriculture, mariculture, sewage discharge, harbour-related pollution) and natural (confinement, sea urchin grazing) stressors (for more information see Figure 6). FIGURE 8 | Polynomial regression between MALUSI and (A) percentage cover of perennial fucalean species (y = 0.73x 3 − 6.81x 2 + 7.39x + 41.96) and (B) the semi-perennial species Cystoseira compressa, C. compressa var. pustulata (y = 1.24x 3 − 3.93x 2 + 0.4x + 9.52). Devescovi and Iveša, 2007;Falace et al., 2010;Giakoumi et al., 2012;Sala et al., 2012;Kletou et al., 2018).

CONCLUSION
The results of this study on the genera of Cystoseira, Ericaria, and Gongolaria in the Mediterranean Sea provided new insights into: (a) the role of different anthropogenic and natural stressors, which can individually or cumulatively affect these algal forests, (b) the differential responses of species belonging in different trait (functional) groups, i.e., perennial vs. semi-perennial species. However, experimental studies are additionally required to mechanistically identify the drivers for the observed replacement of perennial fucalean species by relatively fast-growing semiperennial species, e.g., C. compressa, or the general seaweed regression along the Mediterranean Sea. Understanding the effect of multiple stressors is particularly challenging because their potential cumulative effects on these habitats cannot be predicted in a single-stressor framework.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.