Spatial and Vertical Benthic Foraminifera Diversity in the Oxygen Minimum Zone of Mejillones Bay, Northern Chile

Introduction

Coastal areas are the most dynamic environments on the planet due to the interaction between the continent and the sea, which allows for a great diversity of high-productivity habitats (McLean and Tsyban, 2001). The benthos is an ecosystem that provides a large number of services, performing important functions in the oceanic element cycle (nutrients and micronutrients) which are necessary for biogeochemical processes (Giblin et al., 1995; Snelgrove, 1997). It also plays a fundamental role in removing elements, including pollutants, and houses a production of macrofauna of commercial interest for human consumption (Snelgrove, 2013). One of the most important environmental variables in this type of ecosystem is oxygen, which directly affects nutrient cycles (Murray, 2006).

In areas with intense upwelling, excess oxygen consumption generates anoxia. When this situation is permanent, it can form an Oxygen Minimum Zone (OMZ) (oxygen levels < 1 to > 3 ml l^–1) (Sen-Gupta and Machain-Castillo, 1993); a low current circulation and ventilation is also required in order to inhibit oxygen exchange (Sarmiento et al., 1988; Levin, 2003; Davis et al., 2021). Where the levels can be undetectable, generating a high productivity and respiration. The expansion of these extreme ecosystems affects the biochemical processes and may cause the intensification of OMZs, being of critical importance to ocean health (Stramma et al., 2012; Levin, 2017).

In the southeast Pacific on the north coast of Chile (23°S), The Mejillones Bay is located in an OMZ, where generates low oxygen concentrations (less than 0.5 ml l^–1) extending from 50 to 200 m depth (upper limit of the euphotic layer) (Rodríguez et al., 1986, 1991; Marín et al., 1993; Morales et al., 1996; Escribano, 1998; Marín and Olivares, 1999; Escribano et al., 2000; Páez et al., 2001). This area is characterized by presenting one of the most productive upwelling centers in the area known as Angamos Point (Páez et al., 2001; Valdés, 2004; Valdés et al., 2008), generated by the presence of winds toward the Equator during most of the year (Letelier et al., 2012). This area is adjacent to the Atacama Desert, being one of the most arid regions in the world, so the incidence of continental contributions is limited only to the atmospheric transport of lithogenic particles (Vargas et al., 2004). For this reason, coastal upwelling strongly influences the chemical composition of the surface waters of the Bay (Valdés et al., 2008). In addition, it is considered an area of high diversity (hotspot) and biological productivity due to its thermal stability (Marín and Olivares, 1999; Marín et al., 2003), yielding annual values of 1,070 g C/m² (Marín et al., 1993). This area is scarcely influenced by currents, generating a sedimentary deposit (Ortlieb et al., 2001; Valdés et al., 2005), the extreme conditions of which favor the presence of species capable of tolerating low oxygen concentrations. Despite being an area with a great diversity of species, it is highly intervened from the twentieth century. Due to the development of industrial activities such as aquaculture, fisheries, port terminals for the export of minerals and thermoelectric plants (Valdés et al., 2005), which may be affecting the environmental characteristics of the bay and the distribution of marine organisms.

Therefore, this OMZ, when grouping a set of particular ecological characteristics (e.g., hotspot, high productivity, anthropic intervention), it is necessary to apply environmental and oxygenation proxies, which allow to understand and identify the changes of these extreme ecosystems over time and space. One of the proxies that can be used as a tool is foraminifera, one of the most abundant benthic fauna groups (Gooday et al., 2012). Their distribution extends over wide bathymetric ranges (5–5,000 m depth) (Gooday, 1986; Gooday et al., 2000; Bernhard et al., 2008) and they play a fundamental role in the biogeochemical cycles of organic and inorganic compounds (Haynes, 1981; Yanko et al., 1998). Foraminifera generally use oxygen from the sediment-water interface, where large amounts of organic matter accumulate depleting oxygen in the sediment. Endobenthic species go even deeper, entering several centimeters into the seafloor in search of substrates rich in organic matter; however, oxygen there is even more scarce compared to surface sediment (Sen-Gupta and Machain-Castillo, 1993). Optimal oxygen conditions for benthic foraminifera can thus be seen to vary widely, depending on the tolerance limits that allow their distribution even under extreme anoxia conditions. The drastic fluctuations in foraminifera communities generated by the variability of abiotic conditions makes them a good study model to develop palaeoecological reconstructions.

For this reason, an ecological reconstruction of the Mejillones Bay was carried out from benthic foraminifera. Since the characterization of this association of microorganisms allowed to understand the dynamics of this OMZ, detecting changes over time and space (e.g., intensification or expansion of the OMZs), being the first work carried out for this area. This facilitated the comparison of the distribution and diversity of benthic foraminifera with the studies previously conducted for the area (i.e., Páez et al., 2001) and other OMZs. As well as allowing to identify the main differences with respect to other OMZs, expanding the knowledge of these adverse ecosystems and how they can influence in the distribution patterns of the species.

Materials and Methods

Sample Collection

Data sampling was taken from Tavera-Martínez (2019). Collection was carried out in Mejillones Bay between 22°58′32.55″ S–70°19′18.72″ W and 23°01′32.97″ S–70°29′50.79″ W. In order to know the vertical distribution of benthic foraminifera, two sediment cores were collected using a polycarbonate gravity core which dimensions were 68 × 71 mm × 60 cm long: the first core was collected at 75 m depth (ZA), and the second at 100 m depth (Z1A) (Figure 1). The first eight centimeters were analyzed only, for a total of 14 samples. On the other hand, to evaluate the spatial distribution, 18 samples were obtained at 10, 50, 70–75, 90, and 100–110 m depth using a mini box core of 225 cm². The first 0.5 centimeters of these samples were analyzed. Then all the samples went through a set of sieves of 8″ diameter with stainless steel mesh of 63, 150, and 212 μm, and were dried at ambient temperature. Only the fraction of 212 μm was considered in this study. Later, the samples were subdivided into equal parts with a splitter according to their size, in order to count and identify the specimens in a more efficient way. It is important to highlight that the total community of foraminifera was analyzed (living + dead specimens). So, their succession times are different from each other and the dead foraminifera, they would be representing organisms that were alive a few years ago, hundreds or perhaps thousands of years ago (Murray, 2000). Which may be a limitation of this study. However, we have the dating and the recent layer is very thin and irrelevant in this study, so the proportion of living organisms is very low.

Sediment Geochronology

The Z1A core was dated according to the alpha decay of ²¹⁰Po in secular equilibrium with ²¹⁰Pb and using the constant rate of supply (CRS) model (McCaffrey and Thompson, 1980; Turekian et al., 1980; Table 1). This model assumes a constant flow of ²¹⁰Pb to the sediments resulting in unsupported or excess activities, that is, those above the ²¹⁰Pb activities produced by the decay of ²²⁶Ra within the sediments. The chemical procedure to determine radionuclides involved dry sediment digestion (∼ 0.5 g) with concentrated acids (HCl, HNO₃, HF) and the addition of ²⁰⁹Po as an internal tracer. The activities were quantified in a Canberra Quad Alpha spectrometer until an adequate statistical count was achieved. The ages in different sediment sections were determined according to the accumulated inventories of ²¹⁰Pb_xs (unsupported), established from their exponential decay with depth. The supported activity was estimated from the best exponential model established for the downcore activities distribution. Standard deviations of unsupported ²¹⁰Pb and inventories were quantified by propagating the counting uncertainties (Bevington and Robinson, 1992), and the age errors were estimated by first-order analysis (Binford, 1990). The CRS method was selected for the age estimations because it considers the bulk density variations in the sediment composition.

Taxonomic Diversity

The organisms were identified at the species level using a stereoscope and light microscope at the Laboratory of Foraminiferology of the University of Concepción, according to the guides and taxonomic keys of Loeblich and Tappan (1988), Acosta-Herrera (2004), and Holbourn et al. (2013), and research on the Chilean coast by Boltovskoy (1954), Pujol-Rodríguez (1977), Zapata et al. (1995), and Figueroa et al. (2005, 2006). Subsequently, photographic records of each of the species were obtained with an SU 3500 Hitachi electron microscope. The rarest species were stored in foraminifera holder plates.

Morphogroups

Different morphogroups were identified in each layer and station according to the morphology of the shells, as this reflects their different feeding strategies and habitats. The grouping was carried out according to descriptions by Scott et al. (2004): (a) Planispiral (coiling of the shell in a single plane); (b) Trochospiral (evolute chambers on one side of the shell and involute on the opposite side); (c) Uniserial (arrangement of the chambers in a single column); (d) Biserial (arrangement of the chambers in two linear columns); (e) Triserial (arrangement of chambers in three columns, observing three chambers per whorl), and (f) Millioline or Quinqueloculina (coiling of the chambers at 72^°, observing four chambers on one side and three on the other). Thus, changes in the relative abundance of each morphogroup reflected changes in the environment (Jones and Charnock, 1985; Corliss and Chen, 1988; Murray et al., 2011).

Faunal Association

Relative abundance was estimated for each station and/or layer, subdividing the sample with a cutter or splitter, in order to count the specimens more efficiently. Subsequently, the number of individuals of each species found was counted and the relative abundance was calculated in units of individuals/cm³. Species contributing a total abundance greater than 70% corresponded to the dominant species for each sample and/or layer. Finally, species diversity and richness were calculated using the Shannon Wiener index (H‘) as a measure of heterogeneity and the Simpson index (D) as a measure of dominance.

Shannon Wiener index:

Simpson index:

Where, Ni is the number of individuals of the same species and N is equal to total number of individuals.

Multivariate Analyses

The data set used to these analyses were adjusted with square root (Supplementary Annexes E,G). For the vertical distribution, a total of 32 species were analyzed in nine samples or sediment layers, while for the spatial distribution the number of species was 39 and 16 samples. A similarity cluster was carried out based on the Bray-Curtis index, which compares the abundances between species and builds a distance matrix. Subsequently, similarity percentages (SIMPER) were obtained identifying the species that contributed most to the formation of the groups in the cluster. In order to spatially represent and corroborate the groups formed by the cluster, a multidimensional scaling (MDS) analysis was carried out, considering the value of the stress coefficient that measures the degree of success of the graph (<0.05: excellent, < 0.1: good, < 0.2 potentially useful) (Carnahan, 2005). Finally, in order to statistically support the results obtained from the MDS, an analysis of similarities (ANOSIM) was developed, determining whether there were significant differences between the stations regarding the composition of species with a value of p > 0.05.

Results

Geochronology

The activities showed quite good exponential decay that allowed estimating unsupported activities (2.04 ± 0.71 dpm g^–1; r² = 0.95, p < 0.01) and a mean sedimentation rate of 0.11 ± 0.02 cm year^–1 was determined. 210Pbxs (in excess or unsupported) were found in the first 6 cm, taking previous₂₂₆Ra measurements and ²¹⁰Pb inventories into account (3.4 dpm g^–1—50 dpm cm^–2; Muñoz et al., 2004), the activities below 6 cm would be in the range of the supported activities. From this depth, the activities diminished drastically, lower than the predicted values, and were in the limit of the unsupported activities that probably correspond to the tsunami event of 1877 that will explain these lower activities from the 6–7 cm to 13–14 cm, corresponding to one deposition event (Figure 2). Unfortunately, this event occurred in the limit of ²¹⁰Pb dating, and it was not possible to establish the ages under this layer. The errors associated with the age estimations increased with depth, inherently to the CRS dating method (Table 1).

Spatial Distribution

Taxonomic Diversity

A total of 39 species distributed in 18 families and six orders were recorded. The order Rotaliida (Figure 3) presented the highest number of families and species, with 13 and 23, respectively. Orders Astrorhizida, Miliolida, and Textulariia were represented by only 1 species each. The genus Bolivina (Figures 3F–H) presented the largest number of species (seven), followed by Valvulineria (Figures 3C,D) presenting four species (Supplementary Annex A).

On the spatial distribution, Bolivinitidae and Buliminellidae were the most abundant families both with 30%. Consecutively Epistominidae was represented in 15% of the total abundance, followed by Nonionidae and Trochaminidae with 7 and 6%, respectively (Figure 4A). The dominant species in the 10 m isobath corresponded to B. elegantissima (Figure 5A), Bolivina seminuda, Rosalina cora, Rotaliammina squamiformis, Epistominella exigua, Buccella peruviana, Cassidulina crassa, Nonionella stella, Nonionella auris, Haynesina depressula, Quinqueloculina seminula, and Bolivina doniezi (Figure 6). These species contributed between 79 and 92% of the total abundance in each station. H. depressula registered as the dominant species at station T8 (Figure 6). In the 50 m isobath, the dominant species were B. seminuda, Buliminella elegantissima, E. exigua, B. peruviana, C. crassa, and N. stella (Figure 5C) and contributed between 81 and 95% of the total abundance. In the 70 m isobath, the dominance of most of the aforementioned species decreased, showing abundances between 21 and 35% for B. elegantissima, 24% for B. seminuda, and 26 and 28% for N. stella. Finally, in the 90 m isobath Trochammina sp. presented the highest abundance (27%), while B. seminuda represented 30 and 73% of the total abundance in the 100 m and 110 m isobaths, respectively (Figure 6 and Supplementary Annex C).

Morphogroups and Faunal Association

The trochospiral shape was the most predominant (54%), registering a greater percentage of these shells between 10 and 90 m depth. Biserial shells represented 26% of all registered foraminifera, being the most predominant in the isobates of 100 m and 110 m. Planispiral shapes were represented in 15% of abundance, being in greater proportion in the station T8-10 m. While spherical and milliolin shells were the least abundant (3 and 2%, respectively) (Figure 7A), its abundance was not dominant in any of the isobates.

According to the Shannon-Wiener diversity index (H′), the values of the 10 m isobath ranged between 1.07 (T8) and 2.08 (T7), while the 50 m isobath showed a minimum diversity of 1.42 (T7) and a maximum of 2.60 (T8). The 70–75 m isobath yielded the highest diversity at station T8 (2.09) and the lowest value was observed at station T7 (1.78). Similarly, the 100 m isobath yielded the highest value (2.09) compared to depths of 90 and 110 m depth, in which diversities of 1.85 and 1.25, respectively, were estimated (Figure 8A). In contrast, the higher Simpson (D) values reflected the dominance of H. depressula in the 10 m isobath (T8 station) and of B. seminuda at 110 m depth (BMS station). The general trend of spatial diversity of species indicated higher values of H´ and low dominances between 75 and 100 m depth. In contrast, most stations located in the isobath of 50 m, with the exception of T8, showed lower values of H´ and higher dominances.

Multivariate Analyses

Four clusters with the greatest similarity between stations were formed based on the Bray Curtis index values (Figure 9A). The cluster B showed a similarity of 56% between stations of 70–75 m (T6 and T7) and 90 m y 100 m (BMS), where the greatest similarity was found between samples of T6-70 m and T7-75 m (85%). The cluster A registered a 36% of similarity. This was formed by three stations of 10 m isobath (T7, T6, and BMS, 2010), that shared 55% of similarity with each other and two stations of 50 m isobath (T8 and T8 FDI), that showed 65% of similarity. The cluster C reported a similarity of 40% between stations of 75 m (T8, T4, and T1), 50 m (T4 and T7), 110 m (BMS, 2010) and 10 m isobaths (T4), where the highest similarity was found between T8-75 m and T4-75 m (70%). Finally, the station T8-10 m was the most dissimilar compared with the other stations, sharing a similarity of only 20% (Figure 9A). The SIMPER analysis indicated that the species which contributed the most to the formation of the cluster B were Hemisphaerammina sp. (90.84%) and B. patagonica (89.35%). Cibicides aknerianus and Bolivina doniezi contributed to 90.83% and 89.53 of cluster A. In addition to V. inflata that contributed to 90.07 and 90.64% of the subgroups formed in cluster C (50 m and 70–75 m; 10 m and 70–75 m, respectively) (Supplementary Annex H).

However, the MDS diagram indicated that there is no apparent group differentiation, with the exception of T8-10 m (Figure 9B), as corroborated by ANOSIM (R = 0.20 and p = 0.03), evidencing that there are no significant differences between and within the groups.

Vertical Distribution

Taxonomic Diversity

A total of 31 species were reported, distributed in 15 families and 5 orders. Most of the taxa found are part of the order Rotaliida (25 species and 10 families) (Figure 3). Conversely, Lituolida, Miliolida, and Robertinida were present in one family and one genus each. Bolivina and Cassidulina (Figure 3E) presented the greatest number of species with five and four, respectively, followed by Nonionella (Figure 3B) with two species (Supplementary Annex B).

Bolivinitidae was the most abundant family, presenting the greatest dominance proportion in core ZA (64%) (Figure 4B); meanwhile, this family along with Nonionidae were the most abundant in core Z1A with 29 and 30%, respectively. Buliminellidae was the second family with the highest number of individuals, registering 26% (Z1A) and 14% (ZA) (Figure 4C). Bolivina seminuda presented the highest abundance within the first centimeter of both cores (64% in Z1A and 77% in ZA) (Figures 10A,B). However, and despite its abundance in the subsequent strata, this species dominance decreased registering minimal abundances of 20% (2–3 cm layer, Z1A) and 26% (4–5 cm layer, ZA). Fursenkoina sp. (Figure 5D) and Buliminella elegantissima revealed the opposite, with increased abundances in the lower strata: Fursenkoina sp. registered maximum percentages of 28% (4–5 cm layer, Z1A) and B. elegantissima of 18% (1–2 cm layer, ZA). Similarly, the remaining species tended toward increasing abundances in lower substrates (e.g., N. auris and Nonionoides grateloupi). Finally, E. exigua abundance showed scarce variation throughout the core Z1A strata (Figure 10A and Supplementary Annex D).

Morphogroups and Faunal Association

In general, the trochospiral shells were the most abundant (36% Z1A and 38% ZA) and uniserial/biserial shapes represented 31 and 32% of the species. Meanwhile, biserial shells represented 26% on the spatial distribution and planispiral shells were observed between 24% (ZA) and 26% (Z1A) of all specimens. In contrast, spherical and millioline shapes were evidenced only in low proportions (3 and 4%) (Figures 7B,C). The vertical distribution of these morphogroups was evidenced with a dominance of biserial species from the most superficial layer of the sediment to 2–3 cm. On the other hand, trochospiral and planispiral species predominated in the deeper layers (3–8 cm).

The Shannon-Wiener diversity index (H′) values were high between 4 and 8 cm depth (2.01 and 2.16), while the Simpson dominance (D) was between 0.14 and 0.18. Conversely, the first centimeter of both cores registered the lowest diversity values (1.01 in ZA and 0.95 in Z1A) and the highest dominance values (0.60 in ZA and 0.47 in Z1A) (Figure 8B). The highest richness corresponded to 20 species in the 4–5 cm layer of core Z1A and 19 species in the 1–2 cm layer of core ZA.

Multivariate Analyses

The cluster analysis reflected the formation of three groups in the cluster, two of which conformed the first 3 cm of cores ZA and Z1A with similarities of 62% (cluster A and B). The cluster A was formed by a subgroup of 1–2 cm layer and 2–3 cm of core ZA, registering 75% of similarity, and the 0–1 cm layer of ZA, which shared a 62% of similarity with this subgropup. The cluster B was constituted by 1–2 and 2–3 cm layers of Z1A, showing 85% of similarity. Finally, the cluster C conformed the last 4 cm of core Z1A and the 4–5 cm layer of core ZA, with a similarity of 64% (Figure 11A). The SIMPER analysis indicated that the species which contributed the most to the formation of groups were B. patagonica (91.04%), Cassidulina sp. (91.22%) and V. fragilis (89.11%) to the Cluster A. N. stella contributed to 91.1% of the Cluster B. Followed by Suggrunda porosa and Bolivina variabilis with 88.91% and 86.39, respectively. B. elegantisima contributed to 91% to the formation of cluster C (91.74%). Similary to E. exigua and Pseudononion japonicus, each contributing 89% of cluster C (Supplementary Annex F). Like cluster analysis, MDS graph showed no differentiation between layers (Figure 11B), and These samples were very scattered among themselves. These results were confirmed by ANOSIM (R = -0.01 and p = 0.41), indicating no significant differences between and within the groups.

Discussion

Spatial Distribution

Taxonomic Diversity

The most abundant and diverse order corresponded to Rotaliida (79%), their high diversity of trochoid forms gives them the advantageous ability to adhere to various types of substrates, inhabiting most marine environments and especially coastal areas. They can therefore be found in neritic areas, both on the coastline and at the end of the photic zone (Haynes, 1981; Boix-Martínez, 2007).

The high abundance of Bolivinitidae encompassed eight species (Figure 7A), with Bolivina seminuda being the most abundant. The presence of these specimens reflects low oxygen conditions, coinciding with descriptions for other OMZs such as the Gulf of California, Arabian Sea (Sen-Gupta and Machain-Castillo, 1993) and central west coast of India (Mazumder and Nigam, 2021). Although the study carried out by Páez et al. (2001) for the Mejillones bay reported some species of the genus Bolivina (B. seminuda and B. costata), the present study recorded five new species of this genus not registered for this area before (e.g., Bolivina ordinaria, Bolivina plicata, Bolivina striatula, Bolivina earlandi, and Bolivina doniezi). So, these species represent the first records for the Mejillones bay, which could infer a change in the composition of species over time. It also indicates a high productivity environment due to seasonal upwelling and Equatorial Subsurface waters (Figueroa et al., 2005), as well as a high accumulation of organic matter.

The high abundance of Buliminellidae was represented only by Buliminella elegantisima (Figure 5A), with 13,976 individuals. This predominance indicates dysoxic environments (0,07 ml O₂ l^–1) (Murray, 2006) and an increase in organic matter pulses as a result of marine pollution (Sen-Gupta and Machain-Castillo, 1993). Epistominidae was represented solely by Epistominella exigua and indicates low oxygen conditions, presence of detritus deposits and eutrophied environments (Tyszka, 1994; Sagasti and Ballent, 2002; Smolén, 2012). These results coincide with the study carried out previously for the bay (Páez et al., 2001), which indicated that these species are abundant and typical for this zone. Nevertheless, B. elegantissima, E. exigua and genus Buliminella are not registered for the OMZ of Concepción in central-southern of Chile (Tapia et al., 2008). Therefore, it is inferred that the taxonomic diversity of benthic foraminifera differs from other OMZs, including those located within Chile.

Due to the dominance and high abundances of opportunistic infaunal and detritivore species (e.g., B. elegantissima and B. seminuda) registered in all isobaths (10, 50, 70–75, and 90–110 m), we can infer that the Mejillones Bay presents dysoxia and microxia conditions (except for the 50 m isobath) with high productivity and organic matter concentrations, indicating seasonal upwelling. Specifically at 10 m depth, the dominance of Haynesina depressula may indicate an increase in the salinity of the sediment (35 and 38‰), as reported in previous studies (e.g., Vanicek et al., 2000). On the other hand, the dominance of Fursenkoina sp. may represent a high sulfur content and organic nitrogen between 70 and 75 m depth, as reported by Sellier de Civrieux and Bonilla-Ruiz (1971), who mentioned concentrations of 0.520‰ and 0.287%, respectively. The proportions found of Nonionella stella would be given by the nitrate content in the water, from which they perform denitrification processes in their cell (Grzymski et al., 2002; Bernhard et al., 2012). Additionally, the predominance of detritivore and bacterivore foraminifera (e.g., Trochammina sp. and Bolivina plicata) in the 90–110 m isobath indicate hypoxic conditions and an increase in phytodetritus (Reolid et al., 2014; Table 2).

Morphogroups and Faunal Association

The trochospiral, uniserial, and biserial shells are infaunal species. They are grouped within an angular-asymmetric morphology (Khare et al., 2017). These types of shells indicate eutrophied silty substrates, a consequence of low current circulation (Bernasconi and Cusminsky, 2005; Bernasconi et al., 2009), high salinity, low oxygen concentrations and low temperatures (Buzas et al., 1993; Nigam et al., 2007; Caulle et al., 2014). This predominant angular-asymmetric morphology also coincides with what is described for other OMZs, such as the coastal margin of western India (Nigam et al., 2007) and east of the Arabian Sea (Kaithwar, 2021).

However, despite reflecting typical OMZ conditions, the record of planispiral, spherical and millioline shells comprise epifaunal species and a rounded-symmetric morphology. According to Manasa et al. (2016) and Khare et al. (2017), the presence of these morphogroups indicates turbulence and disturbance produced by inland water effluents, as well as high sedimentation rates (Saraswat et al., 2018). Therefore, the Mejillones Bay would be reflecting an OMZ with a high energy environment. This may be due to the contributions of wastewater from industrial activity in the area and the presence of cyclonic and anticyclonic whirlpools. These water movements displace upwelling waters (Letelier et al., 2012) and are formed by masses of subantarctic and intermediate Antarctic waters and by the Subtropical Anticyclone of the South Pacific (APS) (Sievers and Silva, 1982; Silva, 1983; Strub et al., 1998).

The highest H′ indices were seen in the stations located to the southwest of the bay, due to the Angamos Point upwelling center, which is one of the highest productivity points in northern Chile (Páez et al., 2001; Valdés, 2004; Valdés et al., 2008). This allows the formation of a large primary productivity deposit in the southern part of the Bay, which has been estimated at a value of 138 mg C m³h^–1 (Rodríguez et al., 1986) and an annual productivity of 1.070 g C m² (Marín et al., 1993). This condition fertilizes the environment and favors the availability of food for foraminifera consumption, indicating that, despite the large amount of organic matter, microxia conditions that could restrict species diversity are not generated.

Multivariate Analyses

The similarity found between the clusters did not indicate a pattern dependent on bathymetry (Figure 11A), and the similarity percentages did not show sufficient support to be reflected in the MDS analysis. This differs from the results found for other OMZs such as the central coast of India (Mazumder and Nigam, 2021) and the central coast of Peru (Cardich-Salazar, 2012), where the formation of clusters reflects a pattern according to depth, differentiating from a shallow water fauna to organisms from deep zones. Evidence of an absence of a depth-dependent pattern could be reflecting turbulence and removal of sediments, caused by the effect of cyclonic and anticyclonic mesoscale whirlpools (formed by intermediate Antarctic waters, the APS and masses of subantarctic and intermediate Antarctic waters). These circular movements avoid a stratification of the waters and therefore a grouping of foraminifera according to the depth.

Vertical Distribution

Taxonomic Diversity

Similar to the results at the spatial level, 93% of foraminifera species belong to the order Rotaliida. This indicates a wide variety of strategies to be able to adhere to various types of substrates. This gives them the advantage of being able to be distributed along the entire sediment column through their trochoid shells (Haynes, 1981; Boix-Martínez, 2007). The most predominant families generally reflected conditions of low oxygen below to 0.3 ml l^–1 and high concentrations of organic matter (Sen-Gupta and Machain-Castillo, 1993; Cardich-Salazar, 2012). Seven species belonged to Bolivinitidae, where Bolivina seminuda had greater representativeness (Figure 5B), mainly in the surface layers (0–3 cm), corresponding to years 1996–2012. This predominance indicated microxia conditions (0,01 ml O₂ l^–1) (Páez et al., 2001), increase in organic matter and lack of affluents from cold and temperate oceanic water masses (Ortlieb et al., 2001). This coincides with the environmental conditions of the seasonal upwelling in the bay (Cardich-Salazar, 2012), and influence of subsurface equatorial currents such as the Perú-Chile countercurrent and the Gunther current (Figueroa et al., 2005), that favors the presence of Bolivinitidae. Therefore, the most representative species of this Family (Bolivina seminuda) has been described for the OMZ located in Concepción (Tapia et al., 2008) and other areas of Chile, such as Coquimbo (Fonseca-Valdebenito, 2018), where it has been reported as the most abundant species in hypoxic conditions. While, the decrease in the predominance of this species (B. seminuda) from the 4–5 cm to the 7–8 cm layer in both cores (Figures 10A,B), may reflect a change in oxygen concentrations (Ortlieb et al., 2001) and a slight oxygenation (Fonseca-Valdebenito, 2018) between years 1877 (7–8 cm) and 1948 (4–5 cm).

The latter coincides with the Buliminella elegantissima abundance, the only species of Buliminellidae. Its increase in lower layers has shown a preference for dysoxic conditions (0,07 ml O₂ l^–1) (e.g., Páez et al., 2001) and infaunal habits, which allows it to inhabit layers with scarce light moving toward the deeper levels (Kitazato, 1981). As well as Fursenkoina sp., which generally inhabits dysoxic substrates with high contents of organic matter related to productivity in the area (Murray, 2006). These species have previously been described and is common in the Mejillones Bay (e.g., Ortlieb et al., 2001; Páez et al., 2001) and the OMZ in the central coast of Peru (Cardich-Salazar, 2012). Unlike the OMZ of conception, where none of these species has been reported. Other OMZs, such as the eastern Arabian Sea (Verma et al., 2018), Pacific Northwest (Sakai et al., 2021) and Gulf of California (Tetard et al., 2017), record Buliminella tenuata, a similar species not recorded for Mejillones Bay. The family Nonionidae, on the other hand, evidenced high abundances mostly represented in the first 2 cm of core Z1A (years 2007 and 2012), while in the lower layers its abundance remained relatively stable (between years 1877 and 1996). This reflects denitrifying activities attributed to severe post-upwelling hypoxic events (Risgaard-Petersen et al., 2006; Cardich-Salazar, 2012; Polovodova-Asteman and Schönfeld, 2015). Finally, Epistominidae indicates a species with detritivore habits, being typical of eutrophied environments (Tyszka, 1994; Sagasti and Ballent, 2002; Smolén, 2012).

Morphogroups and Faunal Association

Compared to the results at the spatial distribution, the percentage of trochospiral shells decreased to 36% (Z1A) and 38% (ZA), while biserial or uniserial species increased to 31% (ZA) and 32% (Z1A). These angular-asymmetric morphologies had a representativeness from the most superficial layer to 4–5 cm. This indicated that between years 1948 (4–5 cm layer) and 2012 (0–1 cm layer), this area represents eutrophied silty sediments with low current circulation. This coincides with the retention zone of the bay, formed by the current pattern of the South Pacific Subtropical Anticyclone (APS) (Strub et al., 1998) and Antarctic subantarctic and intermediate waters (Sievers and Silva, 1982; Silva, 1983), producing a sedimentary deposit. Other OMZs, such as the coastal margin of western India (Nigam et al., 2007) and east of the Arabian Sea (Kaithwar, 2021), show similar conditions and therefore register high angular-asymmetric morphologies.

On the other hand, 26% (Z1A) and 24% (ZA) of species corresponded to planispiral shells, being these dominant rounded-symmetric morphologies between 6 and 8 cm layers. The increase of these morphogroups is related to a Tsunami event in 1877 (6–8 cm layers), caused by an earthquake occurred in the city of Iquique (area close to Mejillones Bay). This Tsunami produced a removal of sediments and strong turbulence that facilitated the presence of epifaunal species with planispiral shells, that commonly live in high energy zones with constant disturbance (Manasa et al., 2016; Khare et al., 2017). Therefore, this natural disaster drastically changed the environmental characteristics of the bay and the benthic foraminifera communities, being this OMZ extremely different from the other OMZs during this year (1877).

The dominance of Bolivina seminuda in 0–1 cm layer of both cores reflected the highest Simpson values and the lowest Shannon diversity values (Figure 8B), indicating that microxic conditions decrease species richness and favor the occurrence of B. seminuda (Sen-Gupta and Machain-Castillo, 1993; Páez et al., 2001). The lowest dominance and highest diversity values between 4 and 8 cm of core Z1A (between years 1877 and 1948) were caused by the increase of most species, among them Buliminella elegantissima, Fursenkoina sp., Nonionella auris, Epistominella exigua, and Nonionoides grateloupi. This evidenced a change in oxygen concentrations, where dysoxic conditions increase diversity in the bay (Páez et al., 2001). Although hypoxic conditions are typical of Mejillones Bay, changes over time are evident, as the layers dated for more recent years (from 0 to 1 cm and 1–2 cm) reveal an intensification in hypoxia. This could be related to anthropic impact on the coast, starting in the late 90s and early 2000s, with the construction of port terminals and thermoelectrics, that continued their development in 2011 (E-Cl, 2010; Universidad de Chile, 2012; ENDESA, 2014).

Multivariate Analyses

The groups formed by the cluster apparently reflected a change in the intensification of hypoxia (microxia: 0,01 ml l^–1) between 2012 and 1996 (groups formed by the 0–3 cm layers), while prior to 1996, the bay presented dysoxic conditions (0,07 ml l^–1) (groups formed by the 4–8 cm layers). This indicates that oxygen is a limiting factor for the vertical distribution of the different populations of foraminifera. Where there is evidence of less richness and abundance of species in microxic sediments (layers of 0–3 cm) than in disoxic sediments (layers of 4–8 cm), where the greater availability of oxygen favors the presence of species. Opposite to other OMZs such as the central coast of Peru (Cardich-Salazar, 2012; Erdem et al., 2019) and Concepción (Tapia et al., 2008), where conditions of higher oxygenation, and higher richness of foraminifera species were found in superficial layers However, the differentiation of clusters between a microxic and dysoxic fauna was not statistically significant. This may be due to local factors that caused removal and turbulence in the bay, such as the mesoscale whirlpools formed by the APS and Antarctic subantarctic and intermediate waters current. The tsunami in 1877 also played a role, which generated a mixture of benthic fauna, avoiding a marked stratification along the sediment column of the bay.

Conclusion

The OMZ of the Mejillones Bay presents differences in its environmental characteristics compared to other OMZs in terms of: (1) taxonomic diversity, being found several species of the genus Bolivina that had not been reported to date in this bay. As well as B. elegantissima, which is not recorded in most OMZs. (2) Although most morphologies in the shells showed eutrophied silty substrates and low current circulation typical of an OMZ, rounded-symmetric morphologies were also recorded, which indicated a turbulent and high-energy environment. (3) There was no evidence of a pattern of spatial and vertical distribution of benthic foraminifera, where no structuring of the communities according to the depth and conditions of microxia or dysoxia was observed. This can be attributed to local bay factors such as the formation of cyclonic and anticyclonic whirlpools, as well as the contributions of wastewater from industrial activity and natural events such as the tsunami of 1877. These factors generate a mixture of benthic fauna and avoid the registration of a spatial zoning and vertical stratification of the sediment. (4) An apparent intensification in the levels of hypoxia (micoxia) in the most superficial layer of the sediment (year 2012) was evidenced, which is possibly related to the industrial development of the bay in the area. These results differ from other OMZs, that show a higher oxygenation in superficial layers compared with deeper sediments that present lower oxygen levels.

The study of the spatial and vertical distribution of the benthic foraminifera of Mejillones Bay, allowed to understand their dynamics and reconstruct their ecology over time. Although it shows typical characteristics of an OMZ, it is an area that presents a set of very particular ecological characteristics different from the other OMZs.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

LT developed the specimen identification of foraminifers and contributed to conceptualization, methodology, software, data analysis, discussion, investigation, writing—original draft, and project administration. MM developed the specimen identification of foraminifers and contributed to data analysis, conceptualization, writing—review and editing, and supervision. PM conducted the methodology and contributed to conceptualization and chemical and data analyses for geochronology and its use for interpretation and discussion. RA contributed to the review, editing, and supervision of the manuscript. All authors listed have made a substantial, direct, and intellectual contribution to the work, and approved it for publication.

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank FONDECYT project 1130511 and Dr. Ruben Escribano from Universidad de Concepción, for obtaining the sediment cores and Prof. Jorge Valdés for obtaining the sediment dredge samples; and the Fondequip EQM 170124-CONICYT project for allowing access to the SU 3500 Hitachi electron microscope at the University of Concepción, Chillán headquarters. And finally, we thank Vicerrectoría de Investigación y Desarrollo-VRID and VRID-investigation project N° 2021000237INV for financing the publication of this manuscript.

Supplementary Material

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

References