Oceanic dispersal barriers in a holoplanktonic gastropod

Abstract Pteropods, a group of holoplanktonic gastropods, are regarded as bioindicators of the effects of ocean acidification on open ocean ecosystems, because their thin aragonitic shells are susceptible to dissolution. While there have been recent efforts to address their capacity for physiological acclimation, it is also important to gain predictive understanding of their ability to adapt to future ocean conditions. However, little is known about the levels of genetic variation and large‐scale population structuring of pteropods, key characteristics enabling local adaptation. We examined the spatial distribution of genetic diversity in the mitochondrial cytochrome c oxidase I (COI) and nuclear 28S gene fragments, as well as shell shape variation, across a latitudinal transect in the Atlantic Ocean (35°N–36°S) for the pteropod Limacina bulimoides. We observed high levels of genetic variability (COI π = 0.034, 28S π = 0.0021) and strong spatial structuring (COI ΦST = 0.230, 28S ΦST = 0.255) across this transect. Based on the congruence of mitochondrial and nuclear differentiation, as well as differences in shell shape, we identified a primary dispersal barrier in the southern Atlantic subtropical gyre (15–18°S). This barrier is maintained despite the presence of expatriates, a gyral current system, and in the absence of any distinct oceanographic gradients in this region, suggesting that reproductive isolation between these populations must be strong. A secondary dispersal barrier supported only by 28S pairwise ΦST comparisons was identified in the equatorial upwelling region (between 15°N and 4°S), which is concordant with barriers observed in other zooplankton species. Both oceanic dispersal barriers were congruent with regions of low abundance reported for a similar basin‐scale transect that was sampled 2 years later. Our finding supports the hypothesis that low abundance indicates areas of suboptimal habitat that result in barriers to gene flow in widely distributed zooplankton species. Such species may in fact consist of several populations or (sub)species that are adapted to local environmental conditions, limiting their potential for adaptive responses to ocean changes. Future analyses of genome‐wide diversity in pteropods could provide further insight into the strength, formation and maintenance of oceanic dispersal barriers.


| INTRODUC TI ON
Anthropogenic carbon emissions are an important cause of perturbations in marine ecosystems, including an overall increase in seawater temperature, ocean deoxygenation and acidification Gruber et al., 2019). These perturbations are expected to have wide-ranging impacts on marine organisms and may drive species range shifts and ecosystem regime shifts due to climate-driven invasions and extinctions (Doney et al., 2012;Kroeker et al., 2013;Poloczanska et al., 2016). It is uncertain whether marine organisms have the short-term plasticity and/or long-term evolutionary potential to adapt, given the ongoing and expected rates of environmental change (Donelson et al., 2019;Miller et al., 2018). To gain insight into the evolutionary potential of marine taxa, it is necessary to resolve their underlying genetic variation, population structure and connectivity (Bell, 2013;Harvey et al., 2014;Munday et al., 2013;Poloczanska et al., 2016;Sunday et al., 2014).
Locally adapted populations are increasingly reported in marine systems, showing that the apparent connectivity of the marine environment is often not fully realized by organisms, and populations may be more spatially constrained than once thought (Barth et al., 2017;Nielsen et al., 2009;Sanford & Kelly, 2011). Several case studies have illuminated the processes involved in local adaptation of marine molluscs, for example habitat preference and partial spawning asynchrony in the mussels Mytilus edulis and Mytilus galloprovincialis (Bierne et al., 2003), ecological selection across microhabitats in the ecotypes of the rocky shore gastropod Littorina saxatilis (Butlin et al., 2014;Johannesson et al., 2017;Westram et al., 2018), or variation in thermal stress tolerance in the intertidal gastropod Chlorostoma funebralis (Gleason & Burton, 2016). Much of the existing literature has focused on dispersal barriers in marine benthic taxa, some of which have pelagic larval stages. However, little is known about evolutionary pressures acting in the pelagic environment. One reason is that it is difficult to differentiate between selection pressures acting on the larval pelagic and adult benthic life stages on the same genome (Marshall & Morgan, 2011). Thus, it could be useful to study wholly pelagic organisms, such as holozooplankton, to better understand the selective forces acting upon organisms in the open ocean.
Population structure in marine holozooplankton appears to be more influenced by their ability to establish and maintain viable populations outside of their core range, rather than dispersal limitations (De Vargas et al., 1999;Norton & Goetze, 2013;Peijnenburg et al., 2004). Given the large population sizes, extensive distribution patterns, high standing genetic diversity, and short generation times of marine zooplankton in general, adaptive responses to even weak selection may be expected as the loss of variation due to genetic drift should be negligible (Peijnenburg & Goetze, 2013). While some plankton communities show rapid range shifts and phenological changes in response to ocean warming (Beaugrand et al., 2009(Beaugrand et al., , 2014Hays et al., 2005), we only have limited understanding of the potential for evolutionary responses of marine zooplankton to future conditions (Dam, 2013;Peijnenburg & Goetze, 2013).
Pteropods are a group of holoplanktonic gastropods that play important roles in planktonic food webs and ocean biogeochemistry, due to their high abundance and production of calcium carbonate shells (Bé & Gilmer, 1977;Buitenhuis et al., 2019;Hunt et al., 2008). Shelled pteropods are especially vulnerable to global change, because of their thin aragonitic shells that are prone to dissolution Lischka et al., 2011;Manno et al., 2017), and they have been suggested as bioindicators of the effects of ocean acidification . Several prior genetic studies on shelled pteropods have focused on resolving species boundaries via an integrative approach, combining genetic analyses and morphometric measurements of shells, to obtain a better understanding of species distribution patterns (Burridge et al., 2015(Burridge et al., , 2019Shimizu et al., 2018). Other recent research efforts have addressed the response of shelled pteropods to acidified conditions to evaluate their acclimation to rapid increases in ocean acidity Bogan et al., 2020;Maas et al., 2018;Moya et al., 2016). Little is known, however, about the spatial distribution of natural genetic and phenotypic variability in pteropod species with widespread distributions.
Planktonic ecosystems in the Atlantic Ocean have been relatively well-sampled owing to spatially extensive (~13,500 km) annual transects of the Atlantic Meridional Transect programme (AMT, https:// www.amt-uk.org), and are thus ideal for studying population structure and dispersal barriers in marine zooplankton. The AMT cruises traverse both the northern and southern subtropical gyres, which are separated by the equatorial upwelling region. The subtropical gyres are oligotrophic systems characterized by clear ocean waters with very low primary productivity in the surface layer, high sea surface temperature, a deep thermocline, and the presence of a deep chlorophyll maximum. Conversely, the mesotrophic equatorial province is characterized by upwelling of nutrient-rich deep waters that stimulate primary production. Several species of copepods show spatial genetic structuring congruent with these Atlantic oceanic provinces.
For example, Pleuromamma abdominalis has clades that are endemic to the equatorial province (Hirai et al., 2015). Other mesopelagic copepods, such as Haloptilus longicornis (Andrews et al., 2014;Norton & Goetze, 2013) and Pleuromamma xiphias  show strong genetic breaks between the northern and southern subtropical gyre populations. These dispersal barriers coincide with regions of low abundance, which may represent areas of suboptimal habitat, supporting an ecological basis for these barriers .
However, it is unclear whether other zooplankton groups show similar patterns of basin-scale genetic structuring. For holoplanktonic gastropods, DNA barcoding has revealed intra-specific genetic differentiation between and within ocean basins, pointing towards undescribed diversity (Burridge, Hörnlein, et al., 2017;Jennings et al., 2010;Wall-Palmer et al., 2018). Using DNA barcoding, two closely related species of Protatlanta (a genus of Atlantidae, commonly known as shelled heteropods) were identified in the Atlantic Ocean: P. souleyeti, which is abundant in the subtropical gyres, and P. sculpta, which is most abundant in the equatorial upwelling region (Wall-Palmer et al., 2016). There is also some evidence that populations of the straight-shelled pteropod Cuvierina atlantica show genetic discontinuity across the equatorial upwelling region in the Atlantic Ocean (Burridge et al., 2015).
In this study, we aim to identify barriers to dispersal and spatial population structure in the coiled-shelled pteropod L. bulimoides (d'Orbigny, 1835) across a latitudinal transect in the Atlantic Ocean. This species has a circumglobal warm-water distribution from ~45°N to ~40°S, a preferred depth range of 80-120 m, and it performs diel vertical migration with higher abundance in surface waters at night (Bé & Gilmer, 1977). In the Atlantic basin, the species is common across several ocean provinces with peaks in abundance in the oligotrophic gyres (Bé & Gilmer, 1977;. By assessing the distribution of genetic and phenotypic variability in L. bulimoides across a latitudinal Atlantic transect, we aim to test the hypothesis that areas of low abundance mark regions of suboptimal habitat that represent barriers to gene flow. To do this, we sequenced fragments of the mitochondrial cytochrome c oxidase I (COI) and the nuclear 28S genes and made geometric morphometric measurements of shell shape. Our objectives were to (a) identify the location of dispersal barriers based on two genetic markers, (b) test for congruence of genetic barriers with differences in shell shape, and (c) assess whether the spatial structure was better explained by shifts in abundance or oceanographic parameters. By identifying dispersal barriers and the possible drivers of population structure in holoplanktonic gastropods, we can better predict their capacity to respond to a rapidly acidifying ocean.

| Species ecology and biology
Limacina bulimoides is a sinistrally coiled holoplanktonic gastropod, well-adapted to life in the open ocean. It has a small (maximum shell length = 2 mm, van der Spoel et al., 1997), thin, aragonitic shell and two parapodia, which it uses to swim through the water column in analogous fashion to how insects fly through the air-a remarkable example of convergent evolution (Murphy et al., 2016).
For feeding, an external mucous web traps particulate matter, including phytoplankton and small protists, which is brought to the mouth by ciliary movement (Lalli & Gilmer, 1989). The radula, a typical molluscan structure used for feeding, has been described from a congener, Limacina helicina, and was reduced, with 10 rows of teeth (Gilmer & Harbison, 1986;Lalli & Gilmer, 1989;Meisenheimer, 1905;Ritcher, 1979). However, in a recent detailed morphological study of the same species, only six rows of teeth were found, which led the authors to suggest that the number of rows could be variable due to continuous regrowth (Laibl et al., 2019). Limacina bulimoides lay free-floating egg strings (see Figure S1/Video S1), from which freeswimming veliger larvae hatch (Lalli & Wells, 1978). Limacina spp.
The generation time of L. bulimoides is estimated to be less than a year, with larvae metamorphosing into juveniles after 2 months, juveniles reaching sexual maturity as males after a subsequent 5 months, and reaching their maximum size as females 2 months later (Wells, 1976). Reciprocal copulation has been observed between males (Morton, 1954), or between individuals that are in transition between male and female (Lalli & Gilmer, 1989). They possess a penis for internal fertilization (Lalli & Wells, 1978), and it is unknown whether self-fertilization is possible.  (Table 1). Additionally, one sample site of L. bulimoides from the Pacific was included as an outgroup to provide a broader perspective on the diversity found for the Atlantic individuals (Table 1, Figure 1). For population-level analyses, the Longhurst province assignments were simplified to four major ocean provinces (North Gyre, Equatorial, South Gyre and Convergence) (Table 1), which also take into account previously reported transitions in planktonic community composition along the AMT transects Goetze et al., 2017).

| Sampling
Abundance of L. bulimoides was obtained from AMT24 (September-October 2014) as quantitative samples from AMT22 were not available. For this equivalent AMT transect, bulk plankton samples were collected using a bongo net (200 µm mesh size) with a General Oceanics flowmeter (2030RC) mounted in the mouth to quantify the volume of seawater filtered (see Burridge, Hörnlein, et al., 2017). Based on the peaks in abundance of L. bulimoides across the AMT24 transect and comparison of oceanographic measurements from both the AMT22 and AMT24 transects, the sampled stations were split into three 'population groups' (North, Equatorial, South) (Table 1).

| Molecular analyses
To resolve spatial genetic structuring of L. bulimoides, 356 individuals were sequenced for a 565 base pair (bp) fragment of the mitochondrial COI gene, and 344 individuals for a 938 bp fragment of the nuclear 28S gene (Table 1). In total, 332 L. bulimoides specimens were sequenced for both COI and 28S (Table S1). The aver- Haplotype (H d ) and nucleotide diversity (π) were calculated using DnaSP v6.12.01 (Rozas et al., 2017) and reported for each station for both COI and 28S fragments. The 28S data were phased in DnaSP using PHASE (Stephens & Donnelly, 2003)  Lischer, 2010). Minimum-spanning networks for both genes were calculated and visualized in POPART (Leigh & Bryant, 2015). The phased allele alignment was used to construct the 28S haplotype network.
To identify population subdivision, genetic diversity between populations from different stations was quantified by pairwise Φ ST (Holsinger & Weir, 2009), using the pairwise difference method in Arlequin with 10,000 permutations, and significant comparisons were identified after strict Bonferroni correction (Rice, 1989).
Hierarchical Analysis of Molecular Variance (AMOVA) tests were conducted for both genes separately with sampling sites partitioned into (a) four 'ocean province' groups and (b) three abundance-based 'population' groups (Table 1). To test for isolation-by-distance (IBD), Mantel tests with 10,000 permutations were conducted in Arlequin to assess correlation between the COI or 28S pairwise Φ ST matrices and the geographic distance matrix. Pairwise geographic distance between stations was calculated using the R package geodist with geodesic measures (Padgham & Sumner, 2020). The distance between stations 62 and 64 was averaged since they were analysed as a single site. The pairwise scatter plots were produced in R with ggplot2 (Wickham, 2016).

| Morphometric analyses
Variation in shell shape was assessed using geometric morphometric measurements of 136 undamaged adult shells (shell length > 0.9 mm).
Specimens were positioned in a standardized apertural orientation and photographed using a Zeiss V20 stacking stereomicroscope with Axiovision software, prior to destructive DNA extraction.
The resulting images were processed and digitized at eight (semi-) landmarks in TpsUtil and TpsDig (Rohlf, 2015) ( Figure S2). The coordinates of the (semi-)landmarks were analysed in TpsRelW (Rohlf, 2015), using a generalized least-squares Procrustes superimposition (Rohlf & Slice, 1990;Zelditch et al., 2004). To ensure that the digitization process was repeatable, a subset of 30 individuals was photographed and digitized twice, with eight landmarks placed on each image. Centroid size and relative warp (RW) scores between the pairs of images per specimen after Procrustes Fit were compared using intra-class coefficient (ICC) in PAST3.0 (Hammer et al., 2001), and ICC values > 0.75 were considered sufficiently repeatable.
We tested for significant variation in shell shape across genetically distinct groups with a nonparametric Permutational Multivariate Analysis of Variance (one-way PERMANOVA) using Euclidean distances and 9,999 permutations in R with vegan (Oksanen et al., 2019). Only centroid size and repeatable RWs were used in the one-way PERMANOVA, and strict Bonferroni corrections were applied for the pairwise PERMANOVAs. The first two RW axes were plotted to visualize shell shape variation for different groups of L. bulimoides. Additionally, a Canonical Variates Analysis (CVA) was conducted in R (R Core Team, 2017) to discriminate shell morphometric differences between groups. A one-way ANOVA with a post hoc Tukey HSD test was also conducted in R to test if the means of the canonical variate for each group were different from the other groups.

| Spatial structure
Overall Φ ST calculated for both genes across the Atlantic supports the inference of significant spatial population structure Pairwise Φ ST values for COI ranged from −0.017 to 0.451 for comparisons along the Atlantic transect and from 0.724 to 0.783 for comparisons between Atlantic and Pacific samples (Table 3) Table 3).   (Figure 3a,b). Surprisingly, this barrier is located in the core of the southern oligotrophic gyre and appears incongruent with any strong oceanographic gradients (Figure 3c,d). A more subtle secondary barrier is located in the region 15°N-4°S, which corresponds to the equatorial upwelling province and is mainly supported by 28S data.
We did not find obvious patterns of IBD across the Atlantic transect ( Figure S3). Although there was a significant correlation between the pairwise Φ ST and geographic distance for COI along the complete Atlantic transect (r = 0.520, p = 0.00120, Table S2 Table S2). Pairwise Φ ST values were uniformly high regardless of geographic distance when comparing stations on opposite sides of this dispersal barrier, and uniformly low when comparing stations within population groups ( Figure S3a). This suggests that the F I G U R E 3 Environmental profile of two oceanic dispersal barriers for Limacina bulimoides in the Atlantic Ocean. The primary dispersal barrier is located at 15-18°S, while the secondary barrier is located in the region between 15°N and 4°S. (a) Abundance of L. bulimoides per 1,000 m 3 seawater along the meridional transect sampled in 2014 (data from Burridge, Hörnlein, et al., 2017). Population groups are coloured as in Figure 2, supported by 28S genetic structuring (see Table 3 The two mitochondrial haplogroups appear diagnostic for individuals on either side of the primary dispersal barrier with most individuals from the North and Equatorial stations belonging to COI haplogroup 1 and most individuals from the South belonging to COI haplogroup 2 (Figure 3b). However, we found 13 exceptions that are probably expatriates, or individuals sampled outside of their core range. Seven expatriates from the South clustered within haplogroup 1, and six expatriates from the North + Equatorial groups clustered within haplogroup 2 (Table S3; Figure 3b). We could verify the life stage for nine of these expatriates, of which three were adults (shell length > 0.9 mm) and six were juveniles (Table S3). We assessed whether the 28S and mitochondrial haplotypes were congruent for 12 of these expatriates, but these results were inconclusive due to a lack of variability in the 28S fragment ( Figure S4a).

| Phenotypic variability
The centroid size and RW axes 1, 2, 4, 6, and 7 were repeatable and included in the final analyses. Based on the complete dataset of adult L. bulimoides, including Pacific individuals (N = 136), the repeatable RWs accounted for 85.21% of the total shell shape variation (Table S4). Shell shape variation for all specimens (N = 136) and Atlantic specimens only (N = 111) was visualized on the first two axes, which explain the most variation: RW1 (All: 61.08%; Atlantic: 62.41%) and RW2 (All: 13.45%; Atlantic: 13.43%) ( Figure 4). The thin-plate spline reconstructions show that a positive RW1 is associated with a taller spire, while a negative RW1 is associated with a shorter, more rounded spire. Along RW2, larger values indicate wider shells.
We did not observe a significant difference in shell shape be-  (Figure 4b). Similarly, there is a significant difference in shell shape across all three groups (F(2, 108) = 29.79, p = 0.00000) based on the CVA. Significant shell shape differences were recorded for both North-South (p = 0.00000) and Equatorial-South (p = 0.00003) comparisons, but not between North and Equatorial groups (p = 0.889) based on Tukey HSD multiple pairwise comparisons ( Figure S5, Table S5). Thus, the spatial structuring of samples in terms of shell shape matches the partitioning based on mitochondrial genetic variation (Figure 3b).

TA B L E 4
Results of hierarchical AMOVA for mitochondrial cytochrome c oxidase I (COI) (n = 334) and nuclear 28S (28S) (n = 324, 2n = 648) sequence variants in Atlantic individuals of Limacina bulimoides, partitioned based on four ocean provinces or three abundancebased population groups (Table 1

| D ISCUSS I ON
We found high levels of genetic variability and significant spatial population structure in the pteropod L. bulimoides along a latitudinal transect (35°N-36°S) in the Atlantic Ocean. Specifically, we have evidence for a primary dispersal barrier located in the southern subtropical gyre at 15-18°S and a more subtle, secondary barrier located in the region 15°N-4°S. This finding is in stark contrast to the historical perception of the pelagic environment as a homogeneous habitat lacking strong isolating barriers, and of holoplanktonic organisms as characterized by high dispersal and broad biogeographic ranges (Norris, 2000;Pierrot-Bults & van der Spoel, 1979). Because shelled pteropods are regarded as sentinel species to assess the impacts of ocean acidification , it is important to document their genetic variation and population structure, as these characteristics will influence their response to natural selection and global change (Barrett & Schluter, 2008;Donelson et al., 2019;Riebesell & Gattuso, 2015;Sanford & Kelly, 2011). This is the first study examining genetic structuring of a shelled pteropod with rigorous sampling at an ocean basin scale. The presence of strong oceanic dispersal barriers may indicate that the potential for range shifts and adaptive responses to a changing ocean could be more limited than expected, and that local adaptation should also be taken into account, for instance when modelling (future) species distribution patterns. Here, we discuss the nature and possible drivers of the oceanic dispersal barriers reported here for L. bulimoides, and more generally for pelagic taxa.
The primary dispersal barrier in the Atlantic Ocean probably represents a species barrier, across which populations on either side are reproductively isolated. This inference is supported by congruent differentiation at independently inherited genetic loci (COI and 28S) and congruent genetic and shell shape differences.
Despite the significant IBD test for the COI gene across the Atlantic transect, the disappearance of this result when populations across the primary barrier were tested separately supports the inference that these populations are spatially discrete evolutionary lineages (Teske et al., 2018). Nevertheless, a direct test of reproductive isolation is not possible as the high mortality of pteropods under laboratory conditions precludes meaningful crossing experiments (Howes et al., 2014). Genome-wide data would be necessary to determine the strength of reproductive isolation by quantifying divergence at many loci across the genome (Butlin & Stankowski, 2020; Gagnaire, 2020).
The observed morphometric variation was associated with the primary genetic barrier and not with any strong environmental transitions. Hence, it is more likely that the shell shape differences are the result of accumulated genetic differences than phenotypic plasticity resulting from distinct environments. Examples of strong plasticity in shells of marine molluscs are generally found across heterogeneous habitats, such as L. saxatilis with discrete 'crab' and 'wave' ecotypes (e.g. Hollander & Butlin, 2010;Hollander et al., 2006) or Nacella concinna with continuous morphological differentiation across a depth gradient (Hoffman et al., 2010). Shell shape variation has also been shown to have a large genetic component, with heritabilities ranging between 0.36 and 0.71 in L. saxatilis The secondary dispersal barrier across the equatorial upwelling province is supported by limited evidence from the 28S gene (Table 3), and is congruent with known barriers in other pelagic taxa.
The equatorial province has been suggested as an ecological barrier in some nekton species, such as Atlantic bluefin tuna (Briscoe et al., 2017) and whales (Holt et al., 2020), and as a barrier separating anti-tropical clades in planktonic foraminifers (Casteleyn et al., 2010). Dispersal barriers across the equatorial upwelling region have also been observed for other zooplankton. For example, in the pteropod genus Cuvierina, two distinct morphotypes are found in the Atlantic basin. One morphotype, 'atlantica', has disjunct populations in the northern and southern subtropical gyres, which are genetically differentiated, while the morphotype 'cancapae' dominates in equatorial waters (Burridge et al., 2015). Similar patterns of equatorial dispersal barriers were observed in the epipelagic gastropods Protatlanta sculpta and P. souleyeti (Wall-Palmer et al., 2016) and the mesopelagic copepod H. longicornis (Andrews et al., 2014;Norton & Goetze, 2013). We additionally observed a pattern of unique diversity in the equatorial region for the 28S gene fragment. This trend of higher equatorial diversity or equatorial endemics was also observed in populations of the copepods P. abdominalis (Hirai et al., 2015) and P. xiphias . It was suggested that this unique diversity resulted from a resident equatorial population, which was supplemented with expatriates by advection from the neighbouring gyre populations . It is as yet unknown if expatriates in L. bulimoides survive to reproduce with the resident population and contribute to gene flow, or whether expatriates represent evolutionary dead-ends due to a mismatch of phenotype and environment (Marshall et al., 2010). Increased sampling intensity in the equatorial upwelling region, by having plankton tows with a greater amount of seawater filtered, would allow for the collection of more individuals in this region of low abundance for L. bulimoides, and a more thorough investigation into the prevalence and genetic identity of expatriates.
We identified two dispersal barriers for L. bulimoides in the Atlantic Ocean, of which one (the secondary barrier) appears more permeable than the other (the primary barrier). The locations of these barriers coincide with regions of low abundance, congruent with the hypothesis that areas of less suitable habitat contribute to the maintenance of population structure in oceanic plankton Norris, 2000;Peijnenburg et al., 2004). The presence of expatriates suggests that Atlantic L. bulimoides populations could meet and potentially interbreed. This raises the question: what are the possible isolating mechanisms? These may include overlooked physical boundaries or ecological differences, including their depth habitat, feeding strategy or reproductive biology. Interestingly, the primary dispersal barrier is congruent with complex mesopelagic circulation patterns in the South Atlantic (~20°S, Sutton et al., 2017) where a shift in shrimp species (Acanthephyra kingsleyi to Acanthephyra quadrispinosa) was observed at 18°S (Judkins, 2014). Little is known, however, about the behaviour or depth habitat of L. bulimoides across the Atlantic transect and thus it is possible that mesopelagic currents are important.
Since internal fertilization is the reproductive mode in L. bulimoides, assortative mate choice could also be a possible reproductive isolating mechanism, with chemical cues facilitating mate location and communication, as observed in a variety of marine gastropods (Ng et al., 2013). Other isolating mechanisms could involve differences in reproductive timing (McClary & Barker, 1998), ecological selection and community interactions (Whittaker & Rynearson, 2017), or post-zygotic barriers, such as reduced hybrid fitness (Lessios, 2007) or reduced hatching success of offspring (Ellingson & Krug, 2015).  (Hellberg & Vacquier, 1999), the divergence between the Atlantic and Pacific mitochondrial clades is estimated at ~4.6 MY, which roughly coincides with the emergence of the Isthmus of Panama (Bacon et al., 2015;O'Dea et al., 2016). The rise of the Isthmus of Panama has been directly linked to allopatric divergence across various marine taxa, due to the shoaling of the seaway, which reduced connectivity between mesopelagic populations at about 7 MY. This was followed by a decrease in seasonal upwelling due to the closing seaway at 4 MY, which cut off connectivity between plankton populations and decreased the input of primary productivity into the Caribbean (Knowlton et al., 1993;Leigh et al., 2014). This major marine biogeographic barrier was also used to calibrate the divergence between Atlantic and Pacific populations of three other pteropod species in a phylogenetic analysis (Burridge, Hörnlein, et al., 2017).
The population structure within the Atlantic basin could be the result of different, not mutually exclusive, evolutionary scenarios.
Dispersal barriers could have resulted from secondary contact following allopatric divergence in the past. The divergence time of the two Atlantic haplogroups across the primary barrier is estimated at 1.10 million years (MY) using a strict COI molecular clock (Hellberg & Vacquier, 1999). This estimate points to the potential effect of Pleistocene glacial cycles on population divergence of L. bulimoides. The effect of Pleistocene cycles as important drivers of phylogeographic structure has also been proposed for many other marine molluscs (e.g. Luttikhuizen et al., 2003;Marko, 2004;Reid et al., 2006), including polar Limacina species (Sromek et al., 2015). Another possible scenario is parapatric divergence, where semi-isolated populations living along an environmental cline diverged due to ecological selection. This would be a more plausible scenario for the secondary dispersal barrier across the equatorial upwelling gradients. The populations would progressively become more distinct, when the effect of divergent selection is greater than that of homogenizing gene flow (Crow et al., 2007;Luttikhuizen et al., 2003;Postel et al., 2020;Schluter, 2009). In order to differentiate between historical divergence followed by secondary contact, ecological selection or a combination of the two, genome-wide analyses are required to investigate if genetic differentiation is randomly distributed across the genome (suggesting a gradual process) or associated with particular genes or environmental clines (indicative of ecological selection).
Holoplanktonic organisms are useful to gain insight into global biogeographical patterns in pelagic dispersal (Álvarez-Noriega et al., 2020;Bradbury et al., 2008), as their dispersal potential is not confounded by settlement processes that involve complex biophysical interactions with coastal or benthic habitats (Pineda et al., 2009;Prairie et al., 2012;Weersing & Toonen, 2009). To gain further insight into the nature and drivers of dispersal barriers in the open ocean, as well as to identify signals of selection across the genome, a broader range of ecological observations and in-depth analyses of genome-wide diversity in zooplankton will be necessary (Bucklin et al., 2018;Choo et al., 2020;Choquet et al., 2019;Gagnaire et al., 2015). This is important because marine zooplankton are key players in pelagic food webs and useful indicators as rapid responders to environmental variation and climate change (Beaugrand et al., 2009). Shelled pteropods could be particularly suitable model organisms to gain insight into the evolutionary potential of marine zooplankton and to deepen our understanding of speciation processes in the open ocean for a number of reasons. First, there is growing interest in the group because of their vulnerability to ocean acidification and their potential use as bioindicators. Second, unlike most holoplanktonic animals, they have shells, which provide a record of the ecology and life history of an individual that can be easily measured. Third, pteropods are the only living planktonic animals with a good fossil record Peijnenburg et al., 2020), which allows for comparisons of morphological diversity in extant as well as extinct species and populations. We demonstrated significant spatial population structuring across and within ocean basins in L. bulimoides, a presumably circumglobal species. This underlying genetic variability and structure will influence the response of populations to global change, and in-depth study of their life-history traits, ecology and historical demography, in combination with experimental exposure to future conditions, will be required to more accurately predict their ability to adapt to a rapidly changing ocean. and is contribution number 352 of the AMT programme.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13735.