Depth-Dependent Diversity Patterns of Rocky Subtidal Macrobenthic Communities Along a Temperate Fjord in Northern Chilean Patagonia

Villalobos, Vicente I.; Valdivia, Nelson; Försterra, Günter; Ballyram, Stacy; Espinoza, Juan Pablo; Wadham, Jemma L.; Burgos-Andrade, Katherine; Häussermann, Vreni

doi:10.3389/fmars.2021.635855

ORIGINAL RESEARCH article

Front. Mar. Sci., 04 May 2021

Sec. Marine Ecosystem Ecology

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.635855

This article is part of the Research Topic Oceanography and Benthic Ecology of Patagonian fjords - 500 years from the Discovery of the Strait Magellan View all 16 articles

Depth-Dependent Diversity Patterns of Rocky Subtidal Macrobenthic Communities Along a Temperate Fjord in Northern Chilean Patagonia

$\r\nVicente I. Villalobos*$ Vicente I. Villalobos^1*

Nelson Valdivia^1,2

Günter Försterra³

Stacy Ballyram^4,5

Juan Pablo Espinoza^4,5

Jemma L. Wadham⁶

Katherine Burgos-Andrade¹

Vreni Häussermann⁷

¹Instituto de Ciencias Marinas y Limnológicas, Facultad de Ciencias, Universidad Austral de Chile, Valdivia, Chile
²Centro FONDAP de Investigación de Dinámicas de Ecosistemas Marinos de Altas Latitudes (IDEAL), Universidad de Concepción, Concepción, Chile
³Escuela de Ciencias del Mar, Facultad de Recursos Naturales, Pontificia Universidad Católica de Valparaíso, Valparaíso, Chile
⁴Huinay Scientific Field Station, Puerto Montt, Chile
⁵Programa de Magíster en Oceanografía, Escuela de Ciencias del Mar, Facultad de Recursos Naturales, Pontificia Universidad Católica de Valparaíso, Valparaíso, Chile
⁶Department of Geographical Sciences, Bristol Glaciology Centre, University of Bristol, Bristol, United Kingdom
⁷Facultad de Economía y Negocios, Universidad San Sebastián, Puerto Montt, Chile

Understanding the distribution of biodiversity along environmental gradients allows us to predict how communities respond to natural and anthropogenic impacts. In fjord ecosystems, the overlap of strong salinity and temperature gradients provides us with the opportunity to assess the spatial variation of biodiversity along abiotic environmental gradients. However, in Northern Chilean Patagonia (NCP), a unique and at the same time threatened fjord system, the variation of macrobenthic communities along abiotic environmental gradients is still poorly known. Here, we tested whether macrobenthic species diversity and community structure followed systematic patterns of variation according to the spatial variation in salinity and temperature in Comau Fjord, NCP. A spatially extensive nested sampling design was used to quantify the abundance of subtidal macrobenthic species along the fjord axis (fjord sections: head, middle, and mouth) and a depth gradient (0–21 m). The vertical structure of the water column was strongly stratified at the head of the fjord, characterized by a superficial (depth to ca. 5 m) low-salinity and relatively colder layer that shallowed and decayed toward the mouth of the fjord. The biotic variation followed, in part, this abiotic spatial pattern. Species richness peaked at high salinities (>27 psu) between 5 and 10 m in the head section and between 15 and 21 m in the middle and mouth sections. Diversity and evenness were also highest at these salinities and depth ranges in the head and middle sections, but at shallower depth ranges in the mouth. Information theory-based model selection provided a strong empirical support to the depth- and section-dependent salinity, but not temperature, effects on the three biodiversity metrics. Erect algae and the edible mussel Aulacomya atra numerically dominated in shallow water (0–3 m) at the head and the middle of the fjord, coinciding with the horizontal extension of the low-density water layer—these taxa were further replaced by the crustose algae Lithothamnion sp. and deep-dwelling suspension filters (e.g., corals, polychaetes, and sponges) along depth gradient. Macrobenthic biodiversity correlated, therefore, with the influence of freshwater inputs and the density-driven stratification of the water column in this ecosystem. The spatially variable (across both, horizontal and vertical fjord axes) thresholds observed in our study question the widely accepted pattern of increasing biodiversity with increasing distance from the head of estuarine ecosystems. Finally, non-linear environmental stress models provide us a strong predictive power to understand the responses of these unique ecosystems to natural and anthropogenic environmental changes.

Introduction

Biodiversity plays a central role in maintaining the functioning and the goods and services that ecosystems provide to human society (Cardinale et al., 2012). Predicting changes in community measures such as species diversity and composition along environmental gradients is, thus, one of the main objectives of ecology (Hutchinson, 1953; Paine et al., 2018). Therefore, defining the relationships between abiotic environmental factors and the spatial patterns of biodiversity, particularly in ecosystems characterized by pronounced environmental gradients, is relevant for fundamental ecology and helps to inform ecosystem management and conservation strategies (e.g., Dayton, 1971; Underwood et al., 2000; Palacios et al., 2021).

Fjords are ideal natural laboratories to study the effect of abiotic environmental variations on the distribution of benthic communities (Brattegard, 1980; McGovern et al., 2020). The characteristic U-shape of a fjord (Howe et al., 2010) provides deep, large, and almost-vertical rock walls that sustain dense and diverse assemblages of marine invertebrates (Gasbarro et al., 2018). Some studies demonstrate the associations between community proxies and environmental factors along southern fjords (e.g., sediment properties: Gutt et al., 1999; salinity and predation: Smith and Witman, 1999; surface inclination; Cárdenas and Montiel, 2015; terrestrial organic matter inputs: Quiroga et al., 2016). However, most studies that treat quantitative aspects of fjord benthic biodiversity have been conducted in the northern hemisphere, and the study of southern fjord benthic communities (except for several emergent studies on Antarctic systems) has been focused mostly on species inventories and the description of distribution patterns (e.g., Häussermann et al., 2013; Betti et al., 2017).

The coastline of Chilean Patagonia extends 1,500 km from 41° to 56° S, with a total of ca. 100,000 km of coastline; it represents one of the longest fjord regions in the world and is characterized by a high hydrographic and geomorphological complexity (Pantoja et al., 2011). The Patagonian fjord region sustains important supporting ecosystem services, such as clean and cold water for aquaculture, fast currents for hydroelectric energy, cultural services such as tourist destinations, and supporting services such as biodiversity and primary productivity (Iriarte et al., 2010; González et al., 2013; Outeiro et al., 2015). This region is a biodiversity hotspot (Fernandez et al., 2000; Gutt et al., 2003) that harbors unique assemblages of marine invertebrates (Försterra et al., 2017).

Along the Chilean Patagonian fjords, the freshwater input into the fjord produces strong vertical and horizontal gradients in density (Dávila et al., 2002). The hydrography of this system is characterized by a highly stratified water column with two layers of different salinities: a shallow low-salinity layer (ca. 2–20 psu) and the deeper high-salinity layer (ca. 30 psu), which limit at 5–10 m depth—this border tends to shallow from the head of the fjord toward the ocean (Valle-Levinson et al., 2007; Castillo et al., 2016). Besides, the upper 30–50 m are influenced by seasonal variability in solar radiation, freshwater inputs, and mixing induced by wind and tides (Sobarzo Bustamante, 2009; Pérez-Santos et al., 2014). In this way, the diversity and structure (i.e., the combination of species incidences and abundances) of benthic communities would be expected to show non-linear patterns of variation along the different sections of the fjord (i.e., from the head to the mouth), depth, and salinity. This prediction is based on the premise that subtidal fjord communities are composed mostly by marine organisms (Dahl, 1956; Bulger et al., 1993), which can be physiologically sensitive to variations in water salinity (Smyth et al., 2014). According to early and recent environmental stress models, the interaction between physiological constrains imposed by abiotic factors, dispersal limitations, and density-dependent biotic interactions have been shown to lead to non-linear diversity-environment relationships in different marine ecosystems (e.g., Menge and Sutherland, 1987; Scrosati and Heaven, 2007; Thompson et al., 2020). In rocky intertidal habitats, for example, the overlap between species adapted to the physiologically demanding conditions of the supralittoral and those adapted to the infralitoral generates unimodal patterns of species richness, diversity, and evenness, in addition to a significant variation in community structure along the vertical intertidal stress gradient (e.g., Menge and Sutherland, 1987; Scrosati and Heaven, 2007; Heaven and Scrosati, 2008). In subtidal fjord habitats, however, local adaptation could lead to environmentally decoupled community patterns (Basset et al., 2013; Elliott and Quintino, 2019). For Chilean Patagonian fjords, therefore, the association between vertical and horizontal abiotic gradients and ecological communities still needs to be determined.

Here, we use a spatially extensive sampling program to examine the relationships between abiotic environmental gradients and the distribution of macrobenthic species in a Northern Chilean Patagonian fjord (NCP). Three hypotheses were tested:

H1) Since the influence of freshwater inputs decreases from the head to the mouth of the fjord and the fjord benthos is colonized mostly by marine species (Bulger et al., 1993), we predict an increase in macrobenthic species richness, diversity, and evenness, and a major change in community structure, from the head to the mouth of the fjord and with increasing salinity.

H2) The sharp vertical (along depth) variation in salinity, which shows a threshold at 5–10 m depth, should have non-linear effects on the diversity and structure of fjord macrobenthic communities: species richness, diversity, and evenness should peak between 5 and 10 m depth because of the overlap between species occurring at shallow (above 5–10 m) and deep (below 5–10 m) substrata, leading to a unimodal environmental-diversity relationship. In addition, and assuming the existence of strong abiotic environmental filtering in these ecosystems (Smyth et al., 2014), then we predict a notable change in community structure between the shallow and deep strata.

H3) Since the vertical structure of salinity varies horizontally from the head to the mouth of the fjord, we expect that richness, diversity, and evenness will peak—and community structure will change—at different depths along the fjord sections. This pattern would be represented by a strong empirical evidence supporting an interactive effect of fjord section, depth, and salinity on univariate and multivariate community metrics.

To test these hypotheses, we analyzed the horizontal (along the fjord) and vertical (along increasing depth) variation of species abundances and salinity of the water column in Comau Fjord, NCP. In addition, we analyzed the spatial variation in water temperature as a potentially relevant environmental covariable in these ecosystems. This fjord harbors a unique subtidal ecosystem, inhabited by speciose communities of deep-water organisms, like corals, that thrive below 20 m depth (e.g., Häussermann and Försterra, 2009; Häussermann et al., 2013). However, rampant anthropogenic direct and indirect pressures, such as the multiple impacts derived from aquaculture activities, industrial and artisanal fisheries, and climate change, pose severe threats to this exceptional and complex ecological system (Försterra et al., 2017; Iriarte, 2018; Soto et al., 2019).

Materials and Methods

Study Site

The study was conducted in the Comau fjord, located in NCP (ca. 42°S; Figure 1). The fjord is 34.3 km long, its widest part is at the mouth (ca. 10 km), and the narrowest (ca. 2 km) is near the head. The geomorphology of the fjord creates an almost north-south axis (along ca. 72.47°W) of environmental variation, oriented toward 346° in relation to the geographical north (Häussermann and Försterra, 2009). The main freshwater inputs are located on the eastern side (Quintupeu and Cahuelmó fjords) and at the head (Vodudahue and Leptepu rivers) of the fjord. The Comau fjord has an over-deepened sill, which could lead to shorter bottom water residence times relative to southernmost fjords (Häussermann and Försterra, 2009). The mouth connects the fjord with the Ancud Gulf, across the Marilmó and Comau passes, and to the Hornopirén channel, across the Cholgo channel. The fjord depths are around 500 m along the entire fjord, to <50 m very close to the head, due to the estuary proximity (Häussermann and Försterra, 2009).

FIGURE 1

Figure 1. Study site and experimental design. (A) Map of the Comau Fjord with the sampling sites (circles) nested within the head, middle and mouth sections of the fjord (colors). (B) Northern Patagonia with Comau fjord highlighted by a red square. (C) Section of the Southern Cone of South America with the Northern Patagonia highlighted by a red square.

The Comau fjord follows an estuarine circulation with a marked vertical stratified structure of the column water mainly determined according to the salinity and temperature regimes. This accentuated structure fluctuates annually in position and intensity mainly influenced by the annual cycles of solar radiation, precipitation, and tides mixing (Försterra, 2009; Sobarzo Bustamante, 2009). The surface layer is very variable: thicker, colder, and lower salinity in winter; warmer, and relatively higher salinity in summer. The main water body below this superficial layer has a relatively constant salinity and temperature throughout the year. On average, the surface low-salinity layer varies between 2 and 20 psu and has temperatures that range between 7 and 20°C; the deep sea-influenced water layer has a salinity > 32.5 psu and a temperature of ca. 11°C. Both layers converge vertically at ca. 10 m (Häussermann and Försterra, 2009). However, anomalous atmospheric events, such as the transition to a positive phase of the Southern Annular Mode and the intensification of the positive phase of the El Niño Southern Oscillation, have been shown to account for the recent decay in freshwater inputs in northern Patagonia (Boisier et al., 2018; Aguayo et al., 2019) and the concomitant weakening of the haline stratification of NCP fjords (León-Muñoz et al., 2018).

Environmental Gradients

The physicochemical characteristics of the Comau Fjord were measured at 64 stations across the fjord during February 19–21, 2020. At each station and between 0 and 100 m depth, temperature (° C) and salinity (psu) were recorded at intervals of 0.5 m with a CTD AML Plus X oceanographic probe. Then, an Inverse Distance Weighting and Nearest Neighbor interpolation was made to extract values at 0, 1, 3, 6, 10, 15, and 21 m at each ecological sample sites. The vertical profiles were carried out. The Brunt-Väisälä frequency or “buoyancy frequency” was calculated from the temperature and salinity values. This parameter was used to determine the vertical stratification of the water column, according to Stewart (2008):

N^{2} = - g E,

where N is the buoyant frequency, g is gravity and E is the stability of the water column, which was calculated according to the expression:

E = - \frac{1 δ ρ}{ρ δ z},

where ρ is the water density (1.025 kg m^–3) and z is the depth. Positive values of N² represents water column stratification (i.e., a low-density layer located above a high-density layer). Therefore, when the values of N² are high, the water column is highly stratified and vice versa. On the other hand, negative values of N² represents water column mixing (Garcés-Vargas et al., 2013). Surface and vertical distributions of environmental variables were plotted using Ocean Data View 5.2.0 software (ODV; Schlitzer, 2020).

Ecological Sampling Design

To test our predictions, we used a nested sampling design composed of three 10-km² sections (head, middle, and mouth) located along the fjord and separated by ca. 10–30 km (Figure 1). In each section, four coastal sites, separated between 1 and 5 km, were randomly chosen, except for the head section where only three sites were located due to logistical limitations. At each site, a 26-m alongshore transect consisting of ten rectangular 0.28-m² (62 × 45 cm) observational units (OUs), was positioned at seven depth strata (0, 1, 3, 6, 10, 15, and 21 m). The experimental units were equidistantly distributed within each transect. This design yielded a total of 741 OUs. To survey exclusively within the subtidal range, the 0 m depth was set at the lowest tide mark relative to the historical zero. All transects consisted of almost vertical hard rocky substrates with some small patches of soft sediments. The design was orthogonal, considering the effect of depth and sections (both crossed), and at the same time hierarchical, by considering the effect of sites (nested in sections) as explanatory variables. The study was conducted during two field campaigns in November–December 2017 and July 2018.

Species Abundance Estimations and Diversity Measures

The abundances of sessile and mobile macrobenthic species (>0.5 mm) were measured in each OU by means of underwater photographs (photo-quadrats) taken by the team of SCUBA divers of the Huinay Scientific Field Station. A frame with graduated strips of 62 × 45 cm was placed over each OU. This plot size was chosen to maximize the image quality associated with the camera’s optimal size resolution. The area covered by the frame was photographed with a still Olympus E-m5 OMD camera, equipped with an Olympus M. Zuiko 9–18 mm F4-5.6 lens and a Nautikam NA-EM5 housing. The camera was mounted on a tripod with two Inon D2000x2 strobe lights.

To quantify the percentage cover of each macrobenthic species on each photo-quadrat, we focused on the primary-space holders (i.e., basal species). However, when superposition of organisms occurred, we included secondary-space holders (sessile species attached to the primary-space holders) and the associated mobile species. The images were processed with the Coral Point Count with Excel extensions (CPCe) v4.1 freeware (Kohler and Gill, 2006) at the Laboratory of Coastal Ecology of the Universidad Austral de Chile (UACh). For each image, the quadrant area was first delimited with a digital border, and then a squared grid of 100 equidistant points was established. Each readily visible element found at each point was identified either as species (invertebrates or algae) or characteristic of the benthos (e.g., substrate type). The number of points associated with each species was used as a measure of relative abundance of species (i.e., percentage cover). Points that intersected elements such as shadows or the frame were excluded from the calculations.

All organisms were identified at the finest possible taxonomic resolution with aid of specialized taxonomic guides (Häussermann and Försterra, 2009; WoRMS, 2020). The main criterion for determining species identity was to achieve at least one morphological feature characteristic of the species.

Species richness (S) was expressed as the total number of taxa in each OU. The Shannon diversity (H’; Shannon, 1948) was expressed as $- \sum_{i = 1}^{S} (p_{i} \times \log_{2} p_{i})$ , where p_i is the proportional abundance of species i and combines both species richness and relative abundances. The Pielou evenness (J’; Pielou, 1969) was expressed as $\frac{H^{^{'}}}{\log_{e} S}$ , where log_eS is the maximum diversity ( $H^{'}_{m a x}$ ). J’ summarize the dominance structure of the community and ranges from 0 to 1 (Magurran, 1988). The Bray-Curtis coefficient (BC; Bray and Curtis, 1957) was expressed as $\frac{\sum_{i = 1}^{S} | y_{i j} - y_{i k} |}{\sum_{i = 1}^{S} (y_{i j} y_{i k})}$ , where y_ij and y_ik are the abundances of species i for the observation units j and k, respectively. BC combines species occurrences and relative abundances and is used as a measure of dissimilarity (or distance) between pairs of OUs. BC ranges from 0 (same species and abundances in both objects, i.e., maximum similarity) to 1 (no species with abundance > 0 in common, i.e., maximum dissimilarity).

Statistical Analysis

Data were analyzed with generalized additive mixed-effect models (GAMM or also Hierarchical GAM; Wood, 2017; Pedersen et al., 2019). These models incorporate smooth functions of one or more covariates and can analyze non-linear relationships between covariates and the response variable including random effects. Each diversity measure (S, H’, and J’) was analyzed in individual models. The effect of temperature (°C), salinity (psu), depth (m), and fjord section were inferred in an information theoretic approach for model selection (Burnham and Anderson, 2002). The focus of these analysis was to assess the inter-dependent (i.e., interactive) effect of temperature, salinity, and depth combinations on each diversity measure across the fjord sections. Thus, we followed a parsimonious approach and constructed a family of GAMMs that included the main effect of fjord section, the interactive effect of temperature, salinity, and depth combinations across fjord sections, and the random effect of fjord sites for each diversity measure (Table 1).

TABLE 1

Table 1. Description of the models constructed for the information-criteria-based model selection.

To avoid problems of overfitting and increased residual variation derived from collinearity between continuous variables (Dormann et al., 2013), the Pearson correlation coefficient was calculated for depth, salinity, and temperature combinations (Supplementary Figure 1). Temperature-depth correlation was high (r = 0.84). Therefore, the models that included temperature and depth combinations were excluded from the analysis, since depth often corresponds with driving temperature gradients and can be used as their proxy in statistical models (McArthur et al., 2010).

Graphical inspections of the residuals suggested that the Poisson (log link) and Gaussian (identity link) structure of errors were appropriate for species richness (S) and for Shannon diversity (H’) and Pielou evenness (J’), respectively. The smooth functions for depth, temperature, and salinity were fitted using thin plate regression splines and the maximum number of basis functions was set to seven. Model parameters were estimated through maximum likelihood. Bias-corrected Akaike Information Criteria (AICc) was used to select the top model accounting for the data (Akaike, 1973). Model selection was based on Δ_i, expressed a AICc_i–AICc_min, Akaike weight, which represents the probability of each model as $w_{i} = P r o b {m o d e l g_{i} | d a t a} = \frac{l_{i}}{\sum_{j = 1}^{R} l_{j}}$ , where l_i = L(g_i|data) = e^−0.5_i, and the evidence ratio of the top model, calculated as w_top / w_j (Burnham et al., 2011).

Permutational multivariate analysis of variance was used to analyze community structure (PERMANOVA; Anderson, 2001; Anderson and Braak, 2003). We used the best model selected in univariate analysis. Therefore, the model included the effect of depth ranges, salinity, fjord section, their interactions, and fjord sites (as a permutation constrain) on the variations in community structure (see section “Model Selection” below). The Bray-Curtis coefficient was used as a measure of dissimilarity with 999 random permutations. The patterns of variation in community structure were examined using a non-metric multidimensional scaling chart (NMDS) based on Bray-Curtis dissimilarities. In addition, the percentage contribution of each taxon to the between-group variation in community structure was calculated with a (di)similarity percentage analysis (SIMPER) based on Bray-Curtis dissimilarities. For all multivariate analyses, the depths were grouped into four depth ranges (0–1 3–6, 10–15, and 21 m).

All statistical analyses and plots were performed in the R software, version 4.0.3 (R Core Team, 2021). The “raster” package was used for the Inverse Distance Weighting and Nearest Neighbor interpolation (Hijmans, 2020); “vegan” for the calculation of S and H’, MDS, PERMANOVA, and SIMPER analyses (Oksanen et al., 2019); “car” for the correlation scatterplot (Fox and Weisberg, 2018); “mgcViz” for the GAMM visualization (Fasiolo et al., 2018); and “mgcv” for the GAMMs (Wood, 2011).

Results

Summer Environmental Conditions of Comau Fjord

Surface salinity varied from 5 psu at the head to 25 psu at the mouth of the fjord (Figure 2A). The vertical distribution of salinity showed a superficial freshwater-influenced layer (ca. 20 psu; < 5 m) placed above a homogeneous ca. 32-psu layer along the fjord (Figure 2B). At the head, lowest salinity values were registered (ca. 5–17 psu) in the first 3 m, followed by a strong halocline at 4–5 m (ca. 19–30 psu), which gradually shallowed along the middle section until almost disappearing at the mouth (Figure 2B). The surface temperature varied from 13°C at the head to ca. 16 at the mouth of the fjord (Figure 2C). The vertical distribution of temperature showed a warmer layer (ca. 14–17°C; < ca. 15 m) placed above a colder layer along the central part of the fjord (ca. 11–13°C; Figure 2D). However, the head section registered a superficial cold-water influence from the riverine input (ca. 13°C; < 3 m) and the head section showed a homothermic- ocean-influenced layer (14–13°C) along the column water (Figure 2D).

FIGURE 2

Figure 2. Environmental variables distribution along Comau Fjord. CTD castings data showing (A) superficial and (B) longitudinal profile of temperature, (C) superficial and (D) longitudinal profile of salinity, and (E) longitudinal profile of Brunt-Väisälä frequency. All the physicochemical characteristics were measured during February 2020.

The Brunt-Väisälä frequency distribution exhibited a highly stratified water column (75–125 cycl/h) at the first ca. 3 m along the fjord (Figure 2E). However, a considerable but minor difference in the vertical water density was still observed at 5 and 10 m (50 and 25 cycl/h, respectively) along the head and middle sections (Figure 2E).

Cumulative Diversity Along the Comau Fjord

A total of 130 taxa were identified in the Comau fjord (among which 49.2% identified at species, 31.5% at genus, 0.1% at family, and 18.4% at class level). One-hundred and fourteen taxa of sessile and mobile macroinvertebrates and 18 taxa of macroalgae were identified (Supplementary Table 1). The cumulative richness was higher at the middle and the mouth of the fjord than in the head section (Figure 3A), which partially support our prediction about the horizontal gradient. In the case of depth strata, the higher values of species richness were registered at 10 and 15 m and the minima at 0 and 1 m (Figure 3B).

FIGURE 3

Figure 3. Species-accumulation curves (SACs) of macrobenthic communities of Comau fjord. Separates SACs were performed for each (A) fjord sections and (B) depth strata. Cumulative number of species is shown as mean from each random combination of OUs (0.28 m²).

The sites that showed the lowest species richness and diversity (i.e., sites 9 and 11 at 21 m at the head section; Supplementary Figure 2) were visually related to a substrate with “poor” quality sediment conditions, such as high homogeneity and grayish or blackish sediments (Supplementary Figure 3). Both sites also coincided with the location of a salmon (cage raft type) and mussel (longline type) farm, respectively. On the other hand, the third and last site of the head section, which did not show that a dramatic decrease in richness and diversity (Supplementary Figure 2), was related to a more hard, vertical, and clean substrate.

Model Selection

For species richness (S’), Shannon diversity (H’), and Pielou evenness (J’), the model with the lowest AICc was that included the separate and interactive effects of depth and salinity for each fjord section, and the random effect of sites nested in fjord sections (g₆ in Table 1 and Supplementary Table 2). The probability of g₆ of being the best model given the data (Akaike weight, w₆) was ca. 0.9999 for the three diversity measures. Thus, the empirical support for model g₆ was more than one million times that of the closest competing model at each set of models (S’: evidence ratio, ER₄ = 1422582; H’: ER₅ = 1.5 ^∗ 10⁸; and J’: ER₇ = 2.0 ^∗ 10⁸). Importantly, g₆ was 1.3 ^∗ 10¹³⁰, 6.6 ^∗ 10⁶⁴, and 3.2 ^∗ 10³⁹ for S, H’, and J’, respectively, more likely than the null model (g₈ in Tables 2–4, respectively); this provides a very strong empirical support to the effect of the depth, salinity, and fjord section on S’, H’, and J’. Model g₆ accounted for a 46, 38, and 26% of the variation in S’, H’, and J’, respectively (r²₈ in Tables 2–4, respectively).

TABLE 2

Table 2. Summary of information-criteria-based model selection for the species richness (S).

TABLE 3

Table 3. Summary of information-criteria-based model selection for the Shannon diversity (H’).

TABLE 4

Table 4. Summary of information-criteria-based model selection for the Pielou evenness (J’).

Species richness (S) increased with increasing depth and peaked at 5–10 m (27–30 psu) in the head section and at 15–21 m (27–30 psu) in the middle and mouth sections (Figures 4A–C). Yet, the environmental effects on S’ were strongest in the head (compare effect magnitudes in Figures 4A–C). Diversity (H’) and evenness (J’) followed a similar pattern for the head and middle fjord sections (Figures 4D,E,G,H). On the other hand, at the mouth section, H’ and J’ presented an increase at the first three meters and highest salinities, and reduced values along higher depths (Figures 4F,I). These results support the hypothesis of interdependent effects of spatial and environmental factors on biodiversity (H3).

FIGURE 4

Figure 4. Estimated smooth interactive effect of depth and salinity on diversity measures (species richness, Shannon diversity, and Pielou evenness’s) as a function of the generalized additive mixed model (GAMM) output. (A–C) Estimated interactive effect of depth and salinity on species richness. (D–F) Estimated interactive effect of depth and salinity on Shannon diversity. (G–I) Estimated interactive effect of depth and salinity on Pielou evenness. The gray zones represent absent depth-by-salinity combinations.

Patterns of Community Structure

The analysis of permutational multivariate variance (PERMANOVA) supported the hypothesis of significant interactive effects of depth, fjord sections, and salinity (p < 0.001) on community structure (H3; Supplementary Table 3). Non-metric multidimensional scaling ordination (NMDS) for the community structure showed that at the head section of the fjord a vertical pattern with a “shallow” group (0–1 m) differed from a “deeper” group (21 m), and an intermediate group (3–15 m) covering the almost complete gradient (Figure 5A). On the other hand, at the middle and mouth sections, a more homogeneous vertical pattern was observed between 3 and 21 m with only the 0–1 stratum differing from the rest of the depth ranges (Figures 5B,C). The same method of ordination showed that within the 0–1 stratum each fjord section presented a group of quadrats that differ between them and a shared group (Figure 6).

FIGURE 5

Figure 5. Non-metric Multidimensional Scaling (NMDS) ordination plots of the observational units based on Bray-Curtis dissimilarities between depth ranges at (A) head, (B) middle, and (C) mouth sections of Comau fjord.

FIGURE 6

Figure 6. Non-metric Multidimensional Scaling (NMDS) ordination plot of the observational units based on Bray-Curtis dissimilarities between fjord sections at 0–1 column water stratum of Comau fjord.

In general, the most dominant taxa in Comau Fjord were the crustose algae Lithothamnion sp. and the mussel Aulacomya atra, contributing ca. 50% to the dissimilarity between sections and depths (Table 5A). The relative abundances of both species appeared to have inverse relationships with depth and fjord sections, where A. atra dominated the first 3 m with the highest abundances at the head section, while Lithothamnion sp. dominated below 6 m with the highest abundances at the mouth section (Figure 7). Subdominant taxa like the limpet Crepipatella sp., the brown macroalgae Macrocystis pyrifera, the filamentous macroalgae Ectocarpus sp. and Acrosiphonia sp., and the foliose macroalga Ulva sp. dominated at 0–1 m. The encrusting anemone Epizoanthus fiordicus, the soft coral Clavulariidae spp., and the polychaetes Spiochaetopterus sp. and Chaetopterus variopedatus dominated below the 10–15 strata, following an intermediate zone (3–6 m) mostly dominated by Lithothamnion sp. (Table 5B). Within the head, Lithothamnion sp. was the main contributor to the dissimilarities between intermediate depths (3–15 m) but with significant decreases in abundance at 0–1 and 21 m depths (Table 6). However, A. atra and a group of green (Ulva sp. and Acrosiphonia sp.) and brown (Ectocarpus sp. Ceramiales sp.) macroalgae dominated above the 3 m stratum, while the anthzoans Calvulariidae sp. and E. fiordicus and the holoturid P. disciformis was characteristic for the benthic community below the 10–15 m stratum (Table 6). See Supplementary Figure 4 for a graphical representation of the subdominant taxa along depth and sections of Comau fjord.

TABLE 5

Table 5. Summary of SIMPER results for sections and depth ranges of Comau fjord.

The average abundances (% cover) of relevant taxa in each section group, their contribution (% ± standard deviation) to the within-group dissimilitude, and cumulative total (%) of contributions (80% cut-off) are shown.

TABLE 6

Table 6. Summary of SIMPER results for depth ranges at the head of Comau fjord.

FIGURE 7

Figure 7. Mean relative abundances (percentage cover) of the most abundant species on Comau fjord along fjord sections and depth: (upper panel) the edible mussel Aulacomya atra and (lower panel) the crustose algae Lithothamnion sp. Vertical lines represent the standard errors (± SE).

Discussion

This study revealed that diversity and structure of the macrobenthic subtidal rocky bottom communities of the Comau Fjord are significantly associated with both horizontal and vertical abiotic environmental gradients, particularly salinity, emphasizing their interdependent relationship. A double-layer structure of the water column (0–21 m) was observed along the fjord: A colder freshwater-influenced layer located above a homogeneous salty layer produced a strong density stratified parcel of water at ca. 5 m in the head that shallowed toward the mouth. Our analyses showed a depth-dependent increase of species richness toward high salinities (>27 psu) along the fjord axis. This corroborates the positive relationship of salinity and the number of species along estuarine environments (Whitfield et al., 2012). Peaks of richness were observed at different depths for the head (5–10 m) and middle and mouth (15–21 m) sections. Diversity and evenness followed the trends of richness along depth and salinities at the head and middle sections, but both metrics peaked at shallow depths in the mouth. Community structure changed from macroalgae in shallow depths (0–1 m) to crustose algae in intermediate depths (3–6), to deep-dwelling suspension feeders below 10–15 m at the head and the middle sections. In contrast, the mouth section maintained a community dominated by mobile predators and crustose algae down to 15 m. At 21 m, the mouth was numerically dominated by dense patches of suspension feeders, mainly Anthozoa. These complex non-linear patterns of variation could be accounted for by the overlap of species with different environmental tolerances (i.e., freshwater and oceanic) and the decay of the low-salinity layer toward the mouth of the fjord. In a broader context, our results could be related to the first half of a unimodal pattern of an environmental stress gradient (e.g., Scrosati and Heaven, 2007; Zwerschke et al., 2013), where, in this case, low salinities (at shallow depths and close to the head) represent the upper stress boundary, and high salinities (at deep depths and close to the mouth) represent the middle to lower stress boundary.

Spatial Patterns of Species Richness and Diversity

The vertical patterns of the three diversity measures were dependent of the fjord sections, with maxima observed at different depths along the fjord sections. Generally, this suggests a non-linear, depth-dependent pattern that can be compared with the well-described pattern of decreasing diversity from the head to the mouth of estuarine ecosystems (e.g., Ysebaert et al., 2003; Meire et al., 2005; Beuchel et al., 2006). Our results contradict therefore previous work in other temperate fjords: for instance, Kuklinski (2013) attributes the absence of horizontal trends in temperate fjord benthic communities to a scarcity of environmental gradients in this axis. In this sense, horizontal patterns of fjord epibenthic organisms can be strongly context- and taxonomic group-dependent (Brattegard, 1966; Rosenberg and Möller, 1979; Hansen and Ingólfsson, 1993). In a more general context, our results agree with a recent model of biogeography, the sub-habitat dependence hypothesis (Scrosati et al., 2020), which predicts that a geographical biodiversity pattern within a given region will differ between different sub-habitats (depths and fjord sections in our case). This model assumes that the main abiotic drivers of sessile species distributions varies across the sub-habitat types.

A low-salinity water layer was clearly associated with the zones of low benthic diversity, mainly at the upper ca. 10 m of the column water. In estuaries, a salinity between 24 and 30 psu has been considered as the physiological limit that distinguishes between marine and “brackish” species (Dahl, 1956; Bulger et al., 1993). An average salinity of ca. 10 psu has been suggested to present the minimum diversity of species (Attrill and Rundle, 2002) and the range between 30 and 40 psu the maximum diversity of species (mostly marine species; Whitfield et al., 2012). The tolerance of marine organisms to variations in salinity can be also temperature-dependent, due to the negative effect of the latter on osmoregulatory capacities (e.g., Navarro et al., 2019; Vargas-Chacoff et al., 2019). In our study, however, the empirical evidence supporting the effects of salinity on the diversity measures was much stronger than that of temperature and the interaction of both variables. Nevertheless, the role of temperature as part of the environmental filters and drivers of population dynamics could well be stronger at different moments of the year. For example, salinity tolerance may vary seasonally according to the annual temperature regime (Smyth and Elliott, 2016). Unfortunately, we were unable to analyze the temporal variation in environmental conditions. In addition, the high temperature-depth correlation indicates that part of the vertical patterns in the diversity proxies would well be accounted for by variations in temperature. Further long-term observational studies—coupled with well-designed manipulative experiments—are needed to fully disentangle the effects of salinity and temperature on these communities.

At the mouth of the fjord, the spatial patterns of species richness differed from those of Shannon diversity and Pielou evenness. This could be attributed to differential environmental effects on species incidences (related to abiotic environmental filtering) and species abundances (related to population dynamics). The outer sections of a NCP fjord are characterized by a relatively homogeneous water column with more oceanic, high-salinity environmental conditions (Castillo et al., 2016). Therefore, other factors, besides salinity, are likely driving the assembly and population dynamics at the fjord’s mouth. In this section, for example, milder abiotic environmental conditions could be associated with a stronger influence of density-dependent biotic interactions, as predicted by early and recent environmental stress models tested elsewhere (Menge and Sutherland, 1987; Scrosati and Heaven, 2007).

Community Structure and Environmental Variability in Shallow Waters

At the head of the fjord, the edible mussel Aulacomya atra peaked in abundance at shallow waters (0–3 m depth; Figure 7). This peak coincided with the absence of carnivorous predators. On the other hand, the mouth was characterized by minimum abundances of A. atra and maximum abundances of predators such as Cosmasterias lurida. Low-salinity waters in the head of the fjord may act as a refuge area for prey like mussels (Witman and Grange, 1998; Wing and Leichter, 2011; Sjøtun et al., 2015). For instance, Wing and Leichter (2011) show that, despite the negative consequences caused by reduced salinities on mussel growth (e.g., Navarro, 1988), these organisms can develop dense patches of individuals within the low-salinity layer at the head of a temperate fjord, where salinity is low enough to restrict the access to predators (e.g., Cosmasterias sp.). Furthermore, in a New Zealand temperate fjord, Witman and Grange (1998) determined that predation success in low-salinity stressful habitats strongly depends on the time before predation is interrupted by vertical fluctuations in the low-salinity layer caused by changes in freshwater discharges.

Similarly, the low or null abundance of foliose macroalgae found at shallow depths (0–1 m) at the mouth of the fjord could well be associated with a negative consumptive effect of the sea urchins Loxechinus albus and Arbacia dufresnii in this section. The almost total decrease of these herbivores toward the middle and head sections of the fjord coincided with a significant increase in the abundance of the large, perennial brown algae Macrocystis pyrifera together with the appearance of small, ephemeral algae (i.e., Ectocarpus sp., Acrosiphonia sp., and Ulva sp.). The latter increased significantly in abundance toward the head, co-dominating the substrate next to A. atra. The exclusion of grazers from the low-salinity fjord sections can allow the development of seaweed-dominated communities (e.g., Ayling, 1981; Husa et al., 2014). As with the invertebrate case, the low-salinity layer observed in this study could also be providing an important refuge for macroalgae in the first few meters of the water column.

At the middle and the mouth of the fjord, the red crustose algae Lithothamnion sp. dominated the community. Species in this genus are very well adapted to low light and temperature regimes (Johansen, 2018), which could result in competitive advantage over other algae and in the ability to thrive at greater depths. The penetration of light into the water column in NCP is strongly limited by the presence of particulate material of terrestrial origin in the head zone (Huovinen et al., 2016), which will define the lower limit of distribution of foliose algae in sites where discharges are high (e.g., Barrett and Edgar, 2010). Indeed, the variation of substrate inclination, which influences the degree of sedimentation and light regimes, has been suggested to influence the co-occurrence patterns of macroalgae elsewhere (e.g., South Patagonia, Cárdenas and Montiel, 2015). In addition, salinity can determine the proportion of red and brown algae in the first few meters of the estuaries (e.g., Munda, 1978; Schubert et al., 2011). Red algae appear to be more sensitive to variations in salinity than brown algae (Cole and Sheath, 1990) and it is suggested that salinity of 20 psu or less is the critical lower limit of tolerance for several red algae (Johansen, 2018). This could explain the reduced abundance of Lithothamnion sp. in the first meters of the head and the subsequent increase (up to 90%) of its abundance toward greater depths (see Figure 7).

Habitat-Forming Species and the Marine Animal Forests

The macrobenthic communities observed below 6 m stood out with a high level of clustering, or “patchiness”, composed of groups of multiple invertebrate species associated with an increase in local diversity. Generally, these patches were found to be surrounded by low diversity areas dominated by the crustose algae Lithothamnion sp. These observations coincide with the shallow rocky bottom communities described by Grange et al. (1981) in the New Zealand temperate fjords, Betti et al. (2017) in a fjord in Southern Patagonia, and Cárdenas and Montiel (2015) in the Strait of Magellan, Chile. Sessile suspension feeders like bivalves, sponges, and corals play a role as ecosystem engineers (sensu Jones et al., 1994) and can maintain complex three-dimensional structures that provide habitats and refuges for other species—these associations have been named as “Marine Animal Forests”, because of their structural and functional similarity to terrestrial forests (Rossi et al., 2017). These communities occur widely around the world, from the tropics (e.g., Soares et al., 2017) to the poles (e.g., Arntz and Gallardo, 1994; Cárdenas and Montiel, 2017) and are characteristic of the rocky walls of the Comau Fjord (Försterra et al., 2017).

This study highlights the occurrence of several habitat-forming species such as mussels (e.g., A. atra), brachiopods (e.g., Magellania venosa), polychaetes (e.g., Chaetopterus variopedatus), sponges (e.g., Scopalina sp.), and cold-water coral forests (e.g., Desmophylum dianthus; Försterra et al., 2017). An increase in structural complexity induced by these habitat-forming organisms may be contributing to the maintenance of the observed diversity thresholds. This idea is supported by Smith and Witman (1999) who observed elsewhere that highly diverse benthic patches are maintained by increased habitat biogenic complexity that enhances larval recruitment. Försterra (2009) describes the extraordinary occurrence of “deep-sea species” in the shallow depths of the fjords of Patagonia (e.g., Comau Fjord), a phenomenon that has also been reported in similar fjords in New Zealand (e.g., Grange et al., 1981) and Australia (Barrett and Edgar, 2010). Overall, the facilitative effects of habitat-forming species could be contributing to the maintenance of the diversity found 10–15 m depth in the Comau Fjord (Reise, 2002; Gili et al., 2006).

Natural and Anthropogenic Sedimentation and Their Impacts on Suspension Feeders’ Communities

The richness and abundance of species below 15 meters decreased notably at the head of the Comau fjord, which coincided with a high quantity of low-quality sediments accumulated on the bottom (Supplementary Figure 3). Here, two of the three sites were located adjacent to the outlet of a tributary rivers. The two sites with the highest levels of accumulated sediments also showed the lowest diversity and abundance of benthic species (Supplementary Figure 2). In pristine non-glacial fjords, freshwater discharges in the form of rivers and inland runoff are usually the most important factors in affecting total sedimentation rates and the proportions of organic (<10%) and inorganic (>90%) particles in the sediments (Pearson, 1980; Syvitski et al., 2012). Epibenthic sessile suspension filter communities are highly susceptible to sediment discharges (Evans et al., 1980; Farrow et al., 1983), because of the obstruction of the filtering organs (e.g., Bell et al., 2015) and enhanced substrate instability (Gili et al., 2006). This could explain the decrease in richness and abundance of species found at the head, in addition to the large variation in the diversity proxies accounted for by the differences among sites. A similar pattern has been described for glacial fjords, where the sedimentation from glacial discharges usually leads to a dominance by buried worms and mollusks (e.g., Wlodarska-Kowalczuk et al., 2005; Voronkov et al., 2013; Valdivia et al., 2020).

An alternative and non-exclusive explanation to the punctual low diversity found below 15 m in the head section would stem from negative effects of salmon and mussel farms operating in the fjord. Our observation of low-quality sediments in the nearby of marine farms (see Supplementary Figure 3) suggests the occurrence of a degraded habitat due to an excess accumulation of organic matter and a possible eutrophication event. Benthic communities react according to the deficiencies in oxygen concentration caused by organic enrichment (Hargrave, 2010). In the fjords of Northern Patagonia (e.g., Comau fjord), the input of organic material of terrestrial origin through the multiple freshwater discharges (Villagrán, 1988; Bustamante, 2009) and the efficient export of algal biomass to sediments in productive times (e.g., González et al., 2010) enhance organic matter flows to the benthic system in the head, increasing the possibility of a eutrophication event. Generally speaking, the frequency of hypoxic conditions in estuarine environments has increased significantly, strongly accelerated by human activity (Diaz and Rosenberg, 1995; Breitburg et al., 2018) and negatively affecting the benthic fauna of fjords (Soto and Norambuena, 2004; Hargrave, 2010). The rapid expansion of the aquaculture industry in Chilean Patagonia represents an environmental risk with important associated regional impacts. For example, the overload of ecosystem-level carrying capacity is currently degrading the endemic fjord communities in Chilean Patagonia (Buschmann et al., 2006; Niklitschek et al., 2013; Soto et al., 2019).

Conclusion

The diversity and structure of the benthic community in the Comau Fjord were significantly associated with both horizontal and vertical abiotic environmental gradients, determined by the influence of freshwater discharges. Species richness, diversity, and evenness peaked under high-salinity conditions, whose positions in depth varied across the horizonal axis of the fjord. Community structure was associated with a superficial water layer of low salinity. Dense patches of suspension feeders, or marine animal forests, may be contributing to the maintenance of diversity through facilitative interactions. Understanding the long-term impact of climate change-induced changes in salinity on the diversity and structure of fjord communities is necessary to generate sustainable management strategies for these ecosystems. Finally, an empirical and permanent assessment of the ecological impacts of anthropogenic sources of sedimentation is necessary to protect current communities and prevent future losses in ecosystem functioning.

Data Availability Statement

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

Author Contributions

VH and GF conceived and designed the study. JW provided the samples. SB and JE collected the samples. VV, KB-A, and SB processed the samples. VV and NV analyzed, interpreted data, and wrote the first draft. All authors reviewed and approved the text.

Funding

The study was financially supported by the FONDECYT #1190529 (NV), FONDECYT #1161699 (VH), and CONICYT/ANID #20150106 (VH and JW). NV was financially supported by FONDAP grant #15150003 (IDEAL). This work is part of VV Honors thesis, submitted to the Universidad Austral de Chile.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank all the members of the Huinay Field and Research team and PISCES project team who realized the sampling and fieldwork, and to all the members of the Coastal Ecology Lab (Lafkenchelab, UACH) who provided considerable help during the processing and analyzing period. José Garcés-Vargas was fundamental to complete the environmental analysis and visualization. Constructive criticism by Luis Miguel Pardo and Humberto González greatly improved the early research. This is publication #187 of Huinay Scientific Field Station.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.635855/full#supplementary-material

References

Aguayo, R., León-Muñoz, J., Vargas-Baecheler, J., Montecinos, A., Garreaud, R., Urbina, M., et al. (2019). The glass half-empty: climate change drives lower freshwater input in the coastal system of the Chilean Northern Patagonia. Clim. Change 155, 417–435. doi: 10.1007/s10584-019-02495-6