Determining the Ecological Status of Benthic Coastal Communities: A Case in an Anthropized Sub-Arctic Area

Dreujou, Elliot; Desroy, Nicolas; Carrière, Julie; Tréau de Coeli, Lisa; McKindsey, Christopher W.; Archambault, Philippe

doi:10.3389/fmars.2021.637546

ORIGINAL RESEARCH article

Front. Mar. Sci., 31 March 2021

Sec. Marine Evolutionary Biology, Biogeography and Species Diversity

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

This article is part of the Research Topic Biodiversity and Distribution of Benthic Invertebrates - From Taxonomy to Ecological Patterns and Global Processes View all 16 articles

Determining the Ecological Status of Benthic Coastal Communities: A Case in an Anthropized Sub-Arctic Area

$\r\nElliot Dreujou,,*$ Elliot Dreujou^1,2,3*

Nicolas Desroy⁴

Julie Carrière^5,6

Lisa Tréau de Coeli^2,3,6

Christopher W. McKindsey^2,7

Philippe Archambault^2,3,6

¹Institut des Sciences de la Mer, Université du Québec à Rimouski, Rimouski, QC, Canada
²Québec-Océan, Université Laval, Québec, QC, Canada
³Takuvik Joint Université Laval/Centre National de la Recherche Scientifique Laboratory, Université Laval, Québec, QC, Canada
⁴Laboratoire Environnement Ressources Bretagne Nord, Institut Français pour la Recherche et l’Exploitation de la Mer, Dinard, France
⁵Institut Nordique de Recherche en Environnement et en Santé au Travail, Sept-Îles, QC, Canada
⁶Département de Biologie, Université Laval, Québec, QC, Canada
⁷Maurice Lamontagne Institute, Fisheries and Oceans Canada, Mont-Joli, QC, Canada

With the widespread influence of human activities on marine ecosystems, evaluation of ecological status provides valuable information for conservation initiatives and sustainable development. To this end, many environmental indicators have been developed worldwide and there is a growing need to evaluate their performance by calculating ecological status in a wide range of ecosystems at multiple spatial and temporal scales. This study calculated and contrasted sixteen indicators of ecological status from three methodological categories: abundance measures, diversity parameters and characteristic species. This selection was applied to coastal benthic ecosystems at Sept-Îles (Québec, Canada), an important industrial harbor area in the Gulf of St. Lawrence, and related to habitat parameters (organic matter, grain size fractions, and heavy metal concentrations). Nearly all indicators highlighted a generally good ecological status in the study area, where communities presented an unperturbed profile with high taxa and functional diversities and without the dominance of opportunistic taxa. Some correlations with habitat parameters were detected, especially with heavy metals, and bootstrap analyses indicated quite robust results. This study provides valuable information on the application of environmental indicators in Canadian coastal ecosystems, along with insights on their use for environmental assessments.

Introduction

Anthropogenic influences on marine ecosystems occur globally, with possible perturbation of habitats and communities (Halpern et al., 2007, 2019). Many international organizations have recognized the importance of biologically diverse ecosystems for humanity and have established objectives and targets for their protection and sustainable use (United Nations, 1992; Secretariat of the CBD, 2010; SDG, 2015). The management of ecosystems requires an understanding of how habitats and communities respond to drivers of change, i.e., forces that affect environmental processes and modify ecosystem state from equilibrium (Boonstra et al., 2015; Beauchesne et al., 2020; Orr et al., 2020). In addition to natural drivers (e.g., temperature anomalies, freshwater inputs, hypoxic events), influences from human activities (e.g., fisheries, chemical pollution, species introductions) are also considered as ecosystem drivers. As natural and anthropogenic drivers may affect ecosystems concomitantly, it is important to understand how both relate to observed effects (Brown et al., 2014). To tackle these questions, environmental assessments rely on the best available knowledge, acquired through ecological groundwork in ecosystems of interest (such as biodiversity surveys, time series monitoring or experimental studies), and on the communication of results to a wide range of stakeholders (Borja et al., 2012; Borja, 2014; Chapman, 2016; Teixeira et al., 2016). Because such assessments are important foundations for decision makers, it is essential to properly account for the inherent complexity and variability of ecological data.

The use of integrative methods, such as indicators, is particularly relevant in this context. An indicator of ecological status is defined as a quantitative measure that synthesizes ecosystem information to infer ecosystem status (Rice, 2003; Rees et al., 2008). Many holistic frameworks, such as ecosystem-based management, marine spatial planning and DPSIR (Driver Pressure State Impact Responses) models, have included indicators in their methodology (Smeets and Weterings, 1999; Niemi and McDonald, 2004; European Commission, 2008; Rees et al., 2008; Levin et al., 2009; Atkins et al., 2011; Borja et al., 2013, 2015, 2016; Santos et al., 2019). However, environmental indicators evaluate specific ecosystem components, perturbations and/or spatiotemporal scales, potentially limiting their applicability in other systems, thus leading to the development of many indicators worldwide (Niemi and McDonald, 2004; Pinto et al., 2009; Teixeira et al., 2016).

One of the ecosystem components most frequently selected for environmental indicators are macrobenthic invertebrates, as they play an important role in the structure and functioning of benthic marine ecosystems (Dauvin and Ruellet, 2007; Pratt et al., 2014). Examples of this include engineering species (e.g., structural features for other species, bioturbation) and interactions with nutrient cycles (e.g., nutrient sequestration in sediments, remineralization, benthic-pelagic coupling) (Largaespada et al., 2012; Link et al., 2013; Belley et al., 2016; Bourque and Demopoulos, 2018). Many macrobenthic species are characterized by a sedentary lifestyle and a relatively long life span, which is particularly interesting when studying human influence as communities will reflect medium-term conditions, resulting in adaptation or local extinction (e.g., Dauer, 1993; Borja et al., 2000; Wei et al., 2020).

As pointed out by Rice (2003) and Salas et al. (2006), environmental indicators may be classed into categories according to their methodological basis, including three main categories used in environmental assessments. Category 1 regroups indicators based on measures of abundance—such as density and biomass of individuals—to infer community status. Relationships between abundance and a community status have frequently been discussed, as species do not have the same tolerance to disturbance (Pearson and Rosenberg, 1978). As such, the use of abundance-biomass curves has been proposed to detect if communities are in a balanced state, where K-selected taxa are dominant, compared to a disturbed state, with a dominance of r-selected taxa (Pearson and Rosenberg, 1978; Gray, 1979; Warwick and Clarke, 1994). Category 2 indicators are biodiversity parameters, i.e., community characteristics such as taxa identity and prevalence, which allow complex information to be aggregated into a unique metric. Finally, indicators in Category 3 are computed based on variations of responses of taxa to disturbance. Pioneer works by Pearson and Rosenberg (1978) proposed a model of benthic community evolution along a gradient of organic enrichment, laying the path toward a set of indicators that relate community structure and ecological status.

Environmental indicators, such as the AZTI Marine Biotic Index or the Infaunal Trophic index, have been applied in a number of North American ecosystems, including Chesapeake Bay, Willapa Bay and the Southern California coast (United States), but efficiency to detect perturbation has been mixed (Word, 1978; Maurer et al., 1999; Ferraro and Cole, 2004; Borja et al., 2008b; Pelletier et al., 2018). Less commonly, studies on the Pacific and Atlantic coasts of Canada have also evaluated the utility of existing indicators, although these studies have most often found poor performance (Sutherland et al., 2007; Burd et al., 2008; Callier et al., 2008; Robert et al., 2013). There is thus a need to test and validate indicators for Canadian ecosystems, in particular by comparing outcomes and efficiency of existing methods, which will greatly benefit to tackle ecosystem management objectives within Canada’s Ocean Act and the Oceans Strategy (Government of Canada, 1996; Department of Fisheries and Oceans, 2002).

To this end, we evaluated various indicators of ecological status in a coastal industrial harbor area, where human activities may significantly impact local benthic ecosystems. Industrial harbor areas are regions regrouping significant industrial activities coupled with harbor platforms linking production with commercial shipping routes worldwide. We selected the region of Sept-Îles (Québec, Canada) for this study. Located in the Gulf of St. Lawrence, one of the management areas designated by Fisheries and Oceans Canada and a major strategic region for Québec (Department of Fisheries and Oceans, 2009; Daigle et al., 2017; Schloss et al., 2017; Ferrario and Archambault, in preparation), Sept-Îles is the fourth largest Canadian port in 2019 in terms of total exchanged goods and the second largest in Québec (Statistics Canada, 2011; Binkley, 2020). Industrial activities at Sept-Îles are largely focused on international shipping of iron ore mined in northern Québec and Labrador, the production of aluminum and various fisheries operate in the bay (Department of Fisheries and Oceans, 2019).

The objectives of this study are to (i) compare outcomes of various environmental indicators on benthic ecosystems of the Sept-Îles region and (ii) understand how these indicators relate to habitat parameters for validation and to select appropriate applications.

Materials and Methods

Study Area

We targeted ecosystems with a sandy-silty sediment in the industrial harbor area of Sept-Îles (Côte-Nord region of Québec, Canada), which considers ecosystems in the Baie des Sept Îles and the archipelago at its entrance (Figure 1A; Dreujou et al., 2018, 2020). Coasts are characterized by sandy beaches, tidal marshes and anthropogenic structures. Mean depth is 35 m in the bay and can reach up to 150 m in the archipelago (Dutil et al., 2012). It is influenced by freshwater inputs from multiple streams and strong tidal currents resulting in a mixed water column and an estuarine circulation (Shaw, 2019). Ecosystems in the Sept-Îles region are considered sub-Arctic due to the formation of ice on the shore in November/December and in the bay in January/February, along with an important freshwater run-off due to snowmelt in April (Demers et al., 2018).

FIGURE 1

Figure 1. Map of the study area. (A) Location of the sampled stations, with light blue triangles and dark blue squares representing shallow (< 15 m) and deep (> 15 m) stations, respectively, (B) Location and identity of human activities present in the area.

This region hosts several human activities, including industrial, commercial and dredging operations located at the City of Sept-Îles and the Pointe-Noire sector (on the southern section of Baie des Sept Îles), along with an aquaculture site and various fisheries throughout the bay (Figure 1B). Many projects have been done in this region to characterize pelagic and benthic communities and habitats in relation to coastal stressors (Canadian Healthy Oceans Network, 2016; Carrière, 2018; Dreujou et al., 2020).

Benthic Ecosystems Sampling

The sampling design and methods used to collect and analyze ecological samples were similar to those presented in Dreujou et al. (2020), with the exception that only one region (Baie des Sept Îles) and one type of community (individuals higher than 0.5 mm) were considered.

A total of 108 stations were selected in the study area, using a randomization algorithm to cover the full extent of the sector, constrained between 0 and 80 m deep, and with increased sampling effort in areas with human activities (Figure 1A). Himmelman (1991) showed that benthic communities in the Northern Gulf of St. Lawrence above and below 15–20 m deep differ. Likewise, preliminary fieldwork in the study region detected a thermocline in the water column at ca. 15 m deep. Consequently, we discriminated two groups of stations in order to ensure habitat homogeneity within depth classes: shallow (<15 m, 26 stations) and deep habitats (>15 m, 82 stations). We sampled the benthic ecosystem in July 2017, using a Ponar grab (0.05 m²) deployed from a boat, with two independent casts at each station.

The first cast collected two subsamples—one for the analyses of organic matter content and another for sediment grain size—stored at -20°C until processing in the laboratory. The percentage of total organic matter (i.e., sum of organic carbon and organic nitrogen) in the sediment was determined using the Loss-on-Ignition method (Davies, 1974). Grain-size analysis was done on a sieving column for the fraction with particles larger than 2 mm and with a Laser Diffraction Particle Size Analyzer for the smaller fractions. Results from both techniques were combined to yield a unified size distribution range from 0.04 μm to 26.5 mm. From this, percentages of gravel, sand, silt and clay were calculated as defined by Wentworth (1922) and Folk (1980).

All sediment obtained from the second cast was sieved on a 0.5 mm mesh size and preserved in a solution of BORAX-buffered formalin (4%) solution for subsequent benthic macrofauna identification (Dreujou et al., 2020). The resulting samples were sorted using a stereomicroscope and taxa identified to the lowest taxonomic level possible with reference manuals and identification guides; names were validated according to the World Register of Marine Species (WoRMS Editorial Board, 2020). Taxon density and biomass per grab were recorded by counting and weighting (blotted wet mass) individuals in each sample, respectively.

In addition to these parameters, we considered estimates of heavy metal concentrations in the sediment. Concentrations at the sampled stations were calculated based on values obtained in the same area in 2014 and 2016, retrieved from a database hosted by Carrière (2018), using Inverse Distance Weighting interpolation (Dale and Fortin, 2014). We focused on metals for which toxicity criteria have been defined in the Biological Effects Database for Sediments (Environment Canada and Ministère du Développement Durable de l’environnement et des Parcs du Québec, 2007; Centre d’Expertise en Analyse Environnementale du Québec, 2014): arsenic, cadmium, chromium, copper, mercury, lead and zinc; we also included iron and manganese to account for possible contamination from local ore industries.

Environmental Indicator Calculation

Indicators of ecological status were selected from Pinto et al. (2009), DEVOTES (2012), and Teixeira et al. (2016), and grouped into three Categories according to their methodology (Table 1). We targeted indicators related to descriptors D1 (biological diversity), D6 (seafloor integrity), and D8 (contaminants) of Good Environmental Status (European Commission, 2008; Borja et al., 2013), choosing those that applied to benthic invertebrates in soft-bottom habitats. We considered each station separately, allowing an assessment of the spatial variability and mean for each indicator, and when possible we pooled all stations together to obtain an estimate for the bay-scale system. We used R v4.0 to perform data manipulations and calculations (R Core Team, 2020).

TABLE 1

Table 1. Summary of the evaluated indicators.

We included in Category 1 the total density (number of individuals collected per grab), total biomass (wet mass of individuals collected per grab), and the W-Statistic Index, calculated based on abundance-biomass curves for the community (Warwick and Clarke, 1994). Those indicators were computed using benthic taxa abundance sampled at each station.

For Category 2, we considered taxa richness (number of collected taxa) and related metrics to describe the community’s structure and the relative prevalence of taxa within it, such as the Shannon index, Margalef index, Simpson index, and Pielou evenness (Legendre and Legendre, 1998; Magurran and McGill, 2011). We also considered taxonomic and functional diversities, based on taxonomic relationships between taxa and information about biological traits, respectively (Warwick and Clarke, 1995; Clarke and Warwick, 1998; Mason et al., 2005; Villéger et al., 2008). Taxa richness, Shannon index, Margalef index, Simpson index, and Pielou evenness were calculated using the benthic community at each station. For taxonomic diversity, we gathered relatedness data for taxa using the WoRMS online database (WoRMS Editorial Board, 2020). To estimate functional diversity, we computed functional richness, functional evenness and functional divergence (Mason et al., 2005; Villéger et al., 2008) by considering five biological traits—body composition, body size, feeding type, mobility and lifestyle—with a total of 26 modalities (Table 2). Because taxa can present several modalities for a trait, we assigned a continuous value between 0 (absence of the modality) and 1 (presence of the modality) for each taxon and each trait (the sum of values for every modality within a trait equals 1). Biological trait data was extracted from WoRMS, SealifeBase, the Encyclopedia of Life, and Arctic Traits databases as well as dedicated articles (Degen and Faulwetter, 2019; EoL, 2020; Palomares and Pauly, 2020; WoRMS Editorial Board, 2020). R Packages vegan and FD were used to calculate indicators in this category (Laliberté and Legendre, 2010; Laliberté et al., 2014; Oksanen et al., 2019).

TABLE 2

Table 2. Summary of the functional traits and modalities.

Finally, indicators in Category 3 included the AZTI Marine Biotic Index (AMBI) and its multivariate version (M-AMBI), which are based on the relative proportion of taxa classified into five ecological groups depending on their tolerance to perturbation (Grall and Glemarec, 1997; Borja et al., 2000; Muxika et al., 2007), BENTIX, where only two ecological groups are considered (Simboura and Zenetos, 2002), and the Benthic Opportunistic Polychaetes Amphipods Index (BOPA), which compares proportions of opportunistic polychaetes and amphipods (Dauvin and Ruellet, 2007). Sampled taxa were assigned to ecological groups, from group I to V, based on the list of Borja et al., version of May 2019 (AZTI, 2019; Supplementary Table S1). Because this list was developed for European taxa, we assigned groups to unregistered taxa based on species physiology studies and taxonomic relationships (Pelletier et al., 2018). We used this list to further regroup taxa to a “sensitive” (groups I and II) and a “tolerant” (groups III to V) metagroup to compute BENTIX (Simboura and Zenetos, 2002), and to obtain the proportion of opportunistic polychaetes (groups III to V) and sensitive amphipods (group I) to calculate BOPA (Dauvin and Ruellet, 2007; Supplementary Table S1). M-AMBI was calculated using the dedicated software AMBI v5.0 (AZTI, 2019), where “bad” and “high” status conditions are required for taxa richness, Shannon index and AMBI (Muxika et al., 2007). Because historical data on benthic invertebrates is scarce in our study area, we used the outcomes of our sampling to establish these values by selecting the 5 and 95 percentiles of the variable distribution (for “bad” and “high” status, respectively, Supplementary Table S2) (Buchet, 2010).

Integration and Statistical Analysis

Results for each indicator were reviewed qualitatively and compared to benthic ecosystem data in the Gulf of St. Lawrence, when available. Robustness for indicators in Categories 1 and 2 was calculated as the 95% confidence interval using a resampling routine (bootstrap, 1000 replicates), and the difference between averages of each indicator and the resampling averages (i.e., bootstrap bias).

We computed the Ecological Quality Ratios for Category 3 indicators. This ratio compares the value of an indicator to a reference, such as a targeted state or unperturbed/pristine ecosystem, so that an Ecological Quality Status can be assigned (five categories: “bad,” “low,” “moderate,” “good,” and “high” status). The formula to compute the Ecological Quality Ratio is the following (Bund and Solimini, 2007):

E Q R = \frac{V_{i n d} - R_{b a d}}{R_{h i g h} - R_{b a d}}

V_ind is the value of an indicator, R_bad is the reference value for a “bad” status and R_high is the reference value for a “high” status. Limits between each Ecological Quality Status class are specific to the indicator used (Borja et al., 2000; Simboura and Zenetos, 2002; Muxika et al., 2005, 2007; Dauvin and Ruellet, 2007).

Finally, we explored covariation between indicators and habitat parameters (organic matter content, grain size distribution and heavy metal concentrations), using scatterplots for each pair of variables. Correlation was assessed with Spearman’s rank coefficients to understand the relevance of each indicator to the computation of ecological status (Quinn and Keough, 2002).

Results

Overview of Benthic Habitats and Communities

Sediment was mostly composed of sand and silt fractions, with concentrations of organic matter rarely surpassing 3% (Supplementary Table S3). Heavy metal concentrations did not reach high toxicity levels as defined by Environment Canada (Environment Canada and Ministère du Développement Durable de l’environnement et des Parcs du Québec, 2007; Centre d’Expertise en Analyse Environnementale du Québec, 2014; Dreujou et al., 2020; Supplementary Table S3). A total of 132 taxa were identified, belonging to eight phyla, with a dominance of arthropods, mollusks, and annelids (Supplementary Table S1). The most abundant taxa were the polychaete Micronephthys neotena, the cumacean Eudorellopsis integra, the amphipod Protomedeia grandimana, Nematoda (adults), and the bivalve Macoma calcarea (Supplementary Table S1). From this list, no species which can be considered as exotic to this region have been reported (Simard et al., 2013).

Indicator Outcomes

Category 1 Indicators

Indicators in this category presented greater mean values in deep than shallow stations, with the exception of total density (Table 3). Shallow stations showed a higher total density than deep stations, but this may be an outlier effect due to a single station close to the City of Sept-Îles (Supplementary Figures S1A–C), where density was 899 individuals.grab^–1 with a dominance of P. grandimana. Overall, shallow and deep stations presented low total biomass, except for a couple of stations due to the presence of the echinoderms Echinarachnius parma and Strongylocentrotus sp. (Supplementary Figures S1A–C). The W-Statistic Index was positive and close to zero at nearly all shallow and deep stations (Supplementary Figures S1A–C) and the abundance-biomass curve presented higher abundance than biomass values when species were ranked (Figure 2).

TABLE 3

Table 3. Values of the mean and standard error (SE) for each indicator, the difference between bootstrapped mean and the true mean (bias) and the 95% confidence interval (CI), for shallow and deep stations.

FIGURE 2

Figure 2. Values of abundance and biomass for ranked species (logarithm scale), for shallow stations and deep stations.

Category 2 Indicators

Category 2 indicators showed similar trends for shallow and deep stations, while being generally higher for the latter (Table 3). In particular, there is a close similarity between the spatial distributions of taxa richness, Shannon and Margalef indices and taxonomic diversity (Supplementary Figures S1D–L). Variability for shallow stations is quite low, except for a station in front of Pointe-Noire where only one taxon was present, while deep stations tend to display the highest values in the archipelago compared to the center of the bay (Supplementary Figures S1D–L). Mean values for the Simpson index and Pielou evenness reached 0.62 (standard error of 0.04) and 0.77 (0.02), respectively, for shallow stations and 0.66 (0.05) and 0.76 (0.02) for deep stations (Table 3). The same relationship between shallow and deep stations is observed for these metrics, even though the distribution for both is skewed with some stations closer to coasts presenting very low values (Supplementary Figures S1D–L). Concerning functional diversity, deep stations presented higher mean functional richness, functional evenness and functional divergence relative to those at shallow stations (Table 3). The most abundant modality for each biological trait was non-calcified tissue for body composition, small individuals for body size, surface deposit-feeders for feeding type, mobile organisms for mobility and burrowers for lifestyle, at both shallow and deep stations.

Category 3 Indicators

Classification of taxa into ecological groups to compute Category 3 indicators yielded 51 taxa in group I (sensitive to disturbance, 38.6% of the taxa), 63 in group II (indifferent to disturbance, 47.7%), 11 in group III (tolerant to disturbance, 8.3%), 1 in each of groups IV and V (second- and first-order opportunists, respectively, 0.8%) and 5 were not assigned due to a too broad taxonomic resolution (Supplementary Table S1). This classified 114 taxa in the “sensitive” group and 13 in the “tolerant” group (Supplementary Table S1). Concerning polychaetes and amphipods, we observed four opportunistic polychaetes (Cossura longocirrata, Eteone sp., Hediste diversicolor, Praxillella praetermissa) and nine sensitive amphipods (Ameroculodes edwardsi, Ampelisca vadorum, Byblis gaimardii, Lysianassidae, Maera danae, Phoxocephalus holbolli, Pontoporeia femorata, Quasimelita formosa, Quasimelita quadrispinosa).

An AMBI score of 1.57 and 1.53 was obtained for the bay-scale estimate at shallow and deep stations, respectively, which corresponds to a “slight imbalance” site classification (Borja et al., 2000). Overall, low AMBI values were obtained at each station, being 1.5 on average (standard error of 0.13) for shallow stations and 1.45 (0.05) for deep stations, and never exceeding 3, and no particular spatial trend can be observed (Table 3 and Supplementary Figures S1M–P). The bay-scale M-AMBI could not be computed with the percentile method, and at the station level, generally high mean values of 0.68 (0.05) and 0.7 (0.03) were observed for shallow and deep stations, respectively (Table 3). Stations outside of the bay tended to be characterized by higher values than those inside it, especially close to the coast and in the northern section of the bay, but this may be related to the spatial distribution of taxa richness and the Shannon index (Supplementary Figures S1M–P). The BENTIX bay-scale estimate was 5.15 for shallow stations and 5.25 for deep stations, while at the station-level mean values were 4.95 (0.23) and 5.31 (0.09), respectively (Table 3). These values correspond to a “normal/pristine” pollution classification for the majority of the area sampled, except for some stations close to coasts (Simboura and Zenetos, 2002). Finally, BOPA produced low scores of 0.002 and 0.004 for shallow and deep bay-scale estimates, respectively, similar to means of 0.0028 (0.0012) for shallow and 0.0067 (0.003) for deep stations, respectively (Table 3), denoting “high status” classifications. Only two stations had a score higher than 0.05, a trend that is not shared with neighboring stations, which may indicate localized low-intensity perturbations (Supplementary Figures S1M–P).

Calculation of Ecological Quality Ratios using Category 3 indicators produced similar results for AMBI, BENTIX and BOPA (Figure 3). The majority of stations (shallow and deep) presented a “high” or “good” ecological status except for a few stations with a “poor” status (Figure 3). In contrast, results for M-AMBI were less uniform, with a high variation among both shallow and deep stations, such that no general trends may be highlighted (Figure 3).

FIGURE 3

Figure 3. Values of Category 3 indicators ranked according to Ecological Quality Ratios, calculated for shallow and deep stations. (A) Calculated with the AZTI Marine Biotic Index (AMBI), (B) Calculated with the Multivariate AZTI Marine Biotic Index (M-AMBI), (C) Calculated with the BENTIX, (D) Calculated with the Benthic Opportunistic Polychaetes Amphipods Index (BOPA). B = “bad” status (red), P = “poor” status (orange), M = “moderate” status (yellow), G = “good” status (green), H = “high” status (blue).

Robustness and Covariation

For Category 1 and 2 indicators, bootstrap bias was low at both shallow and deep stations (less than 0.4), except for functional richness where it reached 3.17 and 7.59, respectively (Table 3), demonstrating a relatively high robustness of the indicators. The true mean was included in the 95% confidence interval for five indicators at shallow stations (taxa richness, total density, total biomass, functional evenness, functional divergence) and eight at deep stations (taxa richness, total density, total biomass, Shannon index, Margalef index, Simpson index, Pielou evenness, taxonomic diversity) (Table 3).

The analysis of covariation between indicators reported moderate to very high Spearman’s coefficients (0.22 < |ρ| < 0.96) (Table 4). Category 2 indicators presented the highest proportion of within-Category significant correlations at both shallow and deep stations (Table 4). The vast majority of these correlations were positive, with the strongest correlations between Shannon and Margalef indices, and were represented by linear proportionality between indicators on the scatterplots. Category 2 indicators were also frequently correlated to indicators from Categories 1 and 3, especially for the W-Statistic Index and the M-AMBI (Table 4). The latter Categories did not present high within-Category correlations, except between AMBI/BENTIX and M-AMBI/BOPA at shallow stations, and the W-Statistic Index and AMBI at deep stations.

TABLE 4

Table 4. Spearman rank correlation coefficients between indicators, for shallow and deep stations.

Relationships With Habitat Parameters

Correlations between Category 1 indicators and abiotic parameters detected non-significant relationships with sediment parameters (except between the W-Statistic Index and gravel and sand contents at deep stations), while they were significant and negative between most heavy metals and total density and total biomass at shallow stations, and the W-Statistic Index at deep stations (Table 5). The absolute value of Spearman’s rank coefficients was high for total density and total biomass at shallow stations (between −0.4 and −0.61), highlighting relatively strong relationships, while they were less for the W-Statistic Index at deep stations (between −0.22 and −0.29).

TABLE 5

Table 5. Spearman rank correlation coefficients between environmental indicators and habitat parameters, for shallow and deep stations.

For Category 2 indicators, correlations with sediment parameters were significant only for some cases involving taxa richness, the Margalef index, taxonomic diversity and functional richness (Table 5). Relationships with heavy metals were detected mainly at deep stations, in particular for cadmium, copper, lead and zinc; at shallow stations, functional richness showed significant correlations with all heavy metals except cadmium, while functional divergence and taxa richness presented marginal correlations. The vast majority of these relationships were moderate to high (between −0.22 and −0.45), except at deep stations for gravel and sand contents and between functional divergence and some heavy metals.

Finally, several significant relationships were observed between Category 3 indicators and sediment parameters (organic matter, sand and silt contents), including at shallow stations for AMBI and BENTIX and at deep stations for BENTIX and BOPA (Table 5). Organic matter was negatively correlated with AMBI values (coefficient of −0.43) at shallow stations and positively with BENTIX values at shallow and deep stations (0.45 and 0.27, respectively); sand and silt contents had the opposite effect at shallow stations for AMBI (0.47 and −0.47, respectively) and at deep stations for BENTIX (−0.26 and 0.23, respectively) and BOPA (−0.31 and 0.34, respectively) values. Many relationships with heavy metals were detected at deep stations for all indicators except AMBI (Table 5). In particular, M-AMBI presented negative correlations with heavy metals (between −0.24 and −0.38), whereas correlations with BENTIX and BOPA were positive (between 0.23 and 0.36).

Discussion

Strengths, Limitations, and Ecological Considerations of Indicators

The analysis of benthic communities using Category 1 indicators relies on abundance relationships (either density or biomass of individuals) without consideration of taxonomic identity. Their calculation requires the least laboratory and analytical time relative to the other calculated indicators. Deep stations present a higher density of benthic organisms than shallow stations, as predicted by patterns of coastal marine biodiversity (Gray and Elliott, 2009; Levinton, 2013; Piacenza et al., 2015). The abundance-biomass curve for shallow and deep stations is characteristic of an unstressed profile (Pearson and Rosenberg, 1978; Warwick and Clarke, 1994), which is further supported by the W-Statistic Index being positive and close to 0 at nearly all stations (Clarke, 1990). Studying communities through abundance relationships thus provides interesting results concerning the status of the ecosystem, but the main assumption behind indicators in this Category is that all species are equivalent and have an identical role in the ecosystem structure and functioning. This, however is not necessarily true, as some species can be considered “key species” in ecosystems due to unique engineering or trophic roles (e.g., Bond, 1994; Lawton and Jones, 1995). Thus, Category 1 indicators should be coupled with ancillary methods focusing on biological characteristics of the species, such as life-history traits and physiological characteristics.

Category 2 indicators focus on community biodiversity, granting additional detail than that provided by Category 1 indicators. The notion of biodiversity can be interpreted along multiple points of view in an ecosystem, such as the diversity of species, genes, habitats or functions (United Nations, 1992; Wilson, 1992; Hooper et al., 2005; Stachowicz et al., 2007). While each targeted component has specific implications for the ecosystem, high richness and high diversity values have generally been interpreted as signs of good ecological status (Covich et al., 2004; Borja et al., 2013). This statement needs to be considered carefully, as it is necessary to discuss results with comparable ecosystems and historical data so that diversity trends are interpreted according to local background patterns (Covich et al., 2004).

Taxa richness indicated a quite diverse community, being nearly twice as great at deep stations, as expected given general trends (Gray and Elliott, 2009; Levinton, 2013; Piacenza et al., 2015). These numbers (132 taxa observed) are comparable to results from available benthic invertebrate surveys done in the study area, such as 27 taxa reported by OBIS (2020) in the Sept-Îles region and 148 taxa by Nozères et al. (2015) in the Gulf of St. Lawrence. Results from diversity indicators (Shannon index, Margalef index, Simpson index, Pielou evenness and taxonomic diversity) showed moderate to high benthic diversity, with no dominance by any taxa (even distribution) and great taxonomic breadth. Few stations differ from this general trend, and those that do are mostly close to coasts where diversity is low and there is no clear evidence of perturbation. Diversity indicators are frequently used in ecological studies to characterize communities and to detect disturbance, which allows discussing their results building on a vast corpus of studies worldwide (Magurran and McGill, 2011). These univariate estimates, however, may mask individual responses arising for example from adaptation or changes in biotic relationships, and they should be coupled with multivariate methods, such as ordination or similarity analysis, so that community-level effects are accurately described (Legendre and Legendre, 1998; Quinn and Keough, 2002; Magurran and McGill, 2011).

Concerning functional diversity, functional richness is generally lower at shallow stations and the two other indicators are in the same range, albeit being slightly greater at deep stations. These results suggest that taxa at shallow stations have more specialized niches, i.e., less diverse functional strategies (Villéger et al., 2008), indicating some redundancy of biological traits. This property is linked to an increased ecosystem stability and resilience, where possible extinctions due to perturbation will not modify the ecosystem structure even if some taxa disappear (Rosenfeld, 2002; Mouillot et al., 2013). However, bootstrap bias was very high for this indicator, making conclusions less robust. Moderate to high functional evenness and divergence denote that values for given biological traits are not evenly distributed and are skewed toward extremes (0 or 1). The consideration of biological traits in addition to species identity allows to study functions of the ecosystem, such as productivity, chemical elements cycling or energy transfers (Somerfield, 2008; Bellwood et al., 2019). This is an important addition to environmental assessments, as biological traits and adaptative responses to disturbance are highly related (Mouillot et al., 2013; Miatta et al., 2021). Indicators of functional diversity thus offer valuable information to characterize communities in complement to other Category 2 indicators, at the expense of increased analytical time to assemble a traits database, for which information may be lacking or difficult to obtain.

For Category 3 indicators, community ecological status is assessed by considering the tolerance of taxa to perturbation. Values for these indicators highlight an overall high status in the study area, where taxa sensitive to perturbation are present without a dominance of opportunists, as illustrated by the number of stations with a “high” Ecological Quality Status. M-AMBI detected greater variability between stations relative to the other indicators, particularly within the bay. A possible interpretation for this result may be the influence of the percentile method used to compute “bad” and “high” status conditions for this indicator, advocating for a careful description of comparison conditions and a fortiori references values (Borja et al., 2012). Furthermore, classification into ecological groups may introduce bias, as the list of Borja et al. (2000) were primarily designed for European coastal ecosystems. Inclusion of taxa found on Canadian coasts has been made based on taxonomic similarity with species already included in the list, by reviewing studies on perturbation tolerance and by expert opinion, but these choices need to be ground-truthed by dedicated ecological works. A wide spectrum of indicators may be included within Category 3, as compiled by Pinto et al. (2009) and Teixeira et al. (2016), and their use in environmental assessments is linked to scientific and management objectives. Such indicators have been used with success in a variety of ecosystems (e.g., Borja et al., 2008a; Gillett et al., 2015), but they require a high volume of data to accurately relate taxa and perturbation status, such as field observations, modeling of species distributions (native and non-indigenous), physiological studies and experimental work.

With these results, we can compare strengths and limitations of the calculated indicators. Category 1 indicators gather relevant baseline information on the ecosystem and requires the least time to be computed. The downside is that it is difficult to discriminate between anthropogenic perturbation and natural variability as other community characteristics may be impacted and, most importantly, they cannot be compared to reference conditions of ecological status. Category 2 indicators, such as the commonly used Shannon index or Pielou evenness, are easy to compute from well-built taxa lists, although taxonomic and functional diversity demand more time to gather complementary information about phylogenetic relationships and ecosystem functions. However, the latter indicators provide more information on the community structure and are backed by ecological literature to infer a certain ecological status (Magurran and McGill, 2011). Finally, Category 3 indicators demand the most time to be calculated, in particular with the classification of taxa relative to their response to disturbance, but they have been specifically designed to determine Ecological Quality Ratios and to consider reference conditions. Bias or uncertainty may be introduced during the classification process as extensive experimental groundwork is needed to properly assign taxa to groups, which is not always available. Many of these indicators are region-specific, with possible poor performance in other ecosystems (e.g., Callier et al., 2008; Robert et al., 2013), so that further research is needed to properly assess ecological status in sub-Arctic regions.

Implications for Sept-Îles and Canadian Ecosystems

Sept-Îles is an important industrial harbor area for Québec, with a variety of economic activities taking place in the bay and archipelago. All calculated indicators except M-AMBI pointed toward diverse benthic communities of generally good ecological status and no particular perturbation patterns have been detected in the study area, which is coherent with previous descriptions of benthic ecosystems in this region (Carrière, 2018; Dreujou et al., 2020). When applying Category 3 indicators on the data from these studies, we obtained similar conclusions with ecological status indicated as “good” to “high” for AMBI, BENTIX and BOPA.

No particular trend has been observed for stations identified as potentially impacted by Dreujou et al. (2020) in coastal regions close to the City of Sept-Îles and Pointe-Noire, further suggesting an overall limited effect of perturbations. Compared to the regions where many of the assayed indicators were developed, e.g., in Atlantic and Mediterranean European ecosystems, the magnitude of human activity is considerably lower at Sept-Îles. As such, it is possible the range of variation induced by anthropogenic perturbation is not sufficient to severely impact benthic ecosystems. This is even more relevant when we compare these results to other industrial harbor areas worldwide, where human influence is more pronounced (e.g., Hewitt et al., 2005; Borja et al., 2006; Chan et al., 2016; Birch et al., 2020). Other hypotheses may explain this, such as (i) high community resilience and resistance, (ii) limitation of effective impacts of activities by the dynamic of the ecosystem (e.g., flushing from tidal currents), and (iii) perturbation effects may be more pronounced on other components (such as phytoplankton or pelagic species).

Ecological indicators represent a valuable method to set conservation targets and to guide sustainable environmental projects. In light of initiatives such as the Marine Strategy Framework Directive, where indicators and descriptors have been identified to monitor the ecological status of European marine waters (European Commission, 2008; Borja et al., 2013, 2015, 2016), local stakeholders have the possibility to build on these works to establish ecosystem-based management adapted to Canadian ecosystems. Further research in other industrial coastal areas, including long-term monitoring, is needed to obtained coherent and robust environmental assessments.

Validation and Limitations

Assessing relationships between indicators highlighted correlations, especially among Category 2 indicators. While this does not necessarily imply causality in the interpretations, covariation indicates that information gathered by some indicators is similar. This was expected for indicators relying on specific ecosystem components to be computed, such as M-AMBI and the Shannon index, both of which being a function of taxa richness. With an environmental assessment perspective, these results show that calculation of some indicators will not provide additional information when the objective is to detect trends in a targeted area. Understanding these links will allow refining methodological protocols and to produce more efficient and accurate assessments.

The use of ecological status indicators requires a validation procedure to ensure that outcomes are relevant (Dauvin et al., 2010; Heink et al., 2016; Burgass et al., 2017; Moriarty et al., 2018). Because the region of Sept-Îles is not frequently represented in the scientific literature, it is then difficult to have baseline data to validate ecosystem assessments. The studies of Dreujou et al. (2018, 2020) represent the only campaigns describing benthic habitats and communities in the Baie des Sept Îles, and increased sampling would greatly improve indicator validation, in particular for Category 1 and 2 indicators where references conditions are currently not available. Inference of ecological status based on Category 3 indicators is relying on explicit reference conditions for “bad” and “high” status. Defining values for these conditions based on contemporary ecosystems will most certainly introduce bias, as most are likely to show some level of degradation (which cannot be assessed), and alternatives, such as historical datasets, are rare (Muxika et al., 2007; Borja et al., 2012).

Overall, Category 1 and 2 indicators were relatively robust, with little difference between mean values calculated from the real and bootstrapped datasets (except for functional richness), indicating quite homogeneous results. The vast majority of significant correlations between indicators and environmental factors were found for heavy metal concentrations and most such correlations were negative. This implies that indicators would successfully detect perturbation due to heavy metal content, thus resulting in reduced ecological status, but fail to detect perturbations affecting other habitat parameters. AMBI and BENTIX were correlated to organic matter content, which was expected as original works of Borja et al. (2000) and Simboura and Zenetos (2002) were based on models predicting community changes in response to an organic enrichment (Pearson and Rosenberg, 1978; Grall and Glemarec, 1997). Concerning grain size variables, only sand and silt contents showed any significant correlations with indicators, mainly at deep stations, which may be due to very low amounts of gravel and clay in the sediments.

It is important to note that indicators summarize complex ecological data into unique values (univariate in a statistical point of view), which may be insufficient to correctly assess perturbation. Category 3 indicators were developed for specific types of disturbance, such as organic matter loading or oil-spill detection (Pearson and Rosenberg, 1978; Borja et al., 2000; Dauvin and Ruellet, 2007), and for specific ecosystems (e.g., European Commission, 2008). Even though we detected significant relationships with heavy metal concentrations, dedicated methods to monitor these types of perturbation would greatly benefit this portrait. Finally, the consideration of cumulative impacts from various sources of disturbance may be a good perspective for environmental assessments (Crain et al., 2008), in order to develop more holistic indicators.

Conclusion

This study provides insight on the use of environmental indicators in Canadian coastal ecosystems, in particular by applying and comparing indicators within an important industrial harbor area. An overall good status for benthic ecosystems was detected in our study area, which is a valuable addition to guide stakeholders and ecosystem management at the local scale. We were able to present strengths and caveats for each Categories of indicators, relating them to ecological implications and possible improvements. Finally, we were able to study the robustness of indicators and we highlighted improvements for ground-throughing and validation.

Further environmental assessments in sub-Arctic coastal areas in the Gulf of St. Lawrence and the Canadian Eastern Atlantic coast are needed to obtain a broader portrait of the region with a more diverse range of environmental conditions. Long-term monitoring will also produce reliable time series data to better understand variability of sub-Arctic benthic ecosystems.

Data Availability Statement

The datasets generated for this study can be found in the online repositories. The names of the repository/repositories and accession number(s) can be found below: https://dataverse.scholarsportal.info/dataverse/chone2, doi: 10.5683/SP2/WDDDMI.

Author Contributions

ED, JC, CM, and PA contributed to the conception and design of the study. ED gathered the data in the field. ED, ND, and LT organized the database. ED and ND performed the statistical analysis. ED wrote the first draft of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.

Funding

This research was sponsored by the NSERC Canadian Healthy Oceans Network and its partners: Fisheries and Oceans Canada and INREST (representing the Port of Sept-Îles and City of Sept-Îles) to PA (grant number: NETGP 468437-14), by the FRQNT Strategic Network Québec-Océan and by the Joint UL/CNRS Laboratory Takuvik.

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 would like to acknowledge and thank all the people that helped during the field campaigns, lab work, and data analysis: Claudy Dechêsne, Serge Gallienne, Dany Lévesque, Aurélie Foveau, Cindy Grant, Laure de Montety, Philippe-Olivier Dumais, and Raphaël Bouchard. This study is a contribution to the research program of FRQNT Strategic Network Québec-Océan, which also provided scientific support.

Supplementary Material

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

References

Atkins, J. P., Gregory, A. J., Burdon, D., and Elliott, M. (2011). Managing the marine environment: Is the DPSIR framework holistic enough? Syst. Res. Behav. Sci. 28, 497–508. doi: 10.1002/sres.1111