Vertical assemblage of the holoplanktonic mollusks (Pteropoda and Pterotracheoidea: Carinaiidae, Pterotracheidae) in the Campeche Canyon, southern Gulf of Mexico, during a “Nortes” season

Study area

The Gulf of Mexico is the largest and deepest interior sea in North America, sharing waters with Mexico, the United States of America, and Cuba (Fig. 1A). The Campeche Canyon in the southern Mexican region of the Gulf is of tectonic origin. It is the deepest geomorphological feature in the region and reaches a depth of more than 3,500 m (Fig. 1B). The Campeche Canyon is highly dynamic and plays an essential role in the marine ecology of the region owing to the occurrence of different hydrodynamic processes, including internal waves (Santiago Arce & Salas de León, 2012), fronts (Aldeco-Ramírez et al., 2009), and eddies (Salas-de-León et al., 2004). These hydrodynamic processes determine the nutrient levels in the euphotic zone that are available for phytoplankton communities (Durán-Campos et al., 2017).

The southern Gulf of Mexico is influenced by northerly winds from September to April (Salas-Monreal et al., 2022). It is called the “Nortes” season. Atmospheric forcing plays a critical role in this region. For example, during the “Nortes” season, extreme (>80 km h⁻¹) and persistent northerly winds cross the gulf, exerting a significant influence on the hydrographic properties of the water column and inducing a deep mixed layer (>70 m depth) (Arriola-Pizano et al., 2022).

Studies on the composition, density, and distribution of zooplankton organisms in the southern Gulf of Mexico began in the middle of the last century because of investigations conducted by the Mexican government, particularly the Mexican Navy. Since then, the importance of holoplanktonic mollusks has been recognized because they are usually the dominant component after the copepods and chaetognaths (Toral Almazán et al., 2022). To date, the study of holoplanktonic mollusks in the Gulf of Mexico has evolved and has been approached in different ways, including taxonomic lists (e.g., Michel & Michel, 1991; Suárez, 1994), studies in neritic areas of the southern portion (e.g., Leal-Rodríguez, 1965; Matsubara-Oda, 1975; Flores-Coto et al., 2013; Lemus-Santana et al., 2014a; Lemus-Santana et al., 2014b), and studies that evaluated the influence of freshwater discharges (e.g., Lemus-Santana, 2011). Most of these studies include pteropods, while heteropods have only been considered in a few cases (e.g., Lemus-Santana et al., 2014b; Espinosa-Balvanera, 2017). Based on these investigations about 53 pteropods and 15 heteropods species were identified.

Sampling

The hydrographic dataset and zooplankton samples considered in this study were recorded during a multidisciplinary scientific expedition “CAÑON-IV” in the Campeche Canyon during February 21–28, 2011, on board of the R/V Justo Sierra owned by the Universidad Nacional Autónoma de México (UNAM). Hydrographic data were acquired at 48 stations equidistantly separated every 14 nautical miles (+ symbols in Fig. 1B: Table S1) using a conductivity-temperature-depth probe (CTD, SeaBird-19 plus) configured with oxygen and chlorophyll-a fluorescence sensors (SBE-43 and ECO-Wet Labs, respectively) which were calibrated pre- and post-cruise and placed in a Rosette (General Oceanics). Each CTD cast was taken from the sea surface to close to the bottom (5 m above the sea floor). The probe was lowered at a rate of 1 m s⁻¹ with the interface configured to store data at 24 Hz.

Immediately following each Rosette-CTD cast, zooplankton samples were collected, during both daytime and nighttime, from a total of 24 stations (O symbols in Fig. 1B), at four depths (10, 50, 100, and 200 m), using a double trip mechanism close-open-close system (General Oceanics) and conical nets (75 cm diameter, 250 cm total length, and 505 µm mesh size) configured with a flowmeter (General Oceanics, model 2030R) placed in the mouth of each of the nets. A vertical array of four systems was placed on a hydrographic winch wire, and the cosine of the wire angle was calculated to ensure the capture of organisms at each depth (Kramer et al., 1972). Each net was opened by manual messengers starting the haul at 2 knots for 15 min; once the hauls finished, the nets were closed again with manual messengers and then recovered. Once on board, the organisms collected were fixed with formalin (added sodium borate) at 4% for 24 h. The final preservation was in 70% ethanol in glass jars with tight lids. The samples were stored under special precautions to avoid degradation of the organisms until analysis; for example, we usually changed both the airtight lids and the ethanol every two months.

Laboratory analyses

Holoplanktonic mollusks were obtained from 65 samples collected at 24 sampling stations in the Campeche Canyon. Complete and undamaged specimens of both groups were picked and separated in a glass Petri dish; approximately 20,000 specimens were sorted. They were identified at the species level using a Carl Zeiss Stemi 508 stereomicroscope equipped with a Zeiss Axiocam 105 color following the standard keys of Richter & Seapy (1999); Van der Spoel & Dadon (1999); Gasca & Janssen (2014), and Wall-Palmer et al. (2018). The identification was compared and confirmed using online open-access databases (e.g., MolluscaBase, Biodiversity Heritage Library, Tree of Life web project, and local repositories, such as the National Commission for the Knowledge and Use of Biodiversity (CONABIO, Spanish acronym)). The number of organisms was finally standardized to density units as number of individuals per 100 cubic meters (ind 100 m⁻³) using registers of flowmeters placed in the mouth of each net, following the standard protocols of Kramer et al. (1972), and Harris et al. (2000). After adult organisms’ separation, we realized taxonomic determination at the lowest possible level using specialized literature as Burridge et al. (2019), Gasca & Janssen (2014), Janssen, Bush & Bednarsek (2019), Richter & Seapy (1999), Tesch (1950); Van der Spoel, Bleeker & Kobayasi (1993). The main characteristics analyzed were the shell shapes and ornamentations (keels, notches, tubercles, reticulations, number of ribs and striation, and spines). For gymnosomes, the anatomy of the gills, pedal lobes, and oral structures (cones, arms, and suckers) were considered. The shape of the visceral nucleus was analyzed for heteropod mollusks of the Pterotracheidae and Carinariidae families, as they are organisms with sexual dimorphism, we analyzed the presence of reproductive structures such as gonads, tentacles, and suckers. The specimens were photographed from various angles using the manufacturer’s software (Zeiss ZEN Core) to visually document the species. Subsequently, the images were processed using Helicon Focus software version 8.2.0 and edited using Adobe Photoshop software version 23.5.1.

Data reduction

The CTD data were initially processed using the routines included in the software provided by the manufacturer (SBE Data Processing v.7.26.7). This involved removing poor-quality data and averaging the remaining data to 1 dbar. Then, conservative temperature (°C), absolute salinity (g kg⁻¹), and density (sigma-t, kg m⁻³) were calculated using the algorithms of the thermodynamic equation for seawater, TEOS-10 (IOC, SCOR and IAPSO, 2010). These data were then used to construct vertical cumulative profiles to identify the thermocline/pycnocline depth, which was later confirmed using the depth of the maximum temperature vertical gradient (∂T/∂z). Geostrophic velocities relative to a depth of 1,000 m or to the bottom at shallower stations were calculated following the standard methods described by Pond & Pickard (1995). Finally, the vertical section shows the hydrographic parameter distributions constructed along transect T-T’ (Fig. 1B).

The dataset generated with the physical-biological parameters was organized in matrices that included the density value of each species at each sampling depth and the corresponding physical variables (e.g., temperature, salinity, and current speed). The organism density values at each sampling depth are represented using box-and-whisker plots. The relative abundance of each taxon was ordered in bar graphs stacked at 100% to analyze the variation in the composition and organism’s density values of holoplanktonic mollusks at each sampling depth (10, 50, 100, and 200 m). These graphs were realized with the “ggplot2” library (Wickham, 2016) in the R Studio software (R Core Team, 2021).

The diversity index (H’) at each sampling station for each sampling depth was calculated following Shannon & Wiener (1949) using the “biodiversityR” (Kindt & Coe, 2005) and “vegan” (Oksanen et al., 2022) libraries in R Studio software (R Core Team, 2021). These values were represented in graphs for each sampling depth using box-and-whisker plots built with the “ggplot2” library (Wickham, 2016) in the R Studio software (R Core Team, 2021). The Kruskal-Wallis test was used to analyze the significance of the differences in density, species richness, and diversity about depth.

Finally, principal component analysis (PCA) was conducted to analyze the physical variables affecting holoplanktonic mollusk assemblages at each sampling depth. Additionally, to identify similarities between the holoplanktonic mollusk communities during sampling, non-metric multidimensional scaling (NMDS) and SIMPROF tests were applied. These analyses were performed following the standard routines of the Plymouth Routines in the Multivariate Ecological Research (PRIMER, v7) software (Clarke, 1993). Finally, we realized a Pearson correlation for to analyze the statistical significance between density of each dominant species and environmental variables included in the PCA.

Hydrography

The vertical distribution of the hydrographic parameters in the upper 200 m shown in cumulative profiles indicated that the mean profiles of conservative temperature ranged from 24 to 14 °C, with the maximum vertical gradient observed at 90 m depth (Fig. 2A). Absolute salinity values varied from 35.8 g kg⁻¹ at the surface and increasing to 36.75 g kg⁻¹ at 70 m depth, then decreasing to 36.0 g kg⁻¹ (Fig. 2B). Density (sigma-t) varied from 24.1 to 26.7 kg m⁻³ (Fig. 2C). The thermocline/pycnocline depth during sampling period was at 90 m (Fig. 2). This was estimated by the maximum vertical temperature gradient (∂T/∂z).

The horizontal distribution of hydrographic parameters at a thermocline/pycnocline depth (90 m) showed features of warm and cold cores in the study area. The conservative temperature ranged from 18.7 to 21.5 °C displaying a large cold tongue in the southern region, a warmer northern region, and the presence of two fragmented cores with different temperatures in the central portion (Fig. 3). The absolute salinity ranged from 36.5 to 36.7 g kg⁻¹ (see Fig. S1). The presence of a northern region with relatively higher salinity, as well as a small tongue of lower salinity in the southern region were observed. However, in general it was very uniform along the domain. The density values ranged from 25.8 to 26.2 kg m⁻³; as expected, their distribution displayed a similar pattern to those of the conservative temperature, with cores of different density values (Fig. 4). The circulation pattern showed westward geostrophic currents in the northern portion and southward currents in the eastern portion, reaching values of approximately 30 cm s⁻¹. The cyclonic circulation in the southern region coincides with the tongue of low temperatures and salinity. Additionally, the presence of two small eddies, one cyclonic eddy (centered at 20.9°N and 93.2°W) and an anticyclonic eddy (centered at 20.3°N and 92.6°W) (Fig. 4), agree with the high and low densities, respectively. In the southwestern study area, horizontal density distribution at 90 m depth, at the beginning of the TT’ transect (Sta. 37; see Fig. 1B), shows a pair of cores low-density and high density, this last with a cyclonic geostrophic circulation. Generally, on the boundary between them there is a thermohaline front. On the other hand, at station 17 (see Fig. 1B) a high-density core is observed (Fig. 4).

The vertical distribution of the conservative temperature, absolute salinity, and sigma-t confirmed that the thermocline, halocline, and pycnocline were at a depth of 90 m with variations along the transect (Figs. 5A–5C). On the other hand, at station 17 (see Fig. 1B), from 200 m to approximately 80 m depth, the uplift of isolines suggests the occurrence of the cyclonic eddy (Figs. 5A–5C). In the upper 40 m layer, the variables show notable changes close to Sta. 37, indicating the front effect.

Holoplanktonic mollusks assemblages

From a total of 19,964 organisms recorded, only 16,235 were identified at the species level belonging to two orders, three suborders, 15 families, 27 genera, and 36 species. Among them, 33 were order Pteropoda, of which four (Pneumodermopsis macrochira Meisenheimer, 1905, Spongiobranchaea intermedia Pruvot-Fol, 1926, Schizobrachium cf. polycotylum Meisenheimer, 1903, and Cliopsis krohnii Troschel, 1854) were recently identified, in a collateral study, as new records for the southern Gulf of Mexico (López-Cabello et al., 2024). Some specimens identified in this study are shown in Fig. 2. Three species were of the superfamily Pterotracheoidea (Carinaiidae, Pterotracheidae) (Table 1). The pteropods Heliconoides inflatus, Creseis conica, and Limacina trochiformis presented the highest density values (2,228.5, 1,541.7, and 1,416.7 ind 100 m⁻³, respectively), representing 69.1% of the total density. In contrast, the pteropods Clio cuspidata, Cymbulia peronii, Clione limacina, and the heteropod Pterotrachea hippocampus presented values of <0.3 ind 100 m⁻³, the lowest among all samples (Table 1). C. conica is dominant at 10 m depth, while H. inflatus is dominant at 50, 100 and 200 m when comparing the density of the species at different depths. The species’ dominance at different depths is such that at 10 m, C. conica has the highest density, followed by H. inflatus and finally L. trochiformis. At 50 m, the species with the highest density is H. inflatus, followed by C. conica and finally L. trochiformis. At 100 and 200 m, the hierarchy in terms of dominance is H. inflatus, L. trochiformis, and C. conica. Considering the total density of each species at each of the four depths, the dominant species is H. inflatus, followed by C. conica, and L. trochiformis (Table 1).

The density, species richness, and diversity of total species differ depending on the depth of the sampling. These differences in depth are significant according to the Kruskal-Wallis test, showing p < 0.05. For example, the highest density values were found at 10 m depth (median = 127.2 ind 100 m⁻³, mean = 205.7 ind 100 m⁻³), followed by the 50 m depth level (median = 68.2 ind 100 m⁻³, mean = 161.6 ind 100 m⁻³), whereas the 100 m and 200 m levels presented lower density values (Fig. 6A). The highest species richness was recorded at a depth of 100 m (28 species), whereas the lowest richness value was observed at a depth of 200 m (21 species). At depths of 10 and 50 m, the richness was 25 and 26 species, respectively (Fig. 6B). The diversity values (H’) also showed variations depending on the sampling depth, with the highest value at 100 and 200 m depth (>1.9 bits ind⁻¹), while 10 m depth presented the lowest value (<1.5 bits ind⁻¹) (Fig. 6C). The vertical distribution of the species that showed the highest density values was explored along T-T’ transect (see Fig. 1B) in relation to the hydrographic parameters. The vertical distribution of H. inflatus along this transect showed that the values varied in a range of 1.3 to 255.2 ind 100 m⁻³, with the maximum one observed at 10 m depth in the station 37. The second-highest relative value was observed at a depth of 200 m at Station 35 (Fig. 7A). C. conica ranged from 0.4 to 255.7 ind 100 m⁻³, with a similar vertical distribution pattern recording the highest values at 10 m depth in station 37 (Fig. 7B). L. trochiformis values ranged from 0.4 to 162 ind 100 m⁻³ (Fig. 7C), showing two maximum peaks located in the surface waters (10 m depth), the maximum in station 17 and a second peak in station 37.

Statistical analyses

The PCA analysis showed that the two first axes explained 53.1% of the accumulated variance. In the first principal component (PC1), dissolved oxygen, conservative temperature, and absolute salinity were the variables with the highest influence, whereas in the second principal component (PC2), longitude had the highest influence (Table 2). Pearson correlation analysis between the density of each dominant species and environmental variables registered in this study showed significant correlation (p < 0.05), between Heliconoides inflatus and longitude, between Lamicina Trochiformis and sampling depth, temperature, and dissolved oxygen and, between Creseis conica and sampling depth, and temperature.

The PCA diagram (Fig. 8A) shows that the data obtained at a depth of 200 m were opposite to chlorophyll-a fluorescence, conservative temperature, absolute salinity, and geostrophic speed. The opposite pattern was observed in the data obtained at a depth of 10 m, where dissolved oxygen, conservative temperature, absolute salinity, and chlorophyll-a fluorescence strongly influenced stations 1, 3, 5, 17, 29, and 31. The PCA diagram also shows that some stations at different depths (31, 39, 41, and 42) in the southeastern region presented a relationship with longitude.

The NMDS (with a stress of 0.13) and SIMPROF analyses (Fig. 8B) showed seven groups, the first two corresponding mainly to the stations at a depth of 10 m. Two additional groups were included at the depth of 50 m. The fifth and sixth groups included stations at 100 m depth, whereas the seventh group included stations at 200 m depth. This suggests that the composition and density of the holoplanktonic community in Campeche Canyon, at the time of our observations, presented differences depending on the depth at which they were collected.