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.

Abstract

Shelled pteropods or 'sea butterflies' are proposed as bio‐indicators to monitor impacts of ocean acidification on open ocean ecosystems. Local adaptation is affected by the presence and geographic extent of distinct populations, however, little is known about the large‐scale population structure of pteropods. We thus examined the spatial distribution of genetic and morphometric variation in the species Limacina bulimoides across the Atlantic Ocean. Two dispersal barriers were identified that were congruent with regions of low abundance, suggesting that regions of suboptimal habitat may create barriers to gene flow in widely‐distributed zooplankton species.

Open in new tabDownload slide

INTRODUCTION

Anthropogenic carbon emissions are an important cause of perturbations in marine ecosystems, including an overall increase in seawater temperature, ocean deoxygenation and acidification (Gattuso et al., 2015; 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, 2014; Hays 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; Bednaršek et al., 2012; 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 (Bednaršek, Tarling, et al., 2012; Lischka et al., 2011; Manno et al., 2017), and they have been suggested as bioindicators of the effects of ocean acidification (Bednaršek, Klinger, et al., 2017). 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, 2019; Shimizu 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 (Bednaršek, Feely, et al., 2017; 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 (Goetze et al., 2017) 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 (Goetze et al., 2015, 2017). 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; Burridge, Goetze, et al., 2017). 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.

MATERIALS AND METHODS

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 free‐swimming veliger larvae hatch (Lalli & Wells, 1978). Limacina spp. are protandrous hermaphrodites, hence, the larvae sexually mature into males before transitioning into females (Lalli & Gilmer, 1989). 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.

Sampling

Bulk plankton samples were collected on AMT cruise 22 (AMT22) in October and November 2012 (Table 1, Figure 1). The samples were collected at night by oblique tows of a bongo net (200 and 333 µm mesh sizes) or RMT1 midwater trawl (333 µm) from a median depth of 323 m to the sea surface (depth range: 132–402 m). Samples were preserved in 95% ethanol, stored at −20°C, and subsequently sorted in the laboratory. Seawater temperature and chlorophyll a concentration in the upper 300 m of the water column were obtained using a Sea‐Bird Electronics 3P Temperature Sensor and Chelsea MKIII Aquatracka Fluorometer, with data calibrated and archived by the British Oceanographic Data Centre (BODC, https://www.bodc.ac.uk). The 14 samples included in this study were collected between 35°N and 36°S in the following Longhurst biogeochemical provinces (Longhurst, 2007; Reygondeau et al., 2013): North Atlantic Subtropical gyre (NAST), North Atlantic Tropical gyre (NATR), Western Tropical Atlantic (WTRA), South Atlantic gyre (SATL), South Subtropical Convergence (SSTL) (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 (Burridge, Goetze, et al., 2017; Burridge et al., 2017; Goetze et al., 2017).

FIGURE 1

Sampling locations of Limacina bulimoides along the basin‐scale Atlantic Meridional Transect (AMT22, left), as well as in the Pacific Ocean off Hawaiʻi (right). Numbers indicate the station number for each site, and symbols indicate the type of data that were obtained (see legend)

Open in new tabDownload slide

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 average number of individuals sampled per station for COI and 28S was 24 and 23, respectively. Individuals from stations 62 and 64 were combined because of their geographic proximity and low numbers of specimens. DNA was extracted with either the NucleoMag 96 Tissue Kit (Macherey‐Nagel GmbH & Co. KG) or the DNeasy 96 Blood & Tissue Kit (Qiagen). DNA extracts were amplified in polymerase chain reaction (PCR) using the primers jgLCO1490 (5′‐TITCIACIAAYCAYAARGAYATTGG‐3′) and jgHCO2198 (5′‐TAIACYTCIGGRTGICCRAARAAYCA‐3′) for COI (Geller et al., 2013), and for 28S using C1‐F (5′‐ACCCGCTGAATTTAAGCAT‐3′) (Dayrat et al., 2001) and D3‐R (5′‐GACGATTCGATTTGCACGTCA‐3′) (Vonnemann et al., 2005). PCR was carried out with a reaction mix consisting of 1 μl template DNA, 17.8 μl milli‐Q H₂O (Ultrapure), 2.5 μl PCR buffer CL (10X) (Qiagen), 0.5 μl MgCL₂ (25 mM) (Qiagen), 0.5 μl BSA (100 mM) (Promega), 1 μl of each primer (10 μM), 0.5 μl dNTPs (2.5 mM) and 0.25 μl Taq polymerase (5 U/μl) (Qiagen). The initial denaturation step of 3 min at 96°C was followed by 40 cycles of 15 s at 96°C, 30 s at 50°C and 40 s at 72°C, ending with a final extension step of 5 min at 72°C. Amplified products were sequenced in both directions at Baseclear B.V. (Leiden, the Netherlands), combined and checked for errors in Geneious v8.1.9 (Kearse et al., 2012) and the final sequences were aligned using MAFFT v7.017 (Katoh & Standley, 2013).

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) with default settings of 100 iterations, thinning interval of 1 and burn‐in of 100. Stations were tested for adherence to neutrality assumptions, with Tajima's D (Tajima, 1989) calculated in Arlequin version 3.5.2.2 (Excoffier & 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.

RESULTS

Genetic variability

The COI locus was highly polymorphic with global haplotype diversity (H_d) of 0.9998 (range: 0.9820–1.000) and nucleotide diversity (π) of 0.05449 (range: 0.0283–0.0402) across all sites (Table 2). In the COI minimum‐spanning network, 349 unique haplotypes were observed across 356 individuals, and there was no shared central haplotype (Figure 2). Instead, Atlantic individuals were clustered into two main haplogroups that were separated by 15 substitutions (haplogroups 1 & 2). A third cluster was comprised of all Pacific individuals, which was separated from Atlantic haplogroup 1 by 62 substitutions. Population samples from stations 13, 19, 25, 43 and 45 showed deviations from neutrality based on Tajima's D (Table 2).

FIGURE 2

Minimum‐spanning networks of Limacina bulimoides gene fragments of mitochondrial cytochrome c oxidase I (COI, top) and nuclear 28S rDNA (28S, bottom). The size of the filled circles represents the number of individuals with each haplotype (COI, 28S), with the smallest circles representing one individual with that haplotype, while colour represents abundance‐based population groups (see map insert and Table 1). For the COI network, the dashed lines around the haplotype clusters indicate the three main haplogroups. Hatch marks on the branches represent the number of mutational steps, with large numbers of mutations between the haplogroups indicated with numbers

Open in new tabDownload slide

The 28S gene was less polymorphic than COI, with a global H_d of 0.6740 (range: 0.3740–0.8550) and a π of 0.0027 (range: 0.0006–0.0035) across all sites (Table 2). This reduced polymorphism is visible in the 28S minimum‐spanning network (Figure 2), with only 63 haplotypes observed across 344 individuals. The three most abundant haplotypes comprised 369, 129 and 45 copies, respectively, and included representatives from both Atlantic and Pacific samples. The second most abundant haplotype was mainly comprised of individuals sampled from the equatorial upwelling region. Only station 55 showed a deviation from neutrality based on Tajima's D (Table 2).

Spatial structure

Overall Φ_ST calculated for both genes across the Atlantic supports the inference of significant spatial population structure (COI Φ_ST = 0.230, p = 0.00000; 28S Φ_ST = 0.255 p = 0.00000). 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). All comparisons between stations 13–49 (North & Equatorial) and stations 51–64 (South) were significant, with high Φ_ST values ranging from 0.250 to 0.451. We did not find any significant Φ_ST values in pairwise comparisons between stations 13–49. For nuclear 28S, pairwise Φ_ST values ranged from −0.001 to 0.621 for comparisons along the Atlantic transect, and from 0.079 to 0.532 for comparisons between Atlantic and Pacific samples (Table 3). The 28S pairwise Φ_ST values show a similar genetic break as for COI, separating stations 13–49 (North & Equatorial) from 51–64 (South), as 35 of the 40 pairwise Φ_ST values were significant (Φ_ST ranging from 0.012 to 0.621). There is additional spatial structuring separating samples from stations 13–29 (North) and 43–49 (Equatorial), as 13 of 17 pairwise Φ_ST comparisons were significant (Φ_ST ranging from −0.001 to 0.497, Table 3).

Hierarchical analyses of molecular variance show that the spatial structure based on the trimodal pattern of abundance of L. bulimoides across a similar Atlantic transect (AMT24, Figure 3a) explains the distribution of genetic variation better than partitioning based on oceanographic parameters (Table 4). For both COI and 28S fragments, the fixation index among groups was higher when partitioned among three abundance‐based 'population groups' (COI: Φ_CT = 0.285, p = 0.00129; 28S: Φ_CT = 0.280, p = 0.00010) than among four 'ocean provinces' (COI: Φ_CT = 0.205, p = 0.00446; 28S: Φ_CT = 0.140, p = 0.0620). Thus, we find a strong primary dispersal barrier separating stations 13–49 (North & Equatorial) from stations 51–64 (South) located at 15–18°S that is supported by both genes (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.

FIGURE 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³ 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) along the Atlantic Meridional Transect 22 (AMT22) sampled in 2012, locations indicated by grey asterisks. (b) Piecharts showing mitochondrial cytochrome c oxidase I (COI) haplogroup frequencies per sampling site with the division into two morphogroups superimposed (see also Figure 4b). (c) Ocean section plots of chlorophyll a concentration and (d) seawater temperature from the upper 300 metres of the water column with latitude. Sampling locations are marked by vertical lines, and major oceanographic provinces (Convergence, South Gyre, Equatorial, and North Gyre) are indicated. Plots b, c and d derive from material sampled in 2012 (AMT22 transect)

Open in new tabDownload slide

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), this relationship breaks down when population groups on either side of the primary barrier at 15–18°S were analysed separately (North + Equatorial: r = 0.117, p = 0.269; South: r = 0.0968, p = 0.418, 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 significant correlation along the entire transect can be explained by the inclusion of high pairwise Φ_ST values between the stations across the barrier. For 28S, correlation of genetic and geographic distance matrices was not observed for the entire Atlantic transect (r = −0.0942, p = 0.770), but only for the North + Equatorial group (r = 0.759, p = 0.003, Table S2; Figure S3b). Again, this correlation is likely driven by the high Φ_ST values for comparisons between North and Equatorial samples on either side of the secondary dispersal barrier. No IBD patterns were found when any of the population groups were analysed separately.

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.

FIGURE 4

Shell shape variation in adult Limacina bulimoides visualized for the first two Relative Warp axes for (a) Atlantic and Pacific individuals (N = 136), and (b) Atlantic individuals (N = 111) coloured according to abundance‐based population groups. Thin‐plate splines of the most positive and negative deformations along the axes are indicated to depict the variation in shell shape. Photographs of typical individuals from the Atlantic and Pacific (a) and Atlantic North and South population groups (b) are also shown

Open in new tabDownload slide

We did not observe a significant difference in shell shape between Atlantic and Pacific L. bulimoides individuals (N = 136, F (1,135) = 1.66, p = 0.204). Shell shapes of the Pacific individuals (N = 25) fall within the variation observed for the Atlantic individuals (N = 111, Figure 4a). When Atlantic individuals were classified according to their spatial structure with North, Equatorial and South population groups, we found a significant difference in shell shape between groups based on the one‐way PERMANOVA (F (2,110) = 14.53, p = 0.0001). Shell shape was different between the North and South (F = 35.9, p = 0.003), whereas comparisons of shell shape between the North and Equatorial or South and Equatorial groups were not significant (F = 2.74, p = 0.113 and F = 2.67, p = 0.116, respectively). This result could be due to the low sample size of the Equatorial group (N = 9). Individuals from the North population group generally had a higher, narrower spire compared to individuals from the South (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).

DISCUSSION

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 (Bednaršek, Klinger, et al., 2017), 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 (Conde‐Padín et al., 2007) and between 0.16 and 0.56 in Nucella lapillus (Guerra‐Varela et al., 2009). Furthermore, shell shape differences in conjunction with genetic differentiation have been used to identify species boundaries in different genera of holoplanktonic gastropods (Burridge et al., 2015, 2019; Wall‐Palmer et al., 2018) using an integrative species concept (McManus & Katz, 2009; Padial et al., 2010). Hence, strong reproductive isolation is probably involved in the maintenance of this primary dispersal barrier in the southern Atlantic gyre.

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 (Goetze et al., 2017). 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 (Goetze et al., 2017). 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 (Goetze et al., 2017; 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). A broader range of ecological observations are required to determine the exact nature of dispersal barriers for L. bulimoides. For instance, future research to reveal population differences in feeding ecology or depth habitat could include the metabarcoding of gut contents, examination of the radula, or stable isotope analyses. Because we observe differences in the locations of maximum abundance and genetic transitions in different planktonic taxa sampled across similar Atlantic transects (Andrews et al., 2014; Goetze et al., 2015, 2017; Hirai et al., 2015; Norton & Goetze, 2013), we expect that dispersal barriers are species‐specific and largely driven by the ecological characteristics of the respective species.

Many circumglobally distributed marine species have distinct populations or (sub)species in different ocean basins. These subdivisions generally result from divergence in allopatry, that is separated by current systems or continental land masses, and are often congruent with global biogeochemical provinces (Bowen et al., 2016; Costello et al., 2017; Reygondeau et al., 2012). The significant genetic divergence between Atlantic and Pacific L. bulimoides is in agreement with the patterns of divergence observed between ocean basins in other circumglobally distributed marine taxa, including fish (e.g. Graves & McDowell, 2015), benthic organisms (e.g. Fauvelot et al., 2020; Prazeres et al., 2020) and other zooplankton (reviewed in Peijnenburg & Goetze, 2013). Based on a strict COI molecular clock rate of 2.4%/million years (MY) for gastropods (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 (Janssen & Peijnenburg, 2017; 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.

ACKNOWLEDGMENTS

This research was supported by the Netherlands Organisation for Scientific Research (NWO) Vidi grant 016.161.351 to K.T.C.A.P. Fieldwork also was supported by NSF grants OCE‐1029478 and OCE‐1338959 to E.G. We thank A. Burridge and M. Jungbluth for assistance at sea, and G. Tarran, A. Rees, and the captains and crew of the Atlantic Meridional Transect (AMT) cruises 22 and 24 for their leadership and assistance. We thank A. Tsuda for providing Pacific samples from the 2011/2012 R/V Hakuho‐Maru KH‐11‐10 cruise. We are also grateful towards two anonymous reviewers whose comments have greatly improved the manuscript. The Atlantic Meridional Transect programme is funded by the UK Natural Environment Research Council through its National Capability Long‐term Single Centre Science Programme, Climate Linked Atlantic Sector Science (grant number NE/R015953/1). This study contributes to the international IMBeR project and is contribution number 352 of the AMT programme.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

Peer Review

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13735

DATA AVAILABILITY STATEMENT

The nucleotide sequences for the COI and 28S genes are archived on NCBI GenBank with the accession numbers MN952611–MN952965 and MN950433–MN950775, respectively. Images for morphometric analyses are archived on Dryad (https://doi.org/10.5061/dryad.p8cz8w9nq).

REFERENCES