Tracking Pre- and Post-Industrialization Changes in the Bay of Sept-Îles, Canada, Using Foraminifera as Bioindicators

Coastal ecosystems are vital components of the Earth’s biosphere, but they face escalating threats from anthropogenic activities and climate change, rendering them among the most vulnerable environments globally (Williams et al., 2022). The interplay of human interventions and natural processes has substantially altered coastal landscapes, with repercussions across ecological, social, and economic dimensions. Anthropogenic stressors, ranging from industrial and agricultural pollutants to urbanization and habitat modification, exert pressures on coastal environments, compromising their functionality (Adams, 2005). Activities, including overuse of fertilizers and emission of toxic and nontoxic gases into the atmosphere, contribute to the deterioration of coastal ecosystem health and the loss of marine life, especially in areas where human settlements are dense and widespread (Allegra et al., 2018; Bailey et al., 2020). Additionally, anthropogenic noise from maritime operations, the presence of industrial facilities, and coastal development further exacerbate the ecological stressors confronting these sensitive ecosystems (Popper et al., 2014).

The Bay of Sept-Îles (BSI), located on the North Shore of Québec province, Canada, within the Gulf of St. Lawrence, is an example of a high-use coastal environment that faces multiple potential sources of contamination (Carrière and Le Hénaff, 2018). This region is subject to various anthropogenic pressures, including North America’s largest mineral port, the largest aluminum smelter in the Americas, and an urban settlement. Commercial fishing and tourism are also prominent, leading to sustained environmental exposure to anthropogenic disturbances (Dreujou et al., 2020).

Paleoecological studies offer invaluable insights into long-term environmental trends and variability (Nieto-Lugilde et al., 2021). By examining changes in ecosystems over extended periods, researchers can discern the drivers, extent, and consequences of environmental shifts (Taffs et al., 2017). Various environmental proxies, including elemental geochemistry, magnetic properties, and bioindicators, have been employed to assess anthropogenic impacts on coastal benthic environments (Dasgupta et al., 2018; Frontalini and Coccioni, 2008; Yang et al., 2007). Among these proxies, benthic foraminifera emerge as particularly effective bioindicators of environmental change, owing to their rapid response to environmental fluctuations and their ability to faithfully record variations in their surroundings (Boltovskoy, Scott, and Medioli, 1991; Weber and Casazza, 2006).

Foraminifera possess several characteristics that render them ideal environmental indicators, including short lifespans, rapid growth rates, diverse ecological preferences, and small size (Desrosiers et al., 2013). Moreover, they exhibit high population densities, allowing for statistically robust sample sizes to be collected quickly and inexpensively with minimal environmental disruption (Frontalini and Coccioni, 2008). In adverse conditions, foraminifera assemblage composition will be affected, and foraminifera employ biological defense mechanisms that can manifest as detectable morphological changes, providing tangible evidence of environmental stress (Ben-Eliahu et al., 2020). By analyzing living populations and surface-sediment assemblages, researchers can assess the current state of benthic ecosystems, whereas stratified sedimentary deposits offer insights into past environmental conditions (Bergamin et al., 2019).

In this context, the objective of this study is to combine sedimentological parameters and foraminifera as proxy bioindicators to evaluate changes in abundance and community structures in response to shifting environmental conditions during the pre- and post-industrialization periods (i.e. before and after 1900 Common Era [CE]) in the dynamic environment of the BSI. These results can be further used to reconstruct historical environmental conditions, assess the effect of anthropogenic activities over time, and provide insights for future environmental monitoring in similar coastal regions using forams as bioindicators.

The BSI covers an area of approximately 100 km². It has an average depth of 25 m, with some areas plunging as deep as 80 m, although more than one-third of the bay is shallower than 10 m (Dreujou et al., 2020; Figure 1). Despite its proximity to the sea, the bay exhibits relatively low salinity levels, hovering around 30 psu (Carrière, Le Hénaff, and Aubut Demers, 2018), mostly owing to the significant influx of fresh water from four major rivers that feed into the bay: Hall, des Rapides, aux Foins, and du Poste Rivers (Figure 1). Conversely, the saltwater component originates from the Gulf of Saint-Lawrence, which connects directly to the bay’s entrance, creating a dynamic and varied aquatic ecosystem. Notably, the BSI experiences substantial tidal fluctuations, with high tides surging up to 3.3 m, as observed in previous studies (Carrière, Le Hénaff, and Aubut Demers, 2018; Dreujou et al., 2020). This tidal regime contributes to the bay’s distinctive hydrodynamic characteristics. The bay’s water-circulation patterns follow a cyclonic path, and water residence time ranges from 2 to 12 days (Shaw et al., 2023).

Located on the eastern shore of BSI is the industrial-port city of Sept-Îles, home to approximately 25,000 residents in 2021 (Statistique Canada, 2024). The region is rich in history, with the Innu people being the first to settle there more than 8000 years ago, followed by sustained European presence and colonization from the 16th and 17th centuries on (La Société Historique du Golfe, 2024). In the mid-19th century, more settlers began arriving on the coast of Sept-Îles, attracted by its rich fishing grounds and forest resources. This growth led to the development of a fishing village where fishing activities became the primary livelihood; these changes marked the beginning of the area’s industrial era.

In 1910, the discovery of iron ore deposits marked a significant turning point on the coast that allowed for expansion beyond fishing and whaling. This discovery catalyzed the shift toward industrialization, transforming the area into a burgeoning industrial port. Between 1950 and 2000 CE, a significant influx of large industries and increased production occurred, leading to the establishment of Sept-Îles as an industrial-port city in the mid-20th century (Côté and Dubreuil, 2019). Today, its economy is largely reliant on heavy industries linked to mining, with four main companies dominating the region: the Alouette Smelter, the Iron Ore Company of Canada mining industry, the Société Ferroviaire et Portuaire de Pointe-Noire, and the Port of Sept-Îles. Of these, Aluminerie Alouette, situated on Pointe-Noire, stands out as the largest aluminum smelter in the Americas, boasting an annual production capacity exceeding 630,000 tons of aluminum (Aluminerie Alouette, 2024; Natural Resources Canada, 2023). Additionally, the Port of Sept-Îles, with its deep-water facilities, serves as a crucial ore handling port in Canada, handling 32.1 million tonnes and accommodating numerous vessels annually (Port de Sept-Îles, 2024).

The region surrounding the BSI is not only an industrial hub but also a center for various maritime activities; fishing and mariculture are prominent industries. Furthermore, the area is known for its rich marine biodiversity, which attracts enthusiasts and researchers interested in observing marine mammals and engaging in cruise tourism (Carrière and Le Hénaff, 2018). The port has a dedicated cruise dock that is capable of accommodating large cruise ships from around the world.

The Sept-Îles region, part of the Canadian Shield, is situated within the Grenville geological province, which is renowned for its iron (Fe) and ilmenite mines and its potential for industrial minerals (Namur et al., 2010). The BSI encompasses two main physiographic features: the coastal plain and the Sept-Îles Archipelago (Higgins, 2005). The coastal plain is characterized by marine and fluvial terraces, initially formed during the Goldthwait Sea invasion and later shaped by the Sainte-Marguerite and Moisie Rivers, bordered by the Laurentian Plateau to the north and the Gulf of St. Lawrence to the south. The archipelago includes islands predominantly comprising intrusive rocks.

The coastal region around Sept-Îles, covering roughly 10 to 15 km, is characterized by ancient igneous and metamorphic formations dating back to the Precambrian and Cambrian periods, linked with the Grenville province (Boyer-Villemaire et al., 2013). These formations were intruded around 540 million years ago, creating the Sept-Îles intrusive suite (SIIS). The Sept-Îles intrusive suite comprises three major units: Layered Series, Upper Border Series, and Upper Series (Figure 2). The Layered Series primarily comprises gabbro-troctolite, the Upper Border Series is dominated by anorthosite, and the Upper Series features granite-syenite layers. Iron–titanium (Fe–Ti) oxide layers, which are a notable feature of the Layered Series, are more dominant in the western part of the intrusion.

As these formations extend further offshore, they intersect with the monoclinal St. Lawrence platform, mainly comprising Ordovician limestone (Higgins, 2005). This geological intersection reveals a faulting system that developed post-Ordovician, with a SW-NE orientation. Over time, differential erosion has formed beach terraces along the coastline, with soft cliffs (∼30 m) of postglacial material and igneous rock foreshores around the islands. Offshore, postglacial sediment thickness ranges from 10 to 80 m (Hein et al., 1993). The average sedimentation rate around the Sept-Îles Archipelago is 0.14 cm/y (Lajeunesse et al., 2007). Currently, sandy sediments eroded from the coastline, especially from the Moisie River Delta, are transported offshore to deeper basins through canyons and channel-levee systems (Lajeunesse et al., 2007). Sedimentation rates within the Bay are unknown thus far, probably because of the shallow depth of the bay and sediment remobilization.

The methodology encompasses three main stages: sample collection, laboratory processing, and statistical analysis. Initially, samples are collected from designated sites and prepared for further examination in the laboratory, where specific tests and measurements are conducted. Following this, the data undergo statistical analysis to interpret patterns and relationships.

A 144 cm long sediment core (AMD-OSL2020-03GC; hereinafter referred to as core 03GC) was collected at the mouth of the BSI (50.17° N; −66.39° W; water depth, 80 m) using a gravity corer during the 2020 Odyssée St-Laurent winter expedition onboard the Canadian Coast Guard Ship (CCGS) icebreaker Amundsen. The coring site is located within a channel and was selected based on the available chirp profiles recorded on board the CCGS Amundsen that indicated sediment accumulation not influenced by mass wasting events (Lajeunesse et al., 2007). This location features a low slope (<2°; Lajeunesse et al., 2007). The sediment core was stored in darkness at 4°C until further processing.

In the laboratory, the sediment core was opened, described, and photographed with a GEOTEK multisensor core logger (MSCL) at Institut des Sciences de la Mer de Rimouski–Université du Québec à Rimouski. The diffuse spectral reflectance was measured using a Konica Minolta CM2600d spectrophotometer mounted on the GEOTEK MSCL at 1 cm intervals (St-Onge et al., 2007). The reflectance data are expressed as L*, a*, and b* according to the Commission Internationale de l’Eclairage. L* ranges from 0 (black) to 100 (white), a* ranges from +60 (red) to −60 (green), and b* ranges between +60 (yellow) and −60 (blue; St-Onge et al., 2007). For this study, sediment core 03GC was subsampled only evenly every 1 cm up to 30 cm. This interval includes the end of the pre-industrial period and the entire industrial period.

To investigate the variation in mineralogy and understand the provenance of present sedimentation in the BSI, a total of 10 surface samples was collected, with five from the western BSI (WBSI) and five from the northern and eastern parts of the BSI (NEBSI; Figure 1). These samples were collected using a grab sampler, and the top 2 cm of each sample was carefully scooped out from the top of the grab using a spatula.

The chronology of core 03GC was assessed by combining radiocarbon (¹⁴C) dating with lead (²¹⁰Pb) and cesium (¹³⁷Cs) measurements. For radiocarbon analysis, shell fragments were gathered at a depth of 25 cm. These fragments underwent cleaning with distilled water, followed by drying in an oven at 60°C for 2 hours and subsequent weighing. Next, shell fragments were pretreated at the Radiochronology Laboratory of the Centre d’études nordiques (CEN Centre for Northern Studies, Université Laval) and dated by accelerator mass spectrometry (AMS) ¹⁴C at the Keck-Carbon Cycle AMS facility (University of California, Irvine). For ²¹⁰Pb and ¹³⁷Cs measurements, sediment samples from the top 22 cm of core 03GC (at 1 cm intervals between 0–8 cm and at 2 cm intervals between 8–24 cm) were lyophilized, homogenized, and sealed in glass vials for >21 days to allow secular equilibrium. The measurements of radioisotopes ²¹⁰Pb and ¹³⁷Cs were performed at the Radiochronology Laboratory of the Centre d’Études Nordiques (Université Laval) by gamma spectrometry following Reyss et al. (1995). The ²¹⁰Pb excess activities were calculated by subtracting the ²²⁶Ra-supported activity from the total ²¹⁰Pb activity.

The age-depth model was constructed by combing ²¹⁰Pb and ¹⁴C data in the R package rplum version 0.4.0 (Aquino-López et al., 2018; Blaauw, Christen, and Aquino-López, 2023). Without premodeling the ²¹⁰Pb dates, this program allows the estimation of the best weighted mean age for each depth with a 95% confidence interval using a Bayesian approach (Aquino-López et al., 2018). The conventional ¹⁴C age was calibrated directly in rplum using the Marine20 calibration curve and a regional marine reservoir age correction (ΔR) of −39 ± 40 years, corresponding to the reservoir age correction reported by a modern shell sample from Île Niapiskau (Havre-Saint-Pierre, Québec) located near the study area (Heaton et al., 2020; McNeely, Kyke, and Southon, 2006). Finally, the ¹³⁷Cs profile was used to validate the age model by comparing it with the initial increase of ¹³⁷Cs activity that is assumed to correlate with the first fallout (FF) from atmospheric nuclear testing in the mid-1950s (Appleby, 2001).

Loss on ignition (LOI) was performed at 550°C to assess the amount of organic matter in each subsample using the method outlined in Heiri, Lotter, and Lemcke (2001). Initially, weighed amounts of wet sediment from each subsample were dried in a conventional oven at 105°C for 24 hours. They were then reweighed and placed in a muffle furnace at 550°C for 4 hours before being weighed again. To ensure that no moisture remained in the samples, they were returned to the conventional oven at 105°C for 24 hours and weighed one final time.

Sediment grain-size analyses (<2 mm fraction) were conducted using a Horiba (model LA950v2) laser diffraction grain-size analyzer. Before the analyses, approximately 0.2 g of ignited sediment samples were moistened with distilled H₂O and deflocculated by adding 30 mL of a Calgon solution (sodium hexametaphosphate; 20 g/L⁻¹) and disaggregated using an in-house rotator for 12 hours before measurement. The grain-size distribution (clay, silt, sand) and statistical parameters (e.g., D₅₀) were calculated using GRADISTAT software version 9.1 (Blott and Pye, 2001).

Eleven elements (Al, Si, P, K, Ca, Ti, Mn, Fe, V, Sr, Zr) underwent analysis through energy-dispersive X-ray fluorescence (ED-XRF) spectrometry using the Malvern Panalytical Epsilon 3-XL instrument and following the procedure used by Gamboa et al. (2017). Before ED‐XRF analysis, LOI was determined gravimetrically by heating the dried samples at up to 950°C for 3 hours. Subsequently, ∼0.6 g of ignited samples was treated by borate fusion (pure composition of 49.75% Li₂B₄O₇, 49.75% LiBO₂, and 0.5% LiBr) in an automated fusion furnace (CLAISSE^® M4 Fluxer) to form glass disks before being analyzed with the spectrometer. Following the analysis, acquired ED‐XRF spectra were processed with the standardless Omnian software package (PANalytical) (Malvern Panalytical, 2024). The resulting dataset was expressed in terms of weight percentage (wt.%) for major elements (Al, Si, P, K, Ca, Ti, Mn, Fe) and micrograms per gram (μg/g) for trace elements (V, Sr, Zr). Procedural blanks consistently represented less than 1% of the lowest concentration measured in the sediment samples. Analytical precision and accuracy were ascertained to be superior, with deviations falling within the range of 1%–5% for major elements and 5%–10% for trace elements, as checked against an international standards of United State Geological Survey SDC‐1 (U.S. Geological Survey, 2022a) and BCR-2 (U.S. Geological Survey, 2022b) and analysis of replicate samples.

For XRD analysis, around 1 g of each sediment sample was mixed with 0.25 g of corundum, which served as an internal standard. The samples were then finely ground in a McCrone micronizing mill using 5 ml of ethanol to achieve a uniform powder. This slurry was dried overnight in an oven at approximately 60°C then further homogenized with an agate mortar. To prevent fine particles from clumping, 0.5 ml of Vertrel was added. The dried and homogenized powder was sieved to less than 300 μm, loaded into sample holders, and analyzed using a Malvern Panalytical X’Pert Powder diffractometer.

The XRD scans were performed from 5° to 65° 2θ, with a step size of 0.02° 2θ and counting time of 4 seconds per step. The resulting diffraction data were processed using ROCKJOCK v11 software, which employs a full pattern-fitting method to quantify the major mineralogical components of the samples (Eberl, 2003). This method provides a precision of ±3 wt.%, and the total mineral wt.% was normalized to 100%. The analysis focused on quantifying 15 key minerals, including quartz, K-feldspar, plagioclase, calcite, dolomite, amphibole, Fe-bearing minerals, amorphous silica, kaolinite, chlorite, illite, biotite, muscovite, smectite, and vermiculite. These minerals accounted for more than 96% of the total mineral content in the bulk sediment samples.

For foraminifera analysis, 13 subsamples were selected from the top 22 cm of the core. Adhering closely to the protocols outlined by de Vernal, Henry, and Bilodeau (1999), for each sediment sample, 10 cm³ of wet sediment was measured, dried at ambient temperature, and reweighed. The dried samples were then carefully wet-sieved through 125 and 63 μm mesh sieves, with the resulting fractions subjected to detailed examination under a binocular microscope at magnifications of 40× and 80×.

Agglutinated and calcareous benthic foraminifera were sorted by hand with a paintbrush, with identification and enumeration conducted in accordance with the taxonomic classifications of Rodrigues (1980), Loeblich and Tappan (1987), and the World Foraminifera Database (Hayward et al., 2024). Planktonic foraminifera were counted to mark their emergence but excluded from the analysis (abundance and diversity calculation) to avoid diverting the focus from benthic foraminifera. Also, occurrences of species less than 10 times altogether were excluded from the numerical analyses.

Foraminifera specimens were examined and categorized based on their relative abundance, ranging from rare (<5%) to occasional (from 5.1%–10%) to common (>10%).

This study used the nonlinear unmixing Excel macro SedUnMix and the R package FingerPro to gain a quantitative understanding of the pre- and post-industrial changes in sediment provenance (Andrews and Eberl, 2012; Lizaga et al., 2020). Based on an environmental study of the BSI where the mineralogy of sediments is analyzed (Carrière, Le Hénaff, and Aubut Demers, 2018), it is suggested that sediments at the entrance of the bay derive from the mixing of two main source areas: NEBSI and WBSI. Linear discriminant analysis (LDA) was performed with these potential sediment sources to demonstrate the different mineralogical compositions of these sources. The LDA was conducted with the R software using the packages adegenet and compositions (Jombart, 2008; van den Boogaart and Tolosana-Delgado, 2008). Based on the LDA result, SedUnMix was run on the normalized (100%) data for the six main minerals (quartz, K-feldspar, plagioclase, amphibole, Fe-oxides, clays) that represented more than 91% of the overall mineral concentration in the sediment samples, although the FingerPro modeling was run on centered log-ratio (CLR) of both quantitative XRD (qXRD; same as SedUnMix) and all ED-XRF data. These two sediment-mixing modeling methods are complementary to each other because one uses only qXRD data and the other uses both qXRD and ED-XRF data.

Multivariate statistical analyses were performed on the elemental geochemical and foraminifera data using R software (R Core Team, 2022). A CLR transformation was applied to the data using the coda.base and compositions packages (Comas-Cufí, 2023; van den Boogaart and Tolosana-Delgado, 2008). The CLR transformation addressed the constant-sum constraint inherent in compositional data, making it suitable for classical statistical methods (Aitchison, 1990; Montero‐Serrano et al., 2010). The transformed data was visualized using the tidypaleoand ggplot2 packages, where depth profiles of the elements were plotted with depth displayed from top to bottom (Dunnington et al., 2022; Wickham, 2016). To handle unequal depth intervals, linear interpolation was applied to create evenly spaced depth intervals, ensuring consistency across the dataset for further analysis. This process helped standardize the depth intervals, which is crucial for stratigraphic analysis. For hierarchical clustering, the Constrained Incremental Sums of Squares method was applied using the tidypaleo package (Grimm, 1987). This clustering approach was used to detect stratigraphic zonation and identify significant compositional changes in the dataset. The results of the Constrained Incremental Sums of Squares clustering were visualized alongside the geochemical profiles, with stratigraphic zone boundaries clearly marked. The plots were further annotated with an age-depth model to represent the corresponding ages of the sediments, and the final visualizations were saved for reporting.

Spearman’s rank correlation was computed using R to assess relationships between species abundances and elemental geochemical data, excluding depth. The corrplot package was used to visualize the correlation matrix (Wei and Simko, 2023). Correlations with coefficients more than |0.7| were identified as significant and tabulated. Redundancy analysis (RDA) was performed with the elemental geochemical and foraminifera data using the R package vegan (Oksanen et al., 2020). The RDA offered a visual representation of the direct influence of environmental gradients on species distributions as well as a comprehensive visualization of the relationships between elemental geochemical and foraminifera data with samples.

The results provide detailed insights into the physical properties, mineralogical and geochemical composition, and foraminiferal assemblages of the core. These findings reveal patterns and changes in the sedimentary environment, depositional history, and ecological dynamics of the BSI region.

Total ²¹⁰Pb activity at the surface was around 63.5 Bq/kg, decreasing with depth, indicating consistent sediment deposition over time (Figure 3A). The age-depth model for core 03GC suggests an age of ∼1735 CE at 30 cm up to the year of coring (2020 CE) at the surface (Figure 3C). The mean sediment accumulation rates average ∼0.10 cm/y through the core (Figure 3D). Similar sediment accumulation rates have been obtained in a nearby core (12BC) by Lajeunesse et al. (2007). The mean modeled age at 6–7 cm (∼1955 CE) corresponds with a marked increase in ¹³⁷Cs activity (Figure 3B), which is assumed to correlate to the FF event from atmospheric nuclear testing that occurred in the mid-1950s (Appleby, 2001). Overall, the age-depth model presented here is as robust as possible given all the available radiogenic data. Based on this age-depth model, the boundary between the pre-industrial (pre-1900 CE) and industrial (post-1900 CE) periods is defined in core 03GC at around 12–13 cm.

The core exhibits considerable variation in its L* (luminosity), a* (color), and b*(hue) values (Figure 4), spanning from 20 to 50 for L*, −1 to 1 for A*, and −4 to 10 for B*. From depths of 30 cm to 25 cm, the sediments are darker, as indicated by low L* values and high LOI (550°C). From 25 cm to the surface, the sediments become progressively lighter, with further increases in L*, a* and b* values except at 12 cm where a fine layer of dark-colored sediments (high LOI value at 550°C) occurs. The LOI (at 550°C) values fluctuate between 1% to 3% from 20 cm upward.

Based on the Shepard’s classification of grain-size plot (Figure 5), the sediments from the analyzed depths predominantly comprise sand and silt with minor clay content (<0.7%). The grain-size variation and D₅₀ profile exhibit a certain variability with depth (Figure 5). The sediment composition stabilizes at a depth between 20 cm and 30 cm, with silt content around 60%–70% and relatively low D₅₀ values. At 20 cm, a notable peak in D₅₀ reaches 55 µm because of the highest sand content and lowest silt content, with no clay present. At 16 cm, the sediment is the finest, with the lowest D₅₀ value of 30 µm, marked by a decrease in sand, an increase in silt, and the peak of clay content. Moving upward, the D₅₀ reaches 60 µm at 15 cm, driven by high sand content, followed by fluctuations between 14 cm and 4 cm. At 2 cm, a decrease in sand and an increase in silt and clay occur, leading to a finer sediment size with a D₅₀ value of 30 µm. At the surface (0 cm), the sediment predominantly comprises sand and silt, resulting in the highest D₅₀ value of 60 µm, indicating a coarse grain size. Overall, the highest D₅₀ is at the surface, and the lowest is at 16 cm, with significant grain-size fluctuations between depth of 20 cm and 12 cm.

The stratigraphical distribution of ED-XRF data (Figure 6) reveals that Si (27%–29%), Al (7.5%–8%), Fe (4%–6%), and Ca (2.7%–3%) are the dominant elements found in the core. For the minor elements, Sr is the most abundant, followed by Zr, V, Zn, and Y. Two distinct zones are identified from 30 to 17 cm and characterized by low levels of Fe and Mn with increasing concentrations of Ca, Si, and Al from the bottom, peaking around 17 cm (Figure 6). Zn, K, Y, and Zr show a similar trend, with Mg reaching its maximum at around 19 cm before declining. At 18–17 cm, a transition zone is observed, where Ca decreases significantly and Fe and Mn begin to increase. Other elements also exhibit slight changes. From 16 to 0 cm, the sediment profile becomes more variable, with a continued decrease in Ca, Si, Al, Zn, Sr, and Y and an increase in Fe and Mn. This stratification highlights the varying elemental composition throughout the sediment core.

Table 1 compares the percentage values of various elements, showing their apparent effects threshold value (%) and average concentration (%) from BSI. The apparent effects threshold values are derived from the National Oceanic and Atmospheric Administration Screening Quick Reference Table cards (Buchman, 2015). As per the Table 1, apart from Al, all other elements are under threshold concentration.

The stratigraphical distributions of the mineralogical data reveal that plagioclase (35%–40%), quartz (15%–18%), K-feldspar (13%–17%), amphibole (10%–12%), and clays (12%–14%) collectively constitute the majority of the overall mineral concentration (Figure 7). Additionally, the clay mineral assemblage primarily comprises illite (54%), chlorite (41%), and kaolinite (6%), which highlights a predominant presence of illite and chlorite. Figure 7 shows a significant shift in mineral concentrations during the post-industrial period. Quartz, plagioclase, amphibole, and clay minerals tend to decrease toward the top, whereas K-feldspar, pyroxene, apatite, and Fe-oxides increase. Notably, Fe-oxide, amphibole, and K-feldspar exhibit their minimum values at 12 cm and peak at the surface. Based on these observations, two major clusters are identified: 0–10 cm and 12–30 cm.

The ternary plots for Fe-oxides, K-feldspar, and clay, as well as for Zr, Fe, and Zn, show that the two potential sediment sources from the BSI region have distinct mineralogical compositions, providing a clear basis for sediment differentiation (Figure 8). The WBSI source is characterized by a higher concentration of Fe-oxides and K-feldspar, with a moderate presence of clay. In contrast, the NEBSI source is dominated by K-feldspar, with moderate clay content and minimal Fe-oxides. Similarly, in the Zr-Fe-Zn ternary plot, the WBSI source is enriched in Fe, whereas the NEBSI source displays a more balanced composition of Zr, Fe, and Zn. These clear mineralogical differences between the sources validate their use in sediment provenance analysis, highlighting distinct sediment characteristics from each region.

Based on the ternary plot in Figure 8A, the mineralogical compositions of the sediment core and the potential sediment sources from the BSI region display slight variation. The core primarily comprises 35%–50% clay, 45%–55% K-feldspar, and 0%–10% Fe-oxides, indicating a mixture rich in clay and K-feldspar but with minimal Fe-oxides. In contrast, the source from the WBSI has higher concentration of Fe-oxides (20%–40%), 20%–45% clay, and 60%–80% K-feldspar. Meanwhile, the source from the NEBSI is characterized by a high proportion of K-feldspar (60%–70%), moderate clay content (20%–30%), and minimal Fe-oxides, indicating that NEBSI predominantly supplies K-feldspar-rich sediments with low Fe-oxide content.

Similarly, based on the ternary plot for Zr-Fe-Zn, the core sample contains 20–35% Fe, 35%–40% Zn, and 30%–45% Zr (Figure 8B). In contrast, the source from the WBSI is characterized by a higher Fe content, ranging from 40%–65%, with lower concentrations of Zr (10%–30%) and Zn (30%). Meanwhile, the source from the NEBSI has a more balanced composition with approximately 30% Fe, 40% Zr, and 30% Zn, indicating a relatively even distribution of these elements.

The LDA plot demonstrates that the sediment sources from WBSI and NEBSI have unique compositions, characterized by variations in minerals and elements such as K-feldspar, Fe-oxides, clays, Zr, Fe, and Zn (Figure 8C). These clear distinctions in mineralogical characteristics allow for a high degree of discrimination between the sources. These differentiated sediment sources were subsequently used in SedUnMix and FingerPro modeling to further analyze the sediment composition and provenance in the BSI region, which highlights that the primary source of sediments in the core is from NEBSI (Figure 8D). At the bottom of the core, a clear dominance of sediments from NEBSI is found, with contributions reaching 100%; however, as the core is moved up, minor fluctuations occur, with a noticeable contribution from WBSI around 1900 CE. From this point toward the surface, a gradual increase occurs in the contribution from WBSI, reaching up to 20% according to SedUnMix and up to 40% according to FingerPro. In both cases, a significant change occurs in sediment provenance within the core, which aligns with the observed changes in mineralogy in the post-industrial section of the core.

The foraminiferal analysis of the sediment profile reveals distinct patterns in species composition and diversity across various depths (Table 2). Figure 9 presents the species found throughout the core. Starting from the deepest layers (22-12 cm), Sphaeroidinella dehiscens is the dominant species (Figure 9A, B, C, D), with its highest relative abundance occurring at 14 cm (52%). The Shannon Diversity Index in these layers varies, with moderate values recorded, such as 1.7 at 19 cm to the lowest value of 1.2 at 22 cm. At 16 cm, Adercotryma glomerata (Figure 9E) becomes the dominant species, with a relative abundance of 21% and a Shannon Diversity Index of 1.7.

In the middle layers (10-4 cm), Spiroplectammina biformis (Figure 9O) becomes the most abundant species, particularly at 4 cm where its relative abundance peaks at 50.3%. However, the dominance of S. biformis in these layers is associated with moderate Shannon Diversity Index values, such as 1.8 at 8 cm and 1.3 at 6 cm.

In the top layers (2-0 cm), Miliammina manitobensis (Figure 9K) is the dominant species, with the highest relative abundance reaching approximately 30% at 2 cm. These layers exhibit a relatively high Shannon Diversity Index of 1.9, indicating a more diverse assemblage near the surface.

This study found distinct patterns in the distribution of foraminiferal species within the uppermost 22 cm, with notable changes in relative abundances from the bottom to the top (Figure 10). The core is dominated by agglutinated taxa, with only one calcareous species present: S. dehiscens, which are found exclusively in the lower section, below 12 cm. Rosalina bradyi and Eggerella medius are prominent at the bottom, disappearing around 16 and 20 cm, respectively, before reappearing at the surface. A. glomerata shows a decreasing trend toward the surface, while Reophax nana is present throughout the core but increases steadily with time. E. medius reappears at 6 cm depth. The upper 12 cm of the core is dominated by taxa Trochamminopsis pusilla, Cribostomoides sp., R. nana, M. manitobensis, and S. biformis. A. cassis is present at the bottom, then disappears between 12 and 6 cm depth, but reappears at 4 cm and is the dominant taxon in the core top. Haplophragmoides sp. is restricted to depths below 8 cm, while Reophax subfusiformis is found only in the top 6 cm.

Based on the Spearman correlation analysis, several significant relationships between species and elements were identified. Species with highly correlated elements (>0.7 or < −0.7) are summarized in Table 3.

S. biformis exhibited multiple significant correlations. It showed a strong negative correlation with silicon (Si) (−0.82), Conversely, it displayed a strong positive correlation with iron (Fe) (0.84), suggesting a preference for iron-rich environments. These correlations highlight the specific elemental sensitivities and preferences of S. biformis.

Similarly, M. manitobensis displayed a strong negative correlation with Al (−0.80), suggesting it avoids environments with high aluminum content. S. dehiscens demonstrated a negative correlation with Fe (−0.80) but a positive correlation with Si (0.75), highlighting its elemental sensitivities.

A. cassis exhibited a strong negative correlation with Mg (−0.80), suggesting it is less prevalent in magnesium-rich environments, but a strong positive correlation with Zr (0.85), indicating a preference for zirconium-rich conditions. Lepidodeuteroammina sp. showed a negative correlation with Si (−0.73), further indicating species-element interactions.

Rosalina bradyi showed a strong positive correlation with Zn (0.76), suggesting it prefers zinc-rich environments, and T. pusilla displayed a positive correlation with Mn (0.78), highlighting its adaptability to manganese-rich conditions.

In addition, there were correlations between species themselves, such as a strong positive correlation between Haplophragmoides sp. and Cribostomoides sp. (0.80), further emphasizing ecological or elemental interdependencies.

The two first RDA axes (RDA1 and RDA2; Figure 11) explain 67.9% and 13,9%, respectively, of the variance in the data and the grouping of elements and species at different depths. In the first quadrant, Mg is the dominant element, strongly correlated with S. biformis at depths of 4, 6, and 10 cm. The second quadrant features K, Ca, Al, and Si, with S. dehiscens showing strong positive correlations with these elements at depths of 16, 20, and 22 cm. A. glomerata exhibits moderate positive correlations with the same elements, while Ammotium diversum and Rosalina bradyi show weaker associations. All species in this quadrant are negatively correlated with Fe and Mn. In the third quadrant, A. cassis and E. medius display weak correlations with Sr, Zn, Y, and Ti at 14 cm depth. In the fourth quadrant, Fe and Mn are dominant, with M. manitobensis showing strong positive correlations at depths of 0, 8, and 12 cm, while Haplophragmoides sp. and Cribrostomoides sp. exhibit moderate correlations. At the center, Ammotium diversum, R. nana, Lepidodeuterammina sp., and Trochamminopsis pusilla show weak correlations with the surrounding elements, indicating less distinct elemental preferences.

The multiproxy analysis of a short sediment core from the BSI covering the last 200 years of accumulation revealed notable environmental changes, likely linked to increasing anthropogenic impacts in the Bay starting in the late 19th century and intensifying throughout the 20th century.

Although the grain-size analysis suggests relatively stable deposition energy over time, with a dominance of sandy silt throughout the core, subtle shifts in sediment composition point to changing conditions. The deeper sections of the core, corresponding with the 19th century, were characterized by a prevalence of silt and low D₅₀ values, suggesting a calm depositional environment with minimal external disturbances that is likely driven by natural sedimentation processes. Around the mid-depth sections, corresponding with the late 19th century, changes in sediment composition begin to reflect the growing influence of human activities in the region. The increase in sand content and slight increase in D₅₀ values point to a period of higher energy conditions, potentially linked to development in the catchment, including the construction of a fishing port and a whale oil factory in the early 20th century in the BSI. Increased deforestation and construction of infrastructure likely increased erosion and could have disrupted the natural sediment supply, contributing to the introduction of coarser materials into the bay.

From the second half of the 19th century, a subtle change in elemental composition of the sediment was observed, with sustained increases in Fe and Mn and decreases in Si, Al, K, Ca, Zn, and Sr into the present. These changes appear to be the first indication of a trend of pre-industrial environmental change in the BSI. They were followed in the 1920s by synchronous variations in the bulk mineral and foraminifera composition of the sediment at depth of 11 cm, which indicate a further change in the environment with a shift in provenance of the sediments from eastern to western sources accompanied by the total disappearance of calcareous foraminifera and a hiatus in the presence of certain other taxa. This likely indicates that limited calcium availability became prevalent in the ecosystem, prompting a change in the foram assemblage. This shift in species dominance points to changes in habitat conditions and ecological dynamics between the pre-industrial and industrial periods. The decrease in Ca concentrations from the late 19th century to the present may be attributed to the dissolution of calcium-bearing minerals, although the exact cause remains uncertain. The ecological responses of benthic foraminifera to these changing environmental conditions highlights the resilience and adaptability of certain taxa in the face of anthropogenic influences, notably certain agglutinated species such as R. nana, which appears not to have been affected by these changes over time.

The most significant changes in sedimentology are observed in the upper section of the core (mid-20th century), with an increasing contribution of sediments from the western part of the bay (WBSI). The contribution reaches up to 40% at the surface and signals the sustained growing impact of industrial activities, particularly construction of infrastructure and development of port operations, which continued to grow over the course of the second half of the 20th century and to this day. These activities have likely resuspended sediments and introduced coarser materials into the depositional environment, transforming the natural sedimentation processes within the bay.

The second most significant change in the foraminiferal assemblage occurred in the recent past (early 21st century) at a depth of 3 cm, suggesting that environmental changes are ongoing in the BSI. However, this change was not recorded in the core’s sedimentology, which leads to speculation that it might perhaps have been a new type of perturbation that affected the forams at that time. In fact, in 2013, an unprecedented environmental catastrophe occurred in the bay, whereby the Cliffs Natural Resources company was responsible for an oil spill totaling many hundreds of thousands of liters of heavy fuel oil. Foram assemblages are known to be affected by this type of contamination, which can also lead to a decrease in their numbers (Joshi, Arya, and Saulnier-Talbot, 2024).

The results of our study show that even though the BSI is an ideal setting for industrial activities—including a short water residence time of only a few days, which helps to evacuate any contamination—the ecosystems is not immune to the pressures brought about by industrialization in the region. This study found that elemental and bulk mineral composition, as well as foraminiferal assemblages, were affected by industrial and urban development over the course of the last 150 years. Foraminifera could be used in regular biomonitoring in the region because they react quickly to environmental disturbance.

This analysis reveals that the BSI remains a relatively stable system, with metal concentrations generally below contamination thresholds; however, the observed increase in the relative abundance of certain metals, notably Fe and Mn, warrants attention because it may signify environmental changes. The subtle shift in metal concentrations during the industrial period has notably affected the foraminifera community diversity, underscoring their sensitivity to environmental fluctuations and their potential as bioindicators of anthropogenic stress. The foraminifera analysis vividly portrays the minimal changes in the geochemistry of the bay, with distinct distributions of agglutinated and calcareous foraminifera reflecting varying habitat conditions over time. The disappearance of calcareous species from the first half of the 20th century, accompanied by reduced calcium concentrations, indicates a stressful environment for these types of forams in the BSI’s current conditions, likely influenced by factors such as calcium limitation and carbonate dissolution. Moreover, the prevalence of agglutinated species in the upper layers suggests their greater tolerance to these environmental changes compared to calcareous species. Among these, Spiroplectammina and Reophax emerge as particularly resilient species, potentially serving as key indicators of environmental stress in the region. The shift in sediment provenance from NEBSI to WBSI in the uppermost layers, driven by industrial activities such as construction and mineral port activities, further emphasizes the role of anthropogenic factors in shaping the bay’s sedimentary and ecological dynamics. These findings underscore the vital role of foraminifera as sensitive indicators of environmental change, providing crucial insights into the historical and current conditions of the BSI.

The authors acknowledge the captain, officers, crew, and scientific participants of AMD20-OSL onboard CCGS Amundsen for the recovery of the sediment core analyzed in this study. This research was supported by the Université Laval EcoZone research chair in partnership with the Port of Sept-Îles and INREST; by the RQM through the Odyssée Saint-Laurent research program; by the Natural Sciences and Engineering Research Council of Canada (NSERC) through Discovery Grants provided to É. Saulnier-Talbot and J.-C. Montero-Serrano; and by Québec-Océan, Geotop, and Amundsen Science. NJ extends her gratitude to Dr. Vincent Bouchet for his invaluable guidance on taxonomy and thanks the ISF 2023 organizers for significantly enhancing her knowledge about foraminifera. The authors also express sincere appreciation to the Joseph A. Cushman Foundation, Dr. Christopher Makowski, and the ICS2024 organizers for their support in presenting this research at the symposium. An executive summary of this research was originally published in Joshi et al. (2024).