Introduction

Planktonic dispersal in the marine environment is one of the leading connectivity processes among populations of meroplanktonic and holoplanktonic gastropods (^{Crocetta et al., 2020}; ^{Roegner, 2000}), and is a relevant element of the trophic dynamics of ecosystems and indirect indicators of the health of adult populations (^{Campos & Landaeta, 2016}; ^{Chávez-Villegas et al., 2014}). Planktonic dispersal study will aid in a better understanding of the molluskan community’s behavior under present and future environmental conditions (e.g., ^{Aceves-Medina et al., 2020}; ^{Molina-González et al., 2018}). For this purpose, baseline studies on the species composition that rely on the use of genetic markers are a starting point.

Planktonic gastropod identification is usually made using their morphology through microscope observation. The use of genetic markers (the Cytochrome Oxidase Subunit I of the mitochondrial DNA, COI) for this purpose is more recent. It has been used to study the genetic connectivity of the marine caenogastropod Bursa scrobilator (^{Crocetta et al., 2020}). Due to its high conservation, the 18S ribosomal nuclear DNA (18S rDNA) gene has been proposed for biodiversity studies (^{Ranjithkumar et al., 2018}).

Planktonic gastropods have been quite extensively studied worldwide (^{Bandel et al., 1997}; ^{Campos & Landaeta, 2016}; ^{Chávez-Villegas et al., 2014}; ^{Oliva-Rivera & Navarrete, 2007}). However, in northwestern Mexico, the main efforts have been in the study of the species composition of holoplanktonic mollusks in the Gulf of California (^{Angulo-Campillo et al., 2011}) and the western coast of the Baja California Peninsula (BCP) (^{Aceves-Medina et al., 2020}; ^{Molina-González et al., 2018}; ^{Moreno-Alcántara et al., 2020}; ^{Sánchez-Hidalgo y Anda, 1989}), characterizing these regions as highly diverse. However, there is no available information about meroplanktonic gastropods in the region. Therefore, this work presents the first record of the planktonic gastropod biodiversity at 2 sites of the western coast of the BCP obtained by genetic markers.

Materials and methods

The sampling design originally intended to evaluate the presence of abalone larvae at 2 sites on the western coast of the Baja California Peninsula: La Bocana (LB: 26°45’1.9” N, 113°42’20.9” W) and Cabo Tosco (CT: 24°18’26.2” N, 111°42’36.6” W) (Fig. 1). However, we focused our study on other mollusks, as abalone larvae were not found. Sampling was done on December 1^st, and 3^rd, 2009 for LB and CT, respectively. At LB, the sampling area was located about 4.5 km southward, with an average depth of 10 m (Fig. 1). The rocky bottom floor is mostly flat with several small stone promontories composed of boulders, crevices, and a high abundance of macroalgae, mainly Eisenia arborea (Areschoug, 1876). CT is located at the southernmost part of Isla Margarita, and the sampling site was just 150 m off the coastline (Fig. 1), with a depth of approximately 3 m. The bottom is composed of large rocks with abundant crevices, surrounded by soft bottom but low macroalgae coverage. These sites were selected because they represent a high abalone productive area (LB) and the southern extreme of abalone species distribution (CT) (^{León-Carballo & Muciño-Díaz, 1996}). The 8-day mean surface temperature in the first week of December 2009 at each site, which includes the sampling period and was obtained from the Aqua MODIS satellite -NPP, 0.0125°, West U.S., day time (11 microns), ^{Simons (2022)}- was 21 °C and 24 °C, respectively (Fig. 1).

Figure 1 Sampling sites at La Bocana (LB) and Cabo Tosco (CT) on the Baja California Peninsula. Color scale represents the 8-day sea surface temperature (°C) data from Aqua MODIS satellite (^{Simons, 2022}). 

At LB, there were 3 towing stations: 1) surface and 2) at 3-4 m depth at the same geographical position (LB-01 and LB-10, respectively), and 3) at the surface, approximately 1 km to the west (LB-07) of the previous position. At CT, there was one towing station (CT-01) at the surface, in a semi-protected area on the southwestern side of Isla Margarita. Tows followed a circular trajectory at a speed of 2.5 knots for 5 min. The plankton net was conical (0.45 m mouth diameter and 180 μm mesh size), with a calibrated flowmeter and a collection cup also with a 180 μm mesh size. Samples were transferred to a 500 mL plastic flask and preserved in 80% ethanol.

The biomass from the plankton tows was estimated by displaced volume (^{Beers, 1976}) and standardized to mL of biovolume per 1,000 m³ of filtered seawater (^{Smith & Richardson, 1979}). Due to high biomass at each towing station, the samples were fractionated to ¼ with a Folsom splitter. Using a stereoscope at 35×, following the description of ^{Courtois de Viçose et al. (2007)} for Haliotis tuberculata abalone larvae, gastropod mollusks with egg or globular shapes were separated (^{Bandel et al., 1997}). A total of 384 individuals were separated, set in slides with wells (1 individual per well), and photo-documented (Leica Application Suite EZ, ver. 1.5.0). Individuals were transferred to 96-well PCR plates containing 70% ethanol.

DNA extraction was done following the protocol of ^{Selvamani et al. (2001)} , with slight modifications. Before extraction, individuals were crushed with a plastic pestle and dried at ambient temperature. Samples were incubated in 18 μL of lysis buffer (Tris-HCl 10 mM, pH 8.3, KCl 50 mM, Tween-20 0.05%) and 2 μL de Proteinase-K (Sigma®, 2.5 mg/mL) at 55 °C for 24 h. After this time, the proteinase was inactivated by incubating the samples at 95 °C for 10 min. The plates were preserved at -20 °C until the PCR was performed.

A fragment of the 18S rDNA gene was amplified with primers 18SHal2F (5’-TTGGATAACTGTGGTAATTCT AGAGC) and 18SHal2R (5’-CCGGAATCGAACCCT GAT) (^{Aranceta-Garza et al., 2011}). The reactions contained 0.2 mM de dNTPs (Invitrogen), 0.5 μM of each primer, PCR buffer 1× (Invitrogen), 1.5 mM MgCl₂, 0.05 U Platinum® Taq DNA polymerase (Invitrogen), 2 μL of DNA and completed to 30 μL with MilliQ® water. PCR conditions were: 1 cycle at 94 °C for 4 min, 30 cycles at 94 °C for 45 s, 57 °C for 45 s, 72 °C for 45 s, and a final extension at 72 °C for 10 min in a thermal cycler (C1000, Bio-Rad). PCR products were verified in a 1.5% agarose gel stained with SYBR gold (Promega), visualized in a UV transilluminator (UVP BioDoc-it Imaging System).

PCR products of approximately 220 bp were individually sequenced either in forward and reverse (LB-7) or in forward only (LB-01, LB-10, CT-01) (Macrogen,Korea). The sequences were manually revised for quality and edited with ChromasPro v. 2.1.9 (Technelysium Pty Ltd). The edited sequences were analyzed to determine the number of similar genotypes with the software DnaSP v. 6.12.03 (^{Rozas et al., 2017}). Identical sequences (100% match) were organized into Molecular Operational Taxonomical Units (MOTUs), for which the putative identity was obtained by a local alignment analysis with BLAST (^{Altschul et al., 1990}). MOTUs were identified either to family, genus, or species level based on the hits with the highest similarity percentage (> 95%) and lowest E value (probability of the sequence match being random < 1 × 10^-90). Nomenclature follows the NCBI (National Center for Biotechnology Information) Taxonomy database (^{Schoch, 2011}).

The MOTUs sequences were analyzed as haplotypes (most sequences were homozygotes) to construct a minimum spanning network (^{Bandelt et al., 1999}). The software PopArt (^{Leigh & Bryant, 2015}) was used to visualize genealogical relationships among towing stations and obtain genetic diversity estimators (nucleotide diversity and number of segregating and parsimony-informative sites). The MOTUs’ diversity among towing stations was estimated by the Shannon-Weaver index using the natural logarithm with the vegan package (^{Oksanen et al., 2020}) in R ver. 4.1.1 (^{R Core Team, 2013}), using MOTUs as putatively non-identified species.

Results

At LB, the plankton biomass at the mid-water tow was more than twice as high (LB-10, 1,222 mL 1,000 m^-3) as at the 2 surface stations (LB-01 and LB-07, 491 and 557 mL 1,000 m^-3, respectively). Surface biomass at CT (CT-01) was 883 mL 1,000 m^-3. Planktonic gastropods, similar in shape, varied between 200 - 500 μm in width. As shells were translucent, it was possible to assess the presence of tissue within (Fig. 2).

Figure 2 Gastropod mollusks obtained from a plankton tow from La Bocana. Scale = approx. 200 μm. 

Of a total of 384 individuals analyzed, 302 resulted in reliable DNA sequences matching gastropod mollusks (phylum Mollusca, class Gastropoda). The rest of the sequences were of low quality (n = 75) or matched other taxonomic groups (phylum Cnidaria, order Bivalvulida, and phylum Arthropoda, orders Euphausiacea and Calanoidea, n = 6) or showed no significant match to any species (n = 1). A total of 71 MOTUs were distinguished with matching identities higher than 95% and E-values lower than 1 × 10^-100 in most sequences, putatively belonging to several families and orders of 5 subclasses (Tables 1, 2).

From the similarity analysis to GenBank sequences by BLAST, 23 MOTUs gave 100% identity. Identification to the species level was likely for 7 MOTUs (i.e., MOTUs 02, 05, 13, 33, 48, 61, and 62; Table 2). The remaining 16 are not certain as some of these MOTUs showed matching hits either to several genera (i.e., MOTUs 01, 03, 08, 10, 12, 15, 20, 22, 25, 28, 64, and 68; Table 2) or several species of a certain genus (i.e., MOTUs 29, 65, 66, and 67; Table 2).

When individual images were organized per MOTU (Supplementary material: Fig. S1[https://doi.org/10.5281/zenodo.7685921]), it was observed that there are many similarities among individuals at distinct taxonomic levels, including subclasses. In contrast, a high genetic diversity was found (nucleotide diversity π = 0.045, with 72 segregating sites and 55 parsimony-informative sites; Supplementary material: Table S1. The phylogenetic analysis inferred from the Minimum Spanning Network indicated that most identified subclasses were represented in all locations (Fig. 3). One exception was Neritimorpha, which was mainly observed in CT.

Figure 3 Minimum spanning network from 18S rDNA sequences per MOTUs. Background shades represent gastropod subclasses. 

The MOTUs’ diversity (Shannon index) varied among towing stations with the highest at LB-10 and the lowest at (LB-07) (Fig. 4). This large difference can be inferred by the presence of many MOTU 04 (family Cypraeidae) individuals at LB-07 (Table 2). Even though the diversity was rather similar among CT-01, LB-01, and LB-10 (Fig. 4), the MOTU composition was different with several MOTUs found only in a single towing station (CT-01: 17; LB-01: 16; LB-10: 14; LB-07: 7), i.e., 54 of the 71 MOTUs found were unique (Table 2).

Figure 4 MOTUs’ Shannon-Weaver diversity index, using the natural logarithm, per towing station. 

Discussion

A relatively high planktonic gastropod diversity was observed in the present study as compared to that found in 2 sites of the Caribbean Sea: 31 species (^{Oliva-Rivera & Navarrete, 2007}) and 34 species (^{Chávez-Villegas et al., 2014}), but lower than the 62 holoplanktonic species from the Gulf of California (^{Angulo-Campillo et al., 2011}), the 6 species of the order Thecosomata (^{Sánchez-Hidalgo y Anda, 1989}), and the 18 species of the family Atlantidae off the BCP (^{Aceves-Medina et al., 2020}; ^{Moreno-Alcántara et al., 2020}). Nevertheless, according to a preliminary rarefaction analysis (results not shown due to reduced sample size) the actual species diversity in our study region is probably under-represented. To overcome the sample size issue, which is mainly due to budget constraints, new research technologies such as environmental DNA (eDNA) metabarcoding are broadening the sampling possibilities (^{Mychek-Londer et al., 2020}).

Differences in community composition between the sampling sites of CT and LB (Fig. 1) would be expected as environmental conditions are usually distinct. Personal observations indicated that both sites are mainly composed of rocky bottoms and a few sandy flats, with a higher abundance of macroalgae at LB than at CT, and clearer waters at CT than LB. These differences, which are probably due to a higher productivity at LB, are explained by the position of the sites relative to the Magdalena Biological Activity Center (MBAC) (^{Lluch-Belda, 2000}), with year-round cooler waters at LB (north of MBAC) than at CT (south of MBAC) (^{Jerónimo & Gómez-Valdés, 2006}). In contrast, the differences in MOTUs composition among stations within LB let us hypothesize that planktonic gastropod distribution within a specific area is probably modulated by micro-environmental processes such as currents and biological events (^{Hernández-Trujillo et al., 2001}). For example, many Cypraeidae individuals at LB-07 points to a massive spawning event in that area.

Even though 18S rDNA has been proposed for biodiversity studies (^{Ranjithkumar et al., 2018}), it generally fails to separate species (^{Aranceta-Garza et al., 2011}; ^{Wu et al., 2015}). In addition, the lack of a comprehensive database of gastropod 18S rDNA sequences from the western coast of the Baja California Peninsula limits the capacity of proper species identification. For example, the report of ^{Angulo-Campillo et al. (2011)} indicates the presence of several Atlantidae or Desmopteridae species, which in our case could not be distinguished (Table 2). Thus, future studies will require additional genetic markers, such as the 16S rRNA or COI genes (e.g., ^{Crocetta et al., 2020}).

The putative presence of molluskan MOTUs not recorded in the study region can be explained by either the low number of 18S rDNA sequences in the public databases, the lack of differences among species in this gene sequence of the same genus or family, and the lack of physical records of those specimens in the study region. In any case, efforts should be devoted to increase the number of gastropod 18S rDNA sequences and other genes in the databases and to promote frequent survey campaigns to detect previously undescribed species in this region. This is particularly relevant because mollusks are gaining attention within climate change research, not only to study phenotypic plasticity and genetic variation (^{Matoo & Neiman, 2021}), but also because of their utility to act as indicators of climatic variability (^{Molina-González et al., 2018}). Also, eDNA metabarcoding, which can be used to detect diverse ecological phenomena such as potential species invasions (e.g., ^{Mychek-Londer et al., 2020}), requires solid supporting information.

The presence of non-gastropod sequences within the analyzed individuals is a phenomenon that requires further attention. Even though we cannot rule out the possibility of cross-contamination, these individuals were visualized as planktonic gastropods (Supplementary material: Fig. S1), so it is unlikely that individuals of Cnidaria (order Bivalvulida), and Arthropoda (orders Euphausiacea, and Calanoida) were mistakenly put into the PCR tubes. We also consider that their presence as stomach content or as epibionts was unlikely because when mixed DNA is individually sequenced, it is not possible to separate the signature of each species. Rather, we hypothesize that these sequences might come from eggs or larvae of those other taxa that use gastropod shells as a substratum.