Habitat-complexity Regulates the Intensity of Facilitation along an Environmental Stress Gradient

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is the peer reviewed version of the following article Navarro-Barranco, Carlos & Gribben, Paul & J, Janine & Poore, Alistair. (2022). Habitat ‐ complexity regulates the intensity of facilitation along an environmental stress gradient. Oikos., which has been published in final form at https://doi.org/ 10.1111/oik.08818. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Use of Self-Archived Versions. This article may not be enhanced, enriched or otherwise transformed into a derivative work, without express permission from Wiley or by statutory rights under applicable legislation. Copyright notices must not be removed, obscured or modified. The article must be linked to Wiley’s version of record on Wiley Online Library and any embedding, framing or otherwise making available the article or pages thereof by third parties from platforms, services and websites other than Wiley Online Library must be prohibited." Habitat complexity regulates the intensity of facilitation along an environmental stress gradient Navarro-Barranco C1, Gribben PE2, Ledet J2, Poore AGB2 1 Laboratorio de Biología Marina, Departamento de Zoología, Universidad de Sevilla, Spain 2 School of Biological, Earth and Environmental Sciences, University of New South Wales, Sydney, Australia * Corresponding author: [email protected]s ABSTRACT 1. Positive interactions between ecosystem engineers and their associated species are expected to be influenced by the degree of environmental stress as well as trait variations of the species involved. However, there is scarce empirical evidence about the way in which these two factors interact and shape the intensity of facilitation processes. With this aim, a colonization experiment using artificial units with different variations in a functional trait (morphological complexity) have been conducted at different intertidal height levels. 2. Our results point out how the interplay between the environmental context and facilitator traits determine the intensity of facilitation processes. Complexity of host species play a key role in determining the structure of its associated epifauna, but its relevance differed along environmental gradients of the intertidal zone. 3. The effect of habitat complexity for mobile epifauna can be detected since early stages of ecological succession. Increased shore height level required higher host complexity to maintain stress levels within tolerance threshold of associated species. 4. The isolated effects of shore-height and host-complexity was in accordance to expectations: number of species and abundance values decrease as units were less complex and remains less time underwater. However, in contrast to the initially proposed model, based on Stress Gradient Hypothesis, magnitude of positive interactions increased at lower shore heights (those levels with, presumably, milder environmental conditions). 5. An additional predation exclusion experiment allows us to rule out the effects of predation pressure as a major factor driving the structure of mobile epifaunal assemblages on intertidal habitats. Thus, further studies are needed in order to show up the explaining factors and processes behind the pattern of facilitation intensity found in the present study. Keywords: Facilitation, Stress Gradient Hypothesis, Foundation species, Complexity, Epifaunal assemblages, Amphipoda. INTRODUCTION Both sign and magnitude of species interactions are context dependent, which implies that the relative importance of processes such as competition, predation and facilitation may change along biotic and/or abiotic gradients (Crain & Bertness 2006, Lany et al. 2017). Historically, the study of negative interactions (e.g. competition and predation) have dominated ecology while the interest about the role of positive species interactions as key structuring forces in biotic communities have arisen during the last decades (Bruno et al 2003). Along this period, most of facilitation research has focused on “Stress Gradient Hypothesis” (SGH), a theoretical model formulated by Bertness & Callaway (1994) which predicts an increase in the frequency of facilitation at increasing level of environmental stress. As a result of many attempts in order to test the predictions of this model, consistent evidences have been found in a wide range of terrestrial and aquatic environments, proving that positive effects are common under abiotically stressful conditions but weak or absent under benign conditions (Bulleri 2009, Maestre et al. 2009, He et al. 2013). The role of foundation species and/or ecosystem engineers such as trees, corals, seagrasses or mussels are common examples of positive interactions. By creating extensive habitats, providing competitor- or predator-free space and resources, or ameliorating environmental conditions, those species enhance the fitness and distribution of other species (Jones 1994, Crain & Bertness 2006, Bruno et al. 2003). When interacting species have similar means of resource acquisition (e.g. plan-plant interactions competing for nutrients or marine sessile species competing for space), a shift in the direction of the interaction have often been observed along stress environmental gradients (Callaway et al. 2002, Maestre 2009, Bulleri et al. 2011). However, ecosystem engineers species and their associated fauna often belong to different trophic levels and, in such cases, direct negative effects trough competition are negligible or absent. In this context, SGH would predict an increase in the magnitude of the facilitator effects associated to higher environmental stress levels (Bruno & Bertness 2001). The degree to which a foundation or ecosystem engineer species reduce the biotic or abiotic stress relies on host traits such as size, density, health or architectural complexity. So, it is expected facilitation to be not only influenced by the environmental context but also trait-dependent (Bruno & Bartness 2001, Bateman & Bishop 2017). In this sense, related literature provide many examples: Morphological traits of mangrove forest influence both the establishment and persistence of facilitation interactions (Bishop et a. 2013); Decline of fish populations occur when coral biomass and complexity are reduced (Syms & Jones 2000); Density and height of cordgrass shoots modify the ability to facilitate other species (Irving & Bertness 2009); Slight changes in the morphology of alpine plants are responsible of facilitation intensity discrepancies between previous field studies and SGH predictions (Anthelme et al. 2017). However, little is known about how the importance of trait variation in ecosystem engineers is related with environmental stress and the influence of such interactions in facilitation processes. Bruno and Bertness (2001) proposed a conceptual model of the hypothesized relationship between facilitator traits and species diversity within an environmental stress gradient. When ambient stress is negligible, species diversity do not depend of the facilitator capacity of reducing local stress and changes in facilitator traits do not increase or reduce the magnitude of facilitation effects. In contrast, increased ambient stress may amplify differences in the relevance of the habitat quality provided by ecosystem engineers. Thus, it would be expected that facilitator trait variations determine higher differences in community response (diversity, abundance, etc.) when the ambient stress increased (Figure 1). In the current study, we aim to test those predictions using a common facilitation relationship (sessile ecosystem engineers and their associate community) within a well-known environmental gradient (biotic and abiotic stress variations from high intertidal to shallow subtidal habitats) and a relevant facilitator trait (the role of morphological host complexity). Numerous studies have provided support for the stress gradient hypothesis in a variety of shallow subtidal habitats (e.g. Bulleri et. 2011; Bennett et al. 2015). Nevertheless, the role of macroalgae and sessile invertebrates (e.g. mussels, barnacles) ameliorating thermal and desiccation stress at intertidal communities is, by far, the most studied facilitation interactions in marine environments. Macroalgal canopies expanded the upper species borders of many species and increase the overall species richness and diversity (Bertness & Leonard 1997, Thornber et al. 2017 and references therein). In agreement with the predictions of the SGH, previous studies revealed a greater importance of the canopy (either in terms of presence/absence or canopy density) at higher intertidal elevations (Bertness et al. 1999, Watt & Scrosati 2013a, 2013b). However, variation in predation pressure should be taken into account: It may not follow the same pattern than abiotic stress gradient since intertidal habitat-former species can both increase or decrease predation pressure by acting as refuge for preys or increasing predator abundance (Bertness et al. 1999, Bateman & Bishop 2017). It has been proposed that, while the amelioration of physical conditions by ecosystem engineers is essential at high elevations, those species enhance biodiversity at subtidal areas mainly by acting as competitor or predator refuge providers (Crain & Bertnesss 2006). Nevertheless, some studies suggested that the role of predation pressure as a controlling factor for mobile epifaunal communities on subtidal rocky reef habitats is very limited (Taylor 1998, Leclerc & Viard 2018). Thus, despite the opposite trends showed by predation pressure and physical stressful conditions, it is expected that there would be the latest factor the main responsible of the environmental stress gradient along the spatial scale considered in this study. Despite all that, interactions strength in the intertidal zone seems to be more related with organismal traits (i.e. body size, trophic level) than environmental conditions associated with tidal height (Wood et al. 2010). Morphological complexity of ecosystem engineers is widely recognized as a key factor shaping the structure of the community (Edgar & Robinson 1992, Hauser et al. 2006, Kovalenko et al. 2012 and references therein, Veiga et al. 2014). Nevertheless, other articles argue that the influence of complexity is limited, suggesting than spatial variations in associated epifaunal are mostly affected by other factors (e.g. epiphytic load, chemical composition, adjacent species) rather than by differences in the structural complexity of the facilitator species (Bologna & Heck 1999, Garcia-Sanz et al 2012, Gartner et al. 2013). Finally, responses among different taxonomic groups could also be different. For example, while some authors highlighted the suitability of habitat complexity as a good predictor of the abundance of small crustaceans (i.e. amphipods, isopods, copepods, etc) (Taylor & Cole 1994, Cacabelos et al. 2010) other studies found no differences in abundance between algae of different complexity (Russo 1990, Schreider et al. 2003). Thus, there is no consensus regarding how the role of complexity of habitat forming species could vary among different habitats and taxonomical groups. The main objective of this study is to determine how variations in functional traits modulate the intensity of facilitation interactions along an environmental stress gradient. We hypothesized that the importance of such variations (i.e. morphological complexity in a marine ecosystem engineer) would be greater in the harsh physical conditions of the high intertidal, rather than the relatively benign conditions of subtidal environments. MATERIAL AND METHODS Study site and Artificial units The experiment was conducted at Chowder Bay (Sydney, Australia) (33°50’22’’S, 151°15’16’’) between 9th November 2017 and 3th December 2017. Chowder Bay is an embayment close to the mouth of Sydney Harbour and with relatively low urban development (land use in the area comprise mainly national park and parklands), which ensure a relatively good water renovation and environmental quality. In order to explore the role of complexity in the present study, artificial units made of rope with different complexity were created. Although it is expected that artificial units do not exactly resemble those of natural assemblages, the use of artificial structures has been proved as a valuable tool for monitoring motile fauna in marine environments (Gee & Warwick 1996, Cacabelos et al. 2010). They provide homogeneous habitat and are easy to manipulate, minimizing the level of variation between replicates and enabling researchers to disentangle the role of specific factors (e.g. complexity, size, time of colonization, environmental conditions) affecting organisms’ distribution (Edgard & Klumpp 2003, Hansen et al. 2010). Each artificial unit were made of 50cm of braided polypropylene rope divided into seven “different branches” joined at the base by a cable tie. Thus, volume was standardized for all complexity treatments. Eight different complexity levels were stablished by increasing the number of frayed ropes; from none (lowest complexity, level 1) to all “branches” frayed (highest complexity, level 8). The suitability of such artificial units made of frayed rope for ecological studies of marine mobile epifauna communities have been previously highlighted (O’Connor & Bruno 2009, Lefcheck et al. 2016). Measures of area, perimeter and fractal indices of both area (Da) and perimeter (Dp) were calculated for each artificial unit by image analyses (ImageJ program, Schneider et al. 2012) following the methodology described by McAbendoth et al. (2005). Fractal indices have been proved to be useful parameters to characterize morphological complexity in both natural and artificial substrates and their influence on epifaunal assemblages (Jeffries 1993; Veiga et al. 2014; Suárez-Jiménez et al. 2017). While Da is an estimate of the way the host substrate stands are divided up (higher values are associated with more and smaller divisions), Dp is related with the degree of convolution of its edge (McAbendroth et al. 2005). Experimental design and sampling procedure Four different shore heights were selected: “Low”, “Medium-Low”, “Medium-High” and “High”, separated by 30cm each. “Low” level was established as a very shallow but permanently subtidal level during the study period, according to tide tables for Sydney (Fort Denison) provided by Bureau of Meteorology (Australian Government). Artificial units were hanged from the pier, so that four replicates of each complexity treatment (1-8) were deployed at each height level. Artificial units were deployed at random locations but with enough separation to avoid contact among them or with pier pillars. After 24 days on the field, they were collected during a strong coefficient high tide, when all units remain underwater, and preserved in 70% ethanol. In the laboratory, samples were washed through a 0.5mm mesh sieve to collect all mobile epifauna. All organisms were sorted into higher taxa (phylum) and all amphipod specimens were identified to lowest taxonomic level possible (species or genus) using binocular microscopes and counted. Additionally, another experiment was conducted in order to explore the relevance of predation pressure in the studied epifaunal community. The impacts of predatory fishes and decapods in epifaunal assemblages are most commonly manipulated using exclusion cages (Sala et al. 1997; Hindell et al. 2001; Martinetto et al. 2005). Three different treatments were established; exclusion cages (artificial units enclosed within six-wall cages), experiment control (replicates without cages; hereafter “uncaged samples”) and cage effects control. Unlike exclusion cages (made of six 25cm2 walls with 1 cm mesh size), cage controls showed big holes of 15cm2 in each lateral wall. By this way, cage control mimic some of the physical effects of cages (e.g. decrease of insolation, sediment retention) without exclusion effect. This methodology was also followed by Machado et al. (2019) in order to investigate the effect of predation on an assemblage of associated amphipods. Only two shore heights (“Low” and “Medium-High”) and two complexity levels (either 2 or 6 out of 8 branches were frayed) were selected for this experiment. Four replicates of each combination of factors (predation treatment, shore height and complexity) were deployed below the pier for 29 days (from 5th February to 6th March 2018). After this period, sampling collection and processing was conducted in the same way that previous experiment. Data analyses All morphological parameters (Area, Da and Dp) were contrasted with the predictor variable of complexity treatment using a linear model. Multivariate analyses were conducted, using the abundance of amphipod species per replicate, to look for patterns in community structure of the mobile epifauna. Data were previously squareroot transformed to reduce the importance of the numerically dominant taxa. Bray-Curtis distances were calculated and visualized by multimensional scaling (MDS) using the function metaMDS in the R package vegan (Oksanen et al. 2011). Permutational multivariate analyses of variance (Anderson 2001) was used to explore whether there were differences in the epifauna community associated with factors “Height” or “Complexity” (Adonis test with 9999 permutations). Based in the results previously obtained with morphological parameters considered, complexity was considered as a continuous variable. Further differences between pairs of height levels were explored by pairwise analyses applying Bonferroni method to correct the P-values. Betadisper test was applied to test differences in group homogeneities. Both abundance of amphipod individuals and species richness values were contrasted with the predictor variables of height above sea level and complexity treatments using a negative binomial generalized linear model (Abundance data) and a linear model (Species richness). Assumptions of all tests were checked by plots of the residuals versus fitted values. When categorical factor “Height” showed significant effects, differences among levels within this factor were tested using Tukey’s post hoc test. Significant Height byComplexity interaction were further explore by testing differences of slope of the continuous variable for each pairs of levels of the categorical variable. Univariate statistical analyses were also conducted in the R software environment (R Development Core Team 2011) using the associated package MASS (Venables and Ripley 2002). Finally, the influence of predation pressure on the number of amphipod individuals and species was tested through a three-way Analysis of Variance (ANOVA) conducted with the following factors; “Predation treatment” (with three levels; cages, cage control and uncaged samples), “Shore Height” (with two levels; high and low) and “Complexity” (with two levels; simple and complex). Since only two complexity levels were considered in this experiment, complexity was considered as a fixed factor, similarly to the other factors included in this study design. All factors were orthogonal among them. Prior to analyses, the homogeneity of variance was checked using Cochran’s test (Underwood 1997). Analyses were conducted with the GMAV-5 program (Underwood et al. 2002). RESULTS Influence of complexity and shore height on facilitation intensity All morphological characteristics measured for our artificial units showed a significantly positive lineal relationship with complexity treatments, which support the use of our complexity treatment gradient as a continuous variable (Figure 2). A total of 10.614 individuals of mobile epifauna were collected at the end of the study. Arthropoda was the dominant phylum, accounting 81% of all specimens (99.8% of all arthropods were amphipods), followed by Mollusca (14%) and Nematoda (3%). Eighteen different amphipod taxa were found, being Jassa slatteryi, Caprella aff. penantis, Stenothoe sp., Ericthonius pugnax and Apocorophium acutum the dominant species at all heights considered (Table 1). Despite this, MDS analysis showed a clear gradient in the structure of the epifaunal community following the height and complexity gradient considered in the study (Figure 3). In this sense, the epifaunal community inhabiting artificial units of low complexity in a certain height level resemble those of high complexity deployed in the height level below. These patterns were supported by PERMANOVA analyses, which showed that Height, Complexity and the interaction between them have a significant effect on the mobile epifaunal assemblages (Table 2). Pair-wise comparisons made between heights levels indicate that all levels were significantly different from each other with the exception of MH and H. Non-significant result in betadisper indicate that those differences among heights are due to differences in composition between groups and not by heterogeneity of variances. Univariate analyses conducted for both species richness and abundances of individuals showed significant differences between height levels, as well as a significant positive effect of complexity (Table 3 and 4). Interaction between both factors indicate differences of the slope of the continuous variable (complexity) for some levels of the categorical variable (height). Regarding species richness, “Low” artificial units showed significant higher values than those deployed at all other heights and the slope of their response to different complexity values is also significantly higher (Figure 4a, Table 3). Epifaunal communities at “Medium-Low” units also showed higher number of species than those collected at higher heights but their slope, although steeper, was not significantly different. Finally, there were no differences between “Medium-High” and “High” levels, neither in their mean values or their absence of response to complexity. Post hoc analyses for abundance values provided similar results: “Low” units showed the highest values, followed by “Medium-Low”, “Medium-High” and “High” levels, respectively (Figure 4b). All comparisons were significantly different with the exception of “Medium-High” and “High” abundance values (Table 4). However, in this case, response curve to complexity at “Medium-Low” level showed the steepest slope (significantly steeper than those obtained at “Low” and “High” levels). Complexity also showed a significant effect on abundance values at “Low” and “Medium-High” levels, but the slope of their response curve was not significantly different than those obtained for “High” levels. Predation exclusion experiment All colonizers amphipods found in the previous experiment were also found here, with the addition of two caprellids; Caprella scaura and Caprella danilevskii (Table 1). Overall, dominant species were similar to previous experiment, comprising J. slatteryi, C. aff. penantis and A. acutum more than 70% of all individuals. Exclusion of predators do not increase neither the abundance nor the species richness of amphipods (Figure 5). On the contrary, mean number of individuals and species was higher on “uncaged samples” than “cages”, although those differences were only significant for number of species, according to ANOVAs (Table 5). While height was not a significant factor for abundance data, lower artificial units showed significantly higher number of species than those deployed at higher level. Finally, complexity was a significant factor for both abundance and species richness; higher values were found in artificial units with higher complexity. DISCUSSION Spatial patterns (difference between heights) in abundance, species richness and community composition were in accordance with expectations. Changes in community structure and decrease in both abundance and species richness values as units remains less time underwater showed the evident gradient of environmental harshness. On the other side, positive effects of increasing complexity values (even when response variables were standardized by surface area) reflects the relevance of this trait in facilitation relationships in shallow and intertidal environments. Both multivariate and univariate analyses showed how, for a certain level, higher values of complexity are needed in order reach similar values than those observed at lower heights. This support one of the predictions by Bruno & Bertness (2003) regarding the importance of facilitator traits under ambient stress gradient: Increased ambient stress would require further change in the trait of interest to maintain stress levels within tolerance threshold of their associate species. They also hypothesized than response curve of the associate species to facilitator traits may adopt different shapes and include both lower and upper thresholds. In this case, complexity showed a simple and linear relationship with species richness values while its influence on abundance values decrease at lower complexities. A priori, absence of upper thresholds values seems likely (specially for abundance values) since increasing complexity is always linked to greater surface area. The overwhelming dominance of Amphipoda have been also highlighted in previous colonization studies using artificial collectors with varying structural complexity (Garcia-Sanz et al. 2012; Lavender et al. 2017). Most abundant taxa (such as genus Jassa, Caprella, Stenothoe Apocorophium or Ericthonius) are successful early colonizers of artificial structures and, thus, common inhabitants of fouling communities associated to aquaculture facilities, boats or floating pontoons (Guerra-Garcia et al. 2015; Navarro-Barranco et al. 2015; Fernandez-Gonzalez & Sánchez-Jerez 2017). Some authors proposed that the effect of habitat complexity is less profound at early stages of ecological succession, since many biological processes modulated by complexity (e.g. niche segregation, protection) require more time to be relevant (Hansen et al. 2010, Torres et al. 2015). In this sense, our results point out how a relatively short period of time (approximately 1 month) is enough to prove the relevance of host complexity for mobile epifauna community. According to our results, relevance of complexity in the development of epifaunal communities is higher at lower height levels, reaching its maximum relative importance at “Low” level for species richness and “Medium-Low” levels for abundance values. This results contrast with our initial expectations; Facilitation effects would be intensified alongside with the increasing necessity of amelioration of temperature and desiccation stress toward higher intertidal elevations. Result obtained by Bertness et al. (1999) exploring the effect of intertidal algal canopy on their associate community along with environmental conditions supported that hypothesis. Similarly, Watt & Scrosati (2013b) found that effects of seaweed coverage affected benthic species richness, diversity and composition at the high and middle intertidal zones but not at the low zone. However, posterior studies conducted in the same area revealed that canopy effects are greater at middle than high elevations (Scrosati 2017), suggesting that facilitation intensity peaks at an intermediate level of environmental stress. All our analyses also reveal a lack of complexity influence at the highest intertidal level considered in the study. One of the proposed explanation for facilitation decreasing at high stress is related with changes in some relevant traits of the facilitator (e.j. differences in seaweed morphology at high elevations) (Scrosati et al. 2017). However, this explanation is not applicable to our data, since the use of artificial units ensures the homogeneity of samples among heights. Drop in facilitation intensity at higher elevation levels could be explained by the inability of facilitator to improve the conditions strongly enough at such elevations (i.e. even the most complex artificial units do not ameliorate desiccation stress). Several studies highlighted exceptions to SGH under weak stress gradients or, specially, at the extreme ends of the stress gradient (when the stress range surveyed is bigger than specie’s niches) (Malkinson & Tielbörger 2010; He & Bertness 2014). On the other side, reasons behind the higher community response to complexity values at subtidal or low intertidal levels remains more obscure. Presumed decrease of complexity effects at subtidal level was related to mild conditions (i.e. if there are not harsh conditions to ameliorate, there is no positive effect of increasing complexity). Bateman & Bishop (2017) found greater effects of mussels beds on their associated species at subtidal than intertidal elevations, which could be driven by the higher predator pressure at lower elevations. Under relatively low physical stressors effects, competition and predation become more relevant and epifaunal organisms would associated with habitat forming species that provide new niches and predation refugee. Same explanation was used by Umanzor et al. (2017) to explain the significant effect of intertidal algae as ecosystem engineer even when environmental conditions are milder. However, our predation exclusion experiment does not support the hypothesis of a greater predation pressure at lower elevations; those samples deployed inside predation exclusion cages do not show higher abundance and/or species richness than “uncaged samples”, neither on subtidal or intertidal levels. This is in accordance with previous studies, suggesting that predation has a minor role driving the structure of epifaunal assemblages (Russo 1991; Machado et al. 2019). Although predation pressure by fishes or seabirds on intertidal amphipods may be more intense in certain habitats (e.g. rockpools) or seasons (Nielsen 2001; Hilton et al. 2002; Compaire et al. 2016), there is no evidence of a noteworthy depth gradient in predation pressure on littoral amphipods. The higher number of species in our uncaged samples suggest that any benefit provided by the cages in terms of predation protection is relatively unimportant in comparison to other factors (e.g. cages may hinder the colonization of artificial units by amphipods dispersed by rafting). Seeking for shelter from dislodging wave shock has been proposed as another relevant factor explaining amphipods’ host preference (Sotka 2007). Amphipod density patterns among different host types or along depth gradients may be enhanced by higher wave exposure (Korpinen & Westerbom 2010). Nevertheless, this is a larger scale effect (manifested as an association to different shore exposures) and its influence at small spatial scale such as those considered in our study system seems less likely. Therefore, it is unclear why epifaunal response to increasing complexity levels peaks at lower heights of the intertidal zone, contrary to our initial expectations Figure 3 Figure 4 Figure 5