Context‐dependent functional dispersion across similar ranges of trait space covered by intertidal rocky shore communities

Abstract Functional diversity is intimately linked with community assembly processes, but its large‐scale patterns of variation are often not well understood. Here, we investigated the spatiotemporal changes in multiple trait dimensions (“trait space”) along vertical intertidal environmental stress gradients and across a landscape scale. We predicted that the range of the trait space covered by local assemblages (i.e., functional richness) and the dispersion in trait abundances (i.e., functional dispersion) should increase from high‐ to low‐intertidal elevations, due to the decreasing influence of environmental filtering. The abundance of macrobenthic algae and invertebrates was estimated at four rocky shores spanning ca. 200 km of the coast over a 36‐month period. Functional richness and dispersion were contrasted against matrix‐swap models to remove any confounding effect of species richness on functional diversity. Random‐slope models showed that functional richness and dispersion significantly increased from high‐ to low‐intertidal heights, demonstrating that under harsh environmental conditions, the assemblages comprised similar abundances of functionally similar species (i.e., trait convergence), while that under milder conditions, the assemblages encompassed differing abundances of functionally dissimilar species (i.e., trait divergence). According to the Akaike information criteria, the relationship between local environmental stress and functional richness was persistent across sites and sampling times, while functional dispersion varied significantly. Environmental filtering therefore has persistent effects on the range of trait space covered by these assemblages, but context‐dependent effects on the abundances of trait combinations within such range. Our results further suggest that natural and/or anthropogenic factors might have significant effects on the relative abundance of functional traits, despite that no trait addition or extinction is detected.


| INTRODUCTION
The study of functional diversity can reveal the mechanisms that underpin the assembly of communities (Hooper et al., 2005;Mouchet, Villéger, Mason, & Mouillot, 2010). Functional diversity expresses the diversity of functional traits, which are those components of an organism's phenotype that influence its response to environmental stress and its effects on ecosystem properties (Hooper et al., 2005;Petchey & Gaston, 2006). Accordingly, processes related to niche differences among species, such as niche complementarity, are intimately linked with variations in functional diversity (Mason, de Bello, Mouillot, Pavoine, & Dray, 2013). For instance, steep environmental gradients, like in wetlands and intertidal rocky shores (Berthelsen, Hewitt, & Taylor, 2015;Mudrak et al., 2016), can drive spatial shifts in composition from functionally similar species due to the selection of traits adapted to the local abiotic and biotic conditions (i.e., trait convergence), to functionally dissimilar species that occupy different niches (i.e., trait divergence; Boersma et al., 2016;de Bello, Leps, & Sebastia, 2006). Functional divergence, in turn, allows the coexistence of functionally dissimilar organisms through niche complementarity (Mason, de Bello, et al., 2013). Assessing how different facets of functional diversity vary from convergence to divergence and the spatiotemporal persistence of these patterns could improve our ability to unravel the mechanisms underlying community assembly.
Ecosystems characterized by strong environmental stress gradients provide an opportunity to understand how communities assemble through the spatiotemporal variation in functional diversity.
Environmental filtering occurs when environmental stress exceeds the physiological limits of tolerances of individuals in local populations (Somero, 2010), which are selected for subsets of organisms with similar functional traits that allow them to persist in a given area (Leibold et al., 2004). Such trait convergence can be represented by reduced functional diversity within the assemblage (e.g., Fukami, Bezemer, Mortimer, & van der Putten, 2005). As environmental stress decreases, reduced similarity and increased niche differences among competing individuals allow coexistence, leading to a stronger influence of niche complementarity on community assembly (MacArthur & Levins, 1967;Mason, de Bello, et al., 2013). This functional divergence leads to increased functional diversity during community assembly (e.g., Mudrak et al., 2016). As a consequence, species-rich communities should show reduced functional diversity in areas where environmental stress is high compared to areas where stress is low.
Functional diversity is a multifaceted concept that can be represented on a multidimensional space where the axes are functional traits (recently reviewed in Carmona, de Bello, Mason, & Leps, 2016).
In such trait space, functional richness (Villéger, Mason, & Mouillot, 2008) represents the range occupied by a given community and can be estimated as the number of functional trait combinations.
Functional dispersion (Laliberté & Legendre, 2010) is estimated as the mean distance of all species to the weighted centroid of the community in the trait space, being equivalent to the multivariate dispersion (Anderson, Ellingsen, & McArdle, 2006). When the influence of species richness on these indices is statistically removed (Mason, Lanoiselee, Mouillot, Wilson, & Argillier, 2008), functional richness and dispersion are "pure" estimators of the occupied range of functional space and the dispersion in trait combination abundances, respectively (Mason et al., 2012). Interestingly, functional richness and functional dispersion can show differing spatiotemporal patterns of variation. Environmental stress would not be strong enough to filter out functional roles, but it would be enough to generate significant differences in the number of individuals within functional trait combinations (Boersma et al., 2016). In temperate regions, for instance, the seasonal variation in environmental conditions and recruitment intensity leads to significant changes in the abundance of marine seaweeds (Cavanaugh, Siegel, Reed, & Dennison, 2011;Valdivia, Gollety, Migne, Davoult, & Molis, 2012). Then, these environmentally driven effects would result in a significant variation in multivariate trait dispersion (i.e., varying functional dispersion) despite that all trait combinations are present under the full range of environmental conditions (i.e., similar functional richness).
Intertidal rocky shore communities represent an ideal habitat to investigate how functional diversity varies along environmental stress gradients, because of the strong vertical gradients observed within a few meters for abiotic factors such as temperature, desiccation, irradiance, and osmotic potential (Gómez & Huovinen, 2011;Menge & Branch, 2001;Raffaelli & Hawkins, 1996;Stephenson & Stephenson, 1949). In these systems, competition for settlement space increases as the abiotic environmental stress decreases from high to low vertical elevations, and biotic interactions such as predation and grazing prevent competitive exclusion in the low-intertidal zone (e.g., Connell, 1961;Donahue, Desharnais, Robles, & Arriola, 2011;Paine, 1974). Then, niche complementarity should take a larger role in determining species occurrences and abundances, allowing more species to coexist in the less stressful environment (MacArthur & Levins, 1967). Regarding the horizontal variation in community structure, local and mesoscale factors such as productivity, wave exposure, and recruitment variability modulate species abundances and occurrences to lesser extent (Navarrete, Wieters, Broitman, & Castilla, 2005;Roughgarden, Gaines, & Possingham, 1988;Valdivia, Aguilera, Navarrete, & Broitman, 2015).
For instance, high wave-exposure can hamper the grazing activity of low-intertidal gastropods and sea urchins, allowing for competitive exclusion by algae (e.g., Hawkins & Hartnoll, 1983;Underwood & Jernakoff, 1981. Environmental filters and biotic interactions have therefore interdependent selective effects on trait combinations and thus on functional diversity. Despite the large body of literature arising from the study of intertidal rocky shore communities, less work has been carried out on how different facets of functional diversity vary at different spatiotemporal scales in these ecosystems (but see Berthelsen et al., 2015).
In this work, we analyzed a dataset of 68 intertidal macrobenthic species (macroalgae and both mobile and sessile invertebrates) to test the hypothesis that functional richness and dispersion increase from high-to low-intertidal elevations, because of decreasing effects of environmental filtering on trait combination' occurrences and abundances along the vertical stress gradient. We further assessed the persistence of local patterns of functional richness and dispersion over multiple sites (covering ca. 200 km) and times (seasonal samplings over three years) in a region showing a mesoscale structure in ecological subsidies such as chlorophyll-a concentration and invertebrate recruitment (e.g., Van Holt et al., 2012).

| Studyregionandenvironmentalconditions
Between June 2012 and June 2015, four intertidal shores were sampled every 6 months along ca. 200 km of the shore (Figure 1; 39.5°-40.5°S). The shores were Cheuque, Calfuco, Chaihuín, and Pucatrihue ( Figure 1). The sites are representative of wave-exposed southeastern Pacific shores in a region characterized by a significant mesoscale variation (i.e., between 10s and 100s of km) in chlorophyll-a concentration, related to freshwater inputs and apparently to runoff from forestry and agricultural activities (Van Holt et al., 2012). In addition, a seasonal upwelling center is located in Punta Galera (Letelier, Pizarro, & Nunez, 2009) Satellite-derived MODIS data were downloaded from http://oceancolor.gsfc.nasa.gov/. Air temperature data would have been desirable to better characterize the environmental conditions to which intertidal organisms are exposed. However, air temperature data for southeastern Pacific shores are available at too coarse resolution (1°), which would have not fitted to our study.

| Sampling
At each site and sampling time, we estimated species abundances on ten 0.25 m 2 plots located along ca. 20-m alongshore transects.
Transects were replicated at three intertidal elevations, low-, mid-, and high-intertidal zones. On each bench, we used perennial marine species occurring highest on the shore as indicators of the upper intertidal boundary. Studies on wave-exposed shores conducted elsewhere show that their upper distribution limit on different sites represents a summary of the local wave regime across multiple sites (Harley & Helmuth, 2003). The barnacle J. cirratus is the sessile and perennial species occurring highest on these shores. Once the upper boundary was determined on each shore, we divided the intertidal range into three zones of roughly equal vertical extent (i.e., high, mid, and low zones). Plots were haphazardly positioned along each elevation zone, but positions were restricted to flat and gently sloping surfaces lacking deep crevices and tide pools.  Bulleri et al., 2012). In order to simultaneously analyze sessile and mobile species, percentage cover data were transformed to proportions of the maximum observed for each taxon across the shores. In this way, all cover and density data ranged between 0 (the absence of species in a given quadrat) and 1 (maximum value across all quadrats; Caro, Navarrete, & Castilla, 2010;Valdivia et al., 2015). All species abundance estimations were conducted during diurnal low-tide hours.
The spatial or temporal patterns of FRic and FDis can be largely influenced by collinear patterns of species richness. In order to account for this collinearity, we estimated the standardized effect size of FRic (i.e., SESFRic) and FDis (i.e., SESFDis) according to Mason et al. (2008) and Mason, de Bello, et al. (2013). By this method, the observed FRic and FDis values are contrasted to those expected from matrixswap models (Manly & Sanderson, 2002) that randomize species' occurrences (for FRic) and abundances within observations (for FDis).
The SES estimators were then calculated according to Gotelli and McCabe (2002). Null model analyses used 10,000 randomizations. In this way, SESFRic and SESFDis were "pure" estimators of functional richness and the dispersion in trait combination abundances, respectively (Mason et al., 2012). Negative (positive) values of both estimators indicate that functional richness and dispersion are lower (greater) than expected at random.

| Statisticalanalysis
Separate random-slope mixed models were fitted to SESFRic and SESFDis in order to test our prediction. In these models, intertidal height (high, mid, and low intertidal) was included as fixed factor, and sampling site and time (nested in site) were included as random factors. In addition, the autocorrelation of errors, due to repeated measurements over time of each site, was specified as an autoregressive structure (AR1). Model diagnostics were checked by visually inspecting autocorrelation of residuals and fitted-versus-residual plots. We compared the second-order Akaike information criterion corrected for small sample sizes (AICc) in order to select the more parsimonious model that accounted for the variation in SESFRic and SESFDis (Burnham & Anderson, 2004). A delta AICc > 2 was used as decision rule for selecting two competing models (Burnham & Anderson, 2004). The random-slope models allowed us to determine whether the intercepts and slopes of the effects of intertidal height on both estimators of functional diversity varied between sites and between sampling times.
We used similarity routines (SIMPER) to determine the trait combinations that accounted for most of the differences between intertidal heights, sites, and times. SIMPER was based on a Bray-Curtis dissimilarity matrix. Environmental SST and Chl-a data were visualized by means of multivariate redundancy analyses with site as the grouping factor. All analyses were conducted in the R programming environment version 3.3.0 (R Core Team, 2016). We used the libraries FD, mgcv, picante, nlme, MuMIn, vegan, and sciplot.

| RESULTS
Our sampling protocol provided representative estimations of the number of unique trait combinations (i.e., functional richness, FRic) and taxonomic units (i.e., species richness; Figure 2). Asymptotic values of FRic ranged from ca. 10 to 15 unique combinations for highand low-intertidal heights, respectively (Figure 2a). Species richness ranged from ca. 30 in the high intertidal to 55 in the low intertidal ( Figure 2b). According to these collinear patterns between FRic and species richness, the estimation of standardized effect sizes was appropriate.
The analysis of nearshore environmental conditions revealed that the sites located to the north of the study region were characterized by, in average, higher SST and Chl-a concentrations than the southern sites (Fig. S1). This result agrees with previously published analyses  SESFRic showed relatively similar patterns of variation among sites, excepting for the mid-intertidal height, where two sites displayed mean SESFRic values close to zero (Figure 3a). SESFDis showed a larger amount of geographic variation, with a balance between functionally similar and dissimilar species at every environmental stress level (Figure 3b). The results of random-slope models supported these observations, as intertidal height was always selected in the most parsimonious models of the variation in SESFRic and SESFDis (Table 2).
In addition, the fixed-effects coefficients of intertidal height were all positive (Figure 4), supporting the prediction that functional richness and dispersion increase with decreasing environmental stress.
The fixed-effects coefficients from the random-slope models  (Table 2). Therefore, the results of model selection procedure indicated that the relationship between local environmental stress and functional dispersion, but not functional richness, significantly varied over space and time.
Medium-sized, massive, and suspension feeder planktonic developers were the organisms that accounted for most of the Bray-Curtis dissimilarities between intertidal heights, between sites, and between sampling times ( Figure 5). Grazer species were also important to explain the between-group differences for the three factors ( Figure 5).
Small-sized encrusting macroalgae were ranked in third place to account for the variation across the vertical intertidal stress gradient and sampling times (Figure 5a,c). For the differences between sites, F I G U R E 2 Accumulation curves for functional richness (FRic) and species richness for each level of the vertical intertidal environmental stress gradient (high, mid, and low heights). FRic was the number of unique functional trait combinations and species richness the number of taxonomic units F I G U R E 3 Patterns of variation of the standardized effect sizes of functional richness and functional dispersion (SESFRic and SESFDis, respectively) related to the vertical intertidal environmental stress gradient and sampling sites. SESFRic was calculated after randomizing the species occurrences according to a matrix-swap null model. SESFDis was obtained by means of a null model in which abundances were randomized across species, but within sampling plots (Mason, Wiser, et al., 2013). Values are given as means ± confidence intervals

| DISCUSSION
In this study, we showed that different facets of functional diversity had differing patterns of spatiotemporal variation. A significant increase in functional diversity along the vertical intertidal stress gradient was consistent across varying environmental conditions (in terms of SST and Chlo-a concentration) and over more than 200 km and 36 months when the range of trait space covered by the local assemblages (functional richness) was analyzed. When the relative abundances of individuals within trait combinations were incorporated in the analyses (functional dispersion), the slopes between functional diversity and decreasing environmental stress significantly varied between sampling sites and times. Importantly, these effects were independent of the spatiotemporal variation in species richness.
Suspension-feeding planktonic developers and encrusting and foliose autotrophs were the functional traits that accounted for most of the spatiotemporal variation in abundances. Therefore, local environmental stress gradients might have persistent effects on the range of trait space covered by the local communities, but context-dependent effects on the trait combination abundance distribution. Our results stress the relevance of analyzing multiple facets of functional diversity to understand the response of community assembly mechanisms to the environmental variability.
In agreement with our hypothesis, these results provide evidence for trait convergence and divergence along a steep environmental gradient over very small spatial scales. The negative values of the standardized effect sizes for functional richness (SESFRic) indicated that high-intertidal assemblages encompassed functionally similar organisms, which were characterized by small body sizes. In general, trait convergence has been usually explained as the result of environmental filters that allow a subset of organisms with similar trait combinations to colonize and survive in patchy environments (e.g., Leibold et al., 2004). In central Chile, small-sized barnacles that share similar environmental tolerances in the high intertidal show no F I G U R E 4 Coefficients of random-slope models for the effects of the vertical intertidal environmental stress gradient (average change) on SESFRic and SESFDis at four rocky shore sites and over 3 years. SESFRic and SESFDis were calculated according to Mason, Wiser, et al. (2013). For each estimator of functional diversity, the coefficients represent the net change from high-to low-intertidal heights at each site and sampling time. The coefficients were calculated from the full models in Table 2. Wi and Su are austral winter and summer, respectively  For each competing model, fixed (intertidal height = H) and random (site = S and time = T(S)) components, in addition to the degrees of freedom (df) and log-likelihood function (logLik), are given. Model were selected according to second-order Akaike information criterion, corrected for small sampling sizes (AICc), Δi, which is the difference between the AICc of each model and the smaller AICc, and Akaike weights (ω i ).
T A B L E 2 Summary of results of random-slope model selection for the standardized effect sizes of functional richness and functional dispersion (SESFRic and SESFDis, respectively) evidence for competitive advantages of one species over the other, suggesting that ecological equivalence can be common at this height of the shore (Shinen & Navarrete, 2010, 2014)-complementarity in other traits, such as filtration rates (e.g., Valdivia et al., 2009), should be further analyzed to assess niche partitioning in this model system.
At the low-intertidal fringe, we observed that SESFRic was greater than expected at random, indicating the occurrence of species with dissimilar functional traits and an expanded coverage of the trait space (Pavoine & Bonsall, 2011). Predation of dominant competitors such as medium-sized mussels by shallow subtidal seastars and muricid gastropods could facilitate the establishment of subordinate species (Moreno, 2001;Navarrete & Castilla, 2003;Paine, 1966) with different traits and thus increase functional divergence. In other systems, grazing limpets have been shown to exert a firm control on the upper limit of low-intertidal macroalgae, preventing substrate monopolization (Boaventura et al., 2002). Our results further suggest that the role of niche complementarity in determining species occurrences and abundances increases with decreasing local environmental stress (MacArthur & Levins, 1967;Mason, de Bello, et al., 2013). Testing this hypothesis will require field-based manipulations of functional diversity (e.g., Aguilera & Navarrete, 2012;Griffin, Mendez, Johnson, Jenkins, & Foggo, 2009) across intertidal heights, which will foster our understanding of the degree to which niche complementarity influences the assembly of these communities.
Contrarily to functional richness, the patterns of standardized effect sizes of functional dispersion (SESFDis) showed significant variations between sites and between seasonal samplings. This spatiotemporal variation in functional dispersion would be the result of changing environmental conditions strong enough to influence species abundances, but not enough to produce local or temporal extinctions.
The higher Chl-a concentrations observed at the northern sites could imply a higher food supply for planktotrophic larvae and adult suspension feeders, leading to greater invertebrate colonization rates and abundances (e.g., Menge, Gouhier, Hacker, Chan, & Nielsen, 2015). In our study region, the sites exposed to higher Chl-a concentrations also exhibit higher invertebrate recruitment rates (Van Holt et al., 2012;N. Valdivia, unpblish. data), which would explain the spatial variation in the abundance of planktonic developers and suspension feeders observed in this study. The sites with greater Chl-a concentrations also showed warmer waters. Along Chilean shores, and particularly in northern and central regions, the links between nearshore oceanography and coastal ecological processes are strongly related to upwelling activity (e.g., Lagos, Castilla, & Broitman, 2008;Navarrete et al., 2005).
In central Chile, for example, upwelling-related nutrient inputs are associated with greater abundances of corticated algae, but lower abundances of ephemeral algae and grazing effects (Nielsen & Navarrete, 2004). In southern central Chile, the upwelling center locked to Punta Galera (south of Chaihuín, Figure 1) might well encompass a source of environmental variability for coastal communities. Indeed, the sites in our study showing the higher Chl-a concentrations and SST are located in the lee of this upwelling center (Pinochet-Velásquez, 2015).
Small bays downwind from upwelling centers can be areas where surface waters are retained close to the shore, resulting in warmer waters (Piñones, Castilla, Guiñez, & Largier, 2007) and greater larval retention (e.g., Morgan & Fisher, 2010). Nevertheless, the activity of the upwelling at Punta Galera is highly limited to spring and summer (Letelier et al., 2009), and the SST difference between the upwelling and shadow areas is relatively small (Pinochet-Velásquez, 2015). Further studies are needed to determine the effects of this upwelling center on the spatial variability of functional traits in this region. Upwelling effects on species abundances and several relevant ecological processes have been shown worldwide, with stronger effects localized at areas characterized by more intermittent upwelling activity (e.g., Menge & Menge, 2013). In this way, nearshore oceanographic conditions can have strong effects on the abundance of functional trait related with, at least, dispersal potential and trophic strategies.
In addition to abiotic environmental factors, anthropogenic impacts should be considered to explain the patterns of functional trait variation along Chilean shores. Human harvesting in these shores has significantly altered the assemblages at the landscape scale. In particular, the overexploitation of keystone predators, such as muricid gastropods, has severely reduced their abundances on the shore, despite small-sized individuals can be still observed (Moreno, 2001). This could increase the occurrences of small-sized carnivore species (e.g., Acanthina monodon; Moreno, 2001), leading to trait convergence in terms of similarity in body size. A similar pattern has been observed for keyhole limpets, and small-sized herbivores seem to dominate F I G U R E 5 Proportional contribution of unique functional trait combinations to between-group Bray-Curtis dissimilarities for each source of variation (intertidal height, site, and time). Each contribution was calculated as the mean contribution of all of the average pairwise dissimilarities for each factor. A cutoff of 70% of contribution was used. Trait combination codes are as in Table 1- the assemblage of grazers in these shores (Tejada-Martinez et al., 2016). Moreover, high-intertidal macroalgae, such as Pyropia orbicularis, and low-intertidal kelps, such as D. antarctica and L. spicata, are also under severe harvesting (Moreno, 2001). The region under study also shows spatial heterogeneity in these human impacts. Three of the sites are managed, while one is under an open-access regime.
Between the enforced sites, there are also important differences in the consequences of management for the natural communities, depending on differences among stakeholders and local environmental conditions (Van Holt et al., 2012). Our results therefore could have important applications for conservation. For instance, the fact that the relative abundances within functional traits change-despite no trait combination is lost or added at ecological scales-can be seen as evidence for imminent functional extinctions and turnover (Boersma et al., 2016;Säterberg, Sellman, & Ebenman, 2013). Accordingly, the analysis of the combined variation of functional richness and dispersion can be a useful tool to prevent the extinction of critical functional traits.
In summary, our results indicate that considering different facets of functional diversity allows the identification of differing spatiotemporal patterns of variation of community assembly processes.
Although the mechanisms that generate variation in functional diversity were not experimentally assessed in this study, the well-known effects of the vertical environmental stress gradient in rocky intertidal habitats support the idea that the interplay between environmental filters and niche complementarity generates a continuum between functional convergence and divergence. Both environmental and anthropogenic factors could be analyzed to explain the significant spatial (horizontal) and temporal variation in the patterns of functional dispersion along the vertical intertidal environmental gradient. Future work on this issue could help to disentangle the impact of human disturbances on the functional structure of aquatic and terrestrial communities worldwide. For instance, a trait-based approach would allow us to assess anthropogenic impacts across biogeographic regions with different levels of taxonomic and phylogenetic diversity.
The integrated analysis of different facets of functional diversity can provide useful information about the variation and potential for local extinction of functional traits that ultimately influence ecosystem functioning.

ACKNOWLEDGMENTS
We thank all the enthusiastic undergraduate interns (Biología Marina; UACH) that provided invaluable help during the fieldwork. Carlos Lara helped with satellite data downloading and preprocessing. This study was financially supported by the FONDECYT grant #1141037 to NV.