Assessing species boundaries in the open sea: an integrative taxonomic approach to the pteropod genus Diacavolinia

Disclaimer/Complaints regulations If you believe that digital publication of certain material infringes any of your rights or (privacy) interests, please let the Library know, stating your reasons. In case of a legitimate complaint, the Library will make the material inaccessible and/or remove it from the website. Please Ask the Library: https://uba.uva.nl/en/contact, or a letter to: Library of the University of Amsterdam, Secretariat, Singel 425, 1012 WP Amsterdam, The Netherlands. You will be contacted as soon as possible.


INTRODUCTION
Integrative taxonomy aims to rigorously delimit species and prevent under-or overestimation of species numbers by statistically testing species hypotheses with diverse character and data types (Edwards & Knowles, 2014;Karanovic et al., 2016). An increasing number of species is described each year as a result of the incorporation of new tools for species discovery, including virtual access to museum collections, advances in DNA sequencing, morphometric methods and geographic information systems (Vogler & Monaghan, 2006;Knapp, 2008). Such tools have enabled integrative taxonomic approaches, in which species are described based on congruence between morphological and genetic information, with additional supporting characteristics such as behaviour, ecology or geography (McManus & Katz, 2009;Padial et al., 2010;Smith & Hendricks, 2013;Edwards & Knowles, 2014). Statistical identification of genetic lineages plays a pivotal role in species detection and satisfies multiple species concepts, because it treats species as hypotheses using objective tests (De Queiroz, 2007;Hausdorf, 2011;Morard et al., 2016). However, the sole use of genetic information in species delimitation may fail to detect the same number of species as methods that combine multiple lines of evidence (Edwards & Knowles, 2014;Burridge et al., 2015;Karanovic et al., 2016). Species may go undetected based on a limited set of selected genetic markers, because they may be distinct only in other genes, morphology or ecological niche space, or species numbers may be overestimated by the use of highly variable genetic markers. Geometric morphometric methods, in which quantitative differences in organismal shape and size are used to distinguish between taxa without requiring analysis of soft tissue, may especially strengthen studies based on limited datasets, such as those containing fossil taxa and valuable museum specimens for which genetic information cannot be obtained (Karanovic et al., 2016). Developing integrative taxonomic frameworks for assessing species boundaries enables the inclusion of museum specimens originally used to describe species in studies that also incorporate fresh samples, given that there are enough specimens available for connecting different data types.
Tracking changes in marine biodiversity in response to climate change on a global scale requires accurate assessment of species boundaries and distributions (e.g. Goetze & Ohman, 2010;Churchill et al., 2014a, b;Burridge et al., 2015;Wall-Palmer et al., 2018). Range shifts of planktonic taxa have been among the most rapid, have occurred over the largest spatial scales in comparison to other marine and terrestrial groups, and may affect higher trophic levels in the marine food web (Richardson, 2008;Beaugrand et al., 2012Beaugrand et al., , 2014Beaugrand et al., , 2015Parmesan et al., 2013;Brown et al., 2016). Plankton distributions are often concordant with biogeochemical provinces at the level of species and communities (e.g. Valentin & Monteiro-Ribas, 1993;Longhurst, 1998;Reygondeau et al., 2013;Dolan et al., 2016;Burridge et al., 2017a,b), as well as at the level of population genetic structure within species (e.g. Norton & Goetze, 2013;Goetze et al., 2015Goetze et al., , 2017Hirai et al., 2015). Persistent dispersal barriers may limit range shifts of some taxa in response to changing ocean conditions, while other taxa may be able to adapt and occupy new ecological niches. Integrative taxonomy will improve the accuracy of marine species delimitation, enable the identification of rare taxa and provide insights in current species distributions, with the potential to predict future range shifts.
P t e r o p o d s a r e a g r o u p o f h o l o p l a n k t o n i c gastropods that has been identified as exceptionally vulnerable to ocean acidification (Fabry et al., 2008;Bednaršek et al., 2016;Manno et al., 2017). Pteropods comprise the Thecosomata (Euthecosomata and Pseudothecosomata), shelled or semi-shelled pteropods commonly known as 'sea butterflies' and the Gymnosomata, shell-less pteropods known as 'sea angels'. Thecosome pteropods have aragonitic shells, and are a diverse group of organisms that are common in marine zooplankton from polar to equatorial habitats (Lalli & Gilmer, 1989). Shelled pteropods have an extensive fossil record (Janssen, 2007a(Janssen, , b, 2012Janssen & Peijnenburg, 2017), and are commonly used to examine the effects of ocean acidification on marine life (Roger et al., 2011;Comeau et al., 2012;Bednaršek & Ohman, 2015;Maas et al., 2016;Moya et al., 2016). However, their usefulness as bio-indicators of the effects of ocean acidification is compromised by limited historical context for understanding species-specific long-term exposure to variations in ocean chemistry. Accurate knowledge of their taxonomy, genetic diversity and biogeography is the essential first step to predicting ecological and evolutionary responses to climate change.
We illustrate an integrative taxonomic approach using the genus Diacavolinia, shelled pteropods with particularly problematic systematics that are usually not identified below genus level, or are listed as Cavolinia sp. or Diacavolinia longirostris (Blainville, 1821) in recent studies (Jennings et al., 2010;Roger et al., 2011;Corse et al., 2013). A new taxonomic assessment of species boundaries in Diacavolinia pteropods is important, because they are morphologically diverse, and some taxa occur in low pH ocean regions, including the California Current coastal upwelling ecosystem (Maas et al., 2013). A study of Diacavolinia pteropods from Australian tropical waters found a significant increase in shell porosity along with a 10% local decline in the aragonite saturation level between the 1960s and 2000s (Roger et al., 2011), suggesting sensitivity of this taxon to contemporary changes in the aragonite saturation state of the ocean.
Previously known as a single species of Cavolinia [Cavolinia longirostris (Blainville, 1821)], Diacavolinia was described as a separate genus by Van der Spoel (1987) on the basis of a distinct shell shape and shell growth compared to Cavolinia taxa. Diacavolinia is the most species-rich genus of pteropods with a total of 24 extant species, of which 18 were introduced by Van der Spoel et al. (1993). Species boundaries were primarily based on shell size and small variations in shell shape that were sometimes found among sympatric taxa. Diacavolinia occurs in tropical and subtropical waters between ~17 and 28 °C in the Atlantic, Pacific and Indian Oceans, as well as in the Red Sea, at depths of ~200 m by day and ~75 m at night (Wormelle, 1962;Van der Spoel, 1967;Bé & Gilmer, 1977). They have complex, bilaterally symmetrical shells (0.4-1.2 cm adult size; Van der Spoel et al., 1993). This shell morphology enables detailed geometric morphometric analyses of shell shape, which can be a powerful tool to distinguish taxa with hard parts (Mitteroecker & Gunz, 2009;Klingenberg, 2010;Burridge et al., 2015). Maas et al. (2013) distinguished four Diacavolinia species from the north-east Atlantic and one from the eastern Tropical Pacific (ETP) based on morphological characteristics. However, they observed that Atlantic specimens were not genetically distinct (<3% divergence), whereas specimens from the Atlantic and Pacific were much more divergent (~19%), based on a fragment of the cytochrome c oxidase subunit I mitochondrial gene (COI). Maas et al. (2013) concluded that broader geographic sampling and a combination of genetic and morphological information are needed to resolve species boundaries in this genus.
We apply an integrative approach to assess species boundaries in Diacavolinia, with inferences based on compatibility between genetic, geometric morphometric and geographic data ( Fig. 1; Padial et al., 2010). To identify species as accurately as possible across our global collection of material, we link datasets comprising fresh specimens for which both genetic and morphometric information is available to morphometric information from museum specimens (969 specimens), including holo-and/or paratype specimens for 14 described species. We aim to (1) develop an objective method for identifying species boundaries by combining incomplete and varied datasets, (2) assess species boundaries and distribution patterns of Diacavolinia taxa by applying an integrative framework of genetic, morphometric and geographic information, and (3) examine consistency of results obtained with this framework across the 24 Diacavolinia taxa as described by Van der Spoel et al. (1993). We find evidence to support a reduction in the number of Diacavolinia species, with at least eight of the species described by Van der Spoel likely representing taxonomic over-splitting. We provide systematic and biogeographic descriptions of the global component of species in this complex genus.

SpecimenS
We included 969 Diacavolinia specimens in this study, with collections from 152 locations between 40°N and 35°S in the Atlantic, Pacific and Indian Oceans (Fig.  2). Of these, 263 fresh specimens suitable for genetic analysis were obtained from 40 Atlantic, 27 Pacific and seven Indian Ocean locations (Table 1). Our fresh material was collected during nine oceanographic expeditions between 2001 and 2012 (Supporting Information, Supplementary Information S1). Plankton nets used across expeditions varied, but our study did not require quantitative sampling and identical net mesh sizes. For example, paired bongo (200-μm, 333μm mesh) and Rectangular Midwater Trawl (RMT1, 333-μm mesh) nets were used in the epipelagic and upper mesopelagic zone during night time on the AMT22 (Atlantic Meridional Transect) expedition in 2012. Specimens from the VANC10MV expedition in 2001 were collected using a bongo net (333-μm mesh) and a 1-m ringnet (333 μm) was used during the COOK11MV and COOK14MV expeditions in 2001. This information is not available for all collected material. All fresh specimens were preserved in 96% ethanol and stored at -20 °C. We also examined 706 specimens from museum collections at the Naturalis Biodiversity Center, Leiden, The Netherlands (NBC) and Zoological Museum of the University of Copenhagen, Denmark (ZMUC). These museum specimens were collected at 78 locations during ten expeditions between 1909 and 1993, and stored in 70% ethanol following initial fixation in formaldehyde (Table 1). Most of the museum specimens (N = 425) were collected during the Danish DANA expeditions between 1921 and 1932 (Supporting Information, Supplementary Information S1). The available museum specimens were identified by Van der Spoel et al. (1993) as 23 of the 24 described Diacavolinia taxa, including holo-and/or paratype specimens (N = 79) of 14 taxa. By examining and including all museum specimens used by Van der Spoel (1993) that could be retrieved, we provide a critical link between prior and current work ( Table 2). Use of this historical material enabled us to make direct comparisons of species boundaries as identified by our methods versus those considered by Van der Spoel et al. (1993). All taxa were represented by the NBC and ZMUC collections, except Diacavolinia robusta Van der Spoel et al., 1993. integrative approach to aSSeSSing SpecieS boundarieS We combined genetic, geometric morphometric and geographic observations on single Diacavolinia specimens wherever possible, following the approach outlined in Figure 1. The information obtained for each specimen varied, but this framework allowed for the combination of partial observations from each specimen. Morphometric measurements consisted of a partial shell shape outline of 49 landmarks (LMs) in lateral orientation and/or 23 LMs in a ventral orientation per specimen (for photographs and geometric morphometric data see: Supporting Information, Supplementary Information S2-3). For some specimens it was possible to obtain an additional 15 ventral LMs to outline the shape of the ventral lip, in cases where the soft tissue did not obscure the ventral lip (also see next paragraph; Fig. 3E; Table 1). Phylogenies were inferred from cytochrome c oxidase subunit I mtDNA (COI; 658 bp) and from nuclear large ribosomal subunit 28S (901 bp) gene fragments (Table 1).
Our approach included six steps for the discovery of, and assignment to, species (Fig. 1). The first step consisted of identifying integrative groups by linking genetic and morphometric information using fresh specimens (N = 173). To test for significant morphometric differences between genetic clades with >five specimens, we applied non-parametric permutational multivariate analyses of variance to shell shape or size parameters (PerMANOVA based on Euclidean distances; Anderson, 2001) in PAST v.2.17c (Hammer et al., 2001. The PerMANOVA F-statistic was tested against 10 4 non-parametric permutations. A significantly different result provided evidence for the presence of distinct species, with concordance observed among genetic and morphometric characters. In the second step, we examined the morphospace position of the holo-and paratype specimens identified by Van der Spoel et al. (1993) for which no genetic information is available due to specimen fixation in formaldehyde. Geometric morphometric measurements were obtained for 79 type specimens from 14 described taxa (Tables 1, 2). We applied Linear Discriminant Analysis (LDA) in R v.3.0.1 (R Development Core Team, 2013) to identify morphological species based on types, merge different types into a single morphological species or merge types with integrative species identified in Step 1. We performed separate LDAs for shape and size data of the different orientations to include as many specimens as possible and to limit the presence of correlated lateral and ventral size variables in the same LDA. Morphometric assignment criteria for synonymization with integrative groups and/ or conservation of holo-and paratypes as distinct morphological species were: at least 80% confidence of belonging to a group for lateral 49 LMs and/or 95% for ventral 23 LMs and/or 85% for ventral 38 LMs and no contradictory assignment between orientations of a specimen. These cut-off values were chosen to reflect the relatively higher information content of the shell outline of the lateral compared to the ventral orientation. If types were synonymized, i.e. merged with groups identified in Step 1 or with another group identified in Step 2, we did so for all type specimens of the same species, if the majority of these types were identically assigned, and we included any remaining unassigned types, because they were always from the same location. In this way, we reduced the number of distinct groups identified in Step 2. The third step was to identify morphological species without holo-or paratypes, based on distinct positions in morphospace not covered by specimens from Steps 1 and 2, using LDA results as support for this distinction. The fourth step was the LDA assignment of the remaining specimens for which morphometric information was available to the groups identified in Steps 1-3. Remaining specimens were either non-sequenced fresh or nontype specimens from museum collections. Individuals remained unidentified (Table 1) if they did not meet the assignment criteria or were assigned ambiguously between orientations, which may indicate possible additional species. In the fifth step, we identified possible species based on individuals with genetic, but without morphometric, information. These are treated as separate groups, because their genetic data could not be linked to morphological data from other groups, although they may be synonymous to groups identified in Steps 2-3. Finally, in the sixth step, we Table 1. Overview of Diacavolinia specimens used in this study. For ventral and lateral geometric morphometrics, numbers of specimens for which morphometric data was obtained are indicated per number of landmarks (LM). See Figure 1 for explanation of the identification steps (Steps 1-5). Some museum specimens could not be identified due to lack of morphological data, because shells were too damaged for geometric morphometric analyses. Some fresh specimens could not be identified due to lack of genetic and/or morphological data These bold values are totals of non-bold values below. For example, the row for "Diacavolinia locations" should be in bold and is the sum of "Diacavolinia fresh locations" and "Diacavolinia museum locations", which are listed below "Diacavolinia locations" *Includes 3 Pacific sequences / 1 location from Maas et al. (2013 ). † Outgroup sequences from Burridge et al. (2017 c). plotted sampling locations of all identified species and possible species of Diacavolinia in order to assemble biogeographic distributions for each species based on all extant information.

geometric morphometricS
For quantitative analyses of Diacavolinia shell shapes and sizes, fresh and museum specimens were photographed in lateral (N = 904) and ventral (N = 919) orientations using a Nikon D100 6 mpx camera (Micro-Nikkor lens 55 mm/3.5, aperture f/11, shutter speed 1/1.3 s, ISO 200, fixed zoom) attached to a stand. To standardize the ventral orientation, specimens were mounted on photographic film using 60% methyl glucose. For lateral standardization, we used fine black sand, free from organic material. For geometric morphometric analyses, we selected photographs of all fully developed adults and excluded specimens that were not well-focused or in standard orientation. We also excluded specimens that were damaged or obscured at relevant positions by soft tissue that could not be removed without damaging the shell. Selected photographs were compiled for further analysis using tpsUtil software (Rohlf, 2006). We used a combination of landmarks and semi-landmarks for partially outlining shell shapes in tpsDig (Rohlf, 2006) to cover as much shell-shape variation as possible for as many specimens as possible (Gunz & Mitteroecker, 2013). To assess the shape variation in the laterally photographed shells, two partial outlines were created, which were connected at the caudal joint between the ventral and dorsal parts at the bottom of the shell (for shell anatomy see: Fig. 3). The first partial outline started at the position of maximum width of the ventral part, ended at the caudal joint and was standardized per specimen to 15 LMs separated by equal length. The second partial outline started at the top of the shell rostrum, ended at the caudal joint via the dorsal part of the shell and was standardized to 35 LMs. Because of the mutual landmark at the caudal joint, this resulted in a total lateral outline of 49 LMs. Creating the first partial outline of 15 LMs was possible for 752 specimens and creating both (49 LMs) was possible for a subset of 549 specimens (Table 1). To assess shape variation in a ventral orientation, two LMs were placed at the left and right gutter corners and a third LM was placed at the caudal joint. Subsequently, left and right partial outlines of ten LMs each were generated from the lock areas left and right of the shell aperture downwards to the closest position without overlap between ventral and dorsal shell parts, the upper start of the dorsal spine surface (see Fig. 3). This resulted in a total of 23 LMs for 646 specimens. No landmarks were created on the two spine tips because these were often damaged. For a subset of 314 specimens it was possible to create an additional 15 LM outline of the ventral lip, resulting in a total of 38 LMs. We used tpsRelw software (Rohlf, 2006) to rotate, translate and scale LM coordinates through generalized least-square Procrustes superimposition (GLS; Kendall, 1977). This provided centroid sizes, a size measure depending on the surface area within all LMs and multiple relative warp (RW) axes containing information on shape variation, with the first RW containing most information. The morphospace of the lateral orientation was represented by 26 relative warps (RWs) for 15 LMs and 94 RWs for 49 LMs. In ventral orientation, there were 42 RWs for 23 LMs and 72 RWs for 38 LMs.
To test for repeatability of RWs, a selection of 19 museum specimens was photographed in two subsequent series for lateral 15 LMs and 49 LMs and ventral 23 LMs, of which ten could also be used for ventral 38 LMs. Intra-class correlation coefficients (ICCs) between the two series were calculated for the centroid sizes and first ten RWs in PAST v.3.0 (Hammer et al., 2001). RWs were considered repeatable when ICC > 0.80, and only repeatable RWs were used in further analyses of shell shape. Centroid sizes were always repeatable (ICC > 0.99). For 15 LMs, the first two RWs were repeatable (ICC > 0.94) and contained 95.96% of the shape variation for this part of the shell. For 49 LMs, RWs 1-8 and 10 were repeatable (ICC > 0.91), accounting for 98.29% of the shell-shape variation. In a ventral orientation, repeatable RWs for 23 LMs were 1-5 and 8 (ICC > 0.83; 92.36% of shape variation) and 1-5 (ICC > 0.81; 83.17%) for 38 LMs, respectively.

geneticS
To assess genetic diversity and phylogenetic relationships in Diacavolinia, we obtained 86 COI mtDNA (GenBank accession numbers MF974762-MF974847) and 138 28S DNA (GenBank accession numbers MF974624-MF974761) sequences from a total of 173 specimens, following photography of shells of the same individuals for morphometric measurements. Tissue fragments of one mm 3 for DNA extraction could only be obtained by damaging the shells. DNA extraction was performed using the EZNA Mollusc DNA Kit (Omega Biotek, 2013), as recommended by Maas et al. (2013). We followed the manufacturer's recommended methods without freeze-drying of tissue before DNA extraction.
Forward and reverse COI and 28S sequences were assembled in MEGA v.6.0 (Tamura et al., 2013) and CodonCode Aligner v.4.1 (CodonCode Corporation, USA, 2013). For 28S, double peaks were registered as ambiguous when apparent in both forward and reverse sequences and with a secondary peak that was at least one-third of the height of the primary peak. Assembled sequences were aligned using the MAFFT algorithm (MAFFT v.7, 2013) and their identities as shelled pteropods were checked by BLAST against the NCBI nt database (Altschul et al., 1997). Three Pacific Diacavolinia sequences from Maas et al. (2013; GenBank accession numbers JX183614-JX183616) were included in the COI alignment, as well as two Atlantic and two Pacific Cavolinia uncinata (d'Orbigny, 1836) specimens (Burridge et al., 2017c; GenBank accession numbers MF048915-MF048918 for COI and MF048968-MF048971 for 28S) in both alignments.
Maximum likelihood (ML; Felsenstein, 1981) was used to infer phylogenetic relationships among specimens for both the COI and 28S alignments and to distinguish clades with high bootstrap s u p p o r t . Fo r C O I , w e u s e d a G e n e r a l T i m e Reversible (GTR) substitution model with different evolutionary rates for the three codon positions (CP), because this is a biologically realistic model for protein coding sequences (Shapiro et al., 2006). GTR with a proportion of invariable sites (I) and gamma distributed rate variation among sites (Γ) was selected from 24 models using the Akaike Information Criterion (AIC) in JModelTest v.2.1.7, in which CP-models were not available (Darriba et al., 2012). For 28S, the most appropriate substitution model selected using AIC was GTR + I. Molecular phylogenies were inferred using maximum likelihood followed by non-parametric bootstrap analysis with 1000 replicates in RaxMLGUI v.1.3 (Stamatakis, 2006;Silvestro & Michalak, 2012).
We identified molecular operational taxonomic units (MOTUs) based on COI and 28S sequences separately, using the online version of ABGD (Automatic Barcode Gap Discovery; Puillandre et al., 2012) and the K2P evolutionary model (Kimura, 1980), the most complex model available in the online version. We also applied the Generalized Mixed Yule Coalescent (GMYC) approach to delimiting species (single threshold; Drummond et al., 2012;Fujisawa & Barraclough, 2013), using an ultrametric tree inferred in BEAST v.2.3.3. Trees were inferred using the HKY substitution model, relaxed log normal clock, Yule model and with an MCMC run for 50 million generations (ten initialization attempts, sampling every 2000 steps for COI and every 5000 steps for 28S, 25% burn-in). To quantify differences between and within MOTUs and genetic clades, we calculated pairwise genetic distances between COI haplotypes as well as 28S alleles using the K2P + Γ distance model of evolution that assumes equal evolutionary rates among all transitions as well as among all transversions (Kimura, 1980) in MEGA v.6.0 (Tamura et al., 2013). Among the available models in MEGA v.6.0, this model most closely represented our sequence data according to the AIC for COI and 28S. We reconstructed alleles from 28S genotypes using the PHASE algorithm (Stephens et al., 2001;Stephens & Donnelly, 2003) in DnaSP v.5 (Librado et al., 2009).
Integrative species were separated based on MOTUs with a pairwise genetic distance of >2% divergence for COI, combined with the ability to distinguish groups morphologically using geometric morphometric methods. If COI sequences were not available, additional groups were identified based on the presence of unique mutations in 28S, a much more slowly evolving gene (e.g. Burridge et al., 2015). If morphologies between different MOTUs, as identified through ABGD and GMYC, could not be distinguished, they were treated as a single group with the notice of possible cryptic species within this group.

RESULTS
integrative approach to aSSeSSing SpecieS boundarieS Geometric morphometric data demonstrate large variation in shell sizes and shapes among Diacavolinia specimens, especially in the Indo-Pacific (Figs 4, 5; Supporting Information, Supplementary Information S4). The first three lateral RWs captured 43.38, 21.95 and 14.97% of the total shell-shape variation, respectively (Fig.  Step 1: identification of integrative species We identify six integrative species based on combined genetic and geometric morphometric information of fresh specimens in Step 1: Groups 1 (Atlantic), 3-5 (Indo-Pacific), 6 (Pacific) and 7 (Indian Ocean; Fig. 6; groups numbered sequentially for Atlantic followed by Indo-Pacific throughout). The COI phylogeny includes Groups 1 (N = 56; 1 MOTU based on ABGD), 3 (N = 21; 1 MOTU), 5 (N = 6; 2 MOTUs), 6 (N = 1; 1 MOTU) and 7 (N = 1; 1 MOTU; Fig. 6A). Within Group 5, a single specimen is assigned to a different MOTU based on ABGD of COI sequences (5B in Fig. 6), but because no geometric morphometric differences could be detected, it is not treated as a separate group. Groups 3 and 5 are also separated by GMYC (support of 1 and 0.92, respectively), as well as Group 12 (support of 1; discussed further in Step 5). Group 1 appears to have some internal structure (subclade support of 0.55), but all sequences in Group 1 are separated from any of the Indo-Pacific groups according to GMYC. The two ABGD MOTUs in Group 5 are not supported by GMYC and neither are the single specimens representing Groups 6, 7 and 13. However, this is most likely explained by the fact that the GMYC method relies on distinguishing the transition from inter-to intraspecific branching within an ultrametric tree and in the absence of multiple individuals, the method cannot detect this transition. Hence, based on results from ABGD, pairwise genetic distances and geometric morphometric information we consider these groups to be separate species. For Group 4, COI sequences failed to amplify in PCR. Groups 12 and 13, from which we (Group 13) and Maas et al. (2013;Group 12) obtained COI sequences, but which lack morphological data due to unavailability of shells (Group 12) or damaged shell characters (Group 13), are discussed in Step 5. Pairwise genetic distances of COI haplotypes between these groups are high with averages of 18.6-43.8% (Table 3). Pairwise distances are small in Groups 1 (average 1.4%, range 0.0-2.7%) and 3 (0.4%, 0.0-1.1%) and larger in Group 5 (6.3%, 0.5-15.6%). Without the single specimen assigned to another MOTU, pairwise distances in Group 5 average 2.0% (0.5-3.0%; 5A in Fig. 6). The Atlantic COI sequences from Maas et al. (2013) all correspond to Group 1 (NCBI BLAST search). Other groups with COI sequences, but without morphological data, are discussed in Step 5. The 28S phylogeny includes Groups 1 (N = 35), 3 (N = 17), 4 (N = 78), 5 (N = 6), 6 (N = 1) and 7 (N = 1) and no additional groups without morphological data are identified based on this genetic marker. Two well-supported Diacavolinia clades are present, each representing one 28S MOTU as identified by ABGD (Fig. 6B). GMYC distinguished 16 MOTUs and, although their composition was not incongruent with the ABGD MOTUs, they are probably an artefact of the low levels of polymorphism in this gene. The first clade contains Groups 1, 3, 6 and 7 (N = 54) and spans all three oceans. The second clade consists of Groups 4 and 5 (N = 84) and is restricted to the Indo-Pacific. The average genetic distance between alleles of the two clades was 2.4% (1.7-3.0%). Within the first clade, average genetic distance is 0.2% (0.0-0.6%) and within the second clade it is 0.1% (0.0-0.4%). Based on 28S, the groups in the first clade cannot be distinguished. Groups 4 and 5 of the second clade can always be distinguished from each other (0.3%, 0.2-0.4% genetic distance between alleles; Table 3).
Across all shell orientations, LDA demonstrates a 100% correspondence between geometric morphometric and genetic data for Groups 1 and 5. Accuracy of assignment is 92.9% for Group 3 and 95.9% for Group 4, with the remaining individuals not reaching all assignment criteria for unambiguous identification (Table 4). We obtained geometric morphometric measurements of    Table 4. Overview of Diacavolinia groups identified in this study following the integrative taxonomic steps outlined in Figure 1. The Linear Discriminant Analysis (LDA) accuracy depicts the correspondence between geometric morphometric and genetic data (Step 1, S1), the correspondence between morphospace position of holo-and paratypes and their LDA identification (Step 2, S2), or the correspondence between estimated morphospace position and their LDA identification (Step 3, S3) following the morphometric assignment criteria (see text).
Step 4 (S4) depicts the numbers of remaining specimens assigned to groups identified in Steps 1-3 using LDA. A = Atlantic; P = Pacific; I = Indian Ocean. See Figure 8 for taxonomic implications Assignment Specimens Stations Steps 1 and 5 Steps 2 and 3 Step 4 Group (Evidence: G = Genetic; M = Morphometric) Total Atlantic Pacific Indian Total Ocean (N specimens) LDA accuracy (N) Species according to Van der Spoel et al. (1993;Ocean, N specimens) LDA accuracy (N) Species according to Van der Spoel et al. (1993; Ocean, N specimens)

Atlantic Group 1 (G,M)
Integrative spe- Steps 2 and 3 Step 4 Group (Evidence: G = Genetic; M = Morphometric) Total Atlantic Pacific Indian Total Ocean (N specimens) LDA accuracy (N) Species according to Van der Spoel et al. (1993;Ocean, N specimens) LDA accuracy (N) Species according to Van der Spoel et al. (1993; Ocean, N specimens)  Additionally, Group 1 is significantly smaller than Group 3 in this orientation (P < 0.01; F = 9.318).

Group 4 (G,M)
Step 2: identification of morphological species with holo-and paratypes From the Atlantic Ocean, holo-and paratypes of D. atlantica Van der Spoel et al., 1993 and D. limbata f. africana Van der Spoel et al.i, 1993 are identified as Group 2 with 100% confidence. Holo-and paratypes of D. constricta Van der Spoel et al., 1993, D. deblainvillei Van der Spoel et al., 1993, D. deshayesi Van der Spoel et al., 1993 and D. ovalis Van der Spoel et al., 1993 are placed in Group 1 with a confidence of 96.4% (Figs 1, 4A, 5A; Table 4). We did not use the ventral orientation with 38 LMs for identifying Atlantic specimens because this part of the shell is always obscured by soft tissue for Group 2.
In the Pacific Ocean, D. trangulata Van der Spoel et al., 1993 holo-and paratypes are identified as Group 8 with 100% confidence. Based on geometric morphometric data from the Indo-Pacific, seven taxa are placed into groups from Step 1 (Figs 1, 4B, 5B, C; Table 4; Supporting Information, Supplementary Information S4). Diacavolinia elegans Van der Spoel et al., 1993 andD. vanutrechti f. vanutrechti Van der Spoel et al., 1993 paratypes from the Pacific fit into Group 3 (75% confidence). The Pacific D. bandaensis Van der Spoel et al., 1993, D. grayi Van der Spoel et al., 1993 andD. pacifica Van der Spoel et al., 1993 holo-and paratypes as well as the Indian Ocean D. bicornis Van der Spoel et al., 1993 andD. striata Van der Spoel et al., 1993 paratypes are identified as belonging to Group 5 (64% confidence). The higher fraction of ambiguous type specimens in Group 5 may be due to the relatively large morphospace covered by this group, in which we cannot distinguish any subgroups, and for which at least two MOTUs are identified by ABGD based on COI sequences.  Step 3: identification of morphological species without types We distinguish three additional morphological species based on the morphospace position of non-type museum specimens from the Indo-Pacific (Groups 9-11; N = 26; Figs 1, 4B, 5B,C; Supporting Information, Supplementary Information S4). For all three groups, there is a 100% correspondence between their estimated position in morphospace and their LDA assignment (Table 4). Non-type specimens were previously identified by Van der Spoel et al. (1993) as Pacific D. strangulata (Deshayes, 1823) (Group 9), Pacific D. mcgowani Van der Spoel et al.i, 1993 and D. longirostris (Group 10) and D. longirostris from the Indian Ocean (Group 11). Group 9 (N = 3 identified in Step 3) is distinguished by a strongly ventrally directed shell rostrum and large size, Group 10 (N = 12) is identified by large spines and a triangular appearance in a ventral orientation and Group 11 (N = 11) is recognized by a subtle outer hump combined with a ventrally directed shell rostrum.
Step 4: assignment of remaining specimens to morphogroups Of our fresh specimens with morphometric, but without genetic, data, all Atlantic specimens (N = 33) are assigned to Group 1 by LDA. Pacific specimens (N = 15; 94.8%) are assigned to Group 4 (1 specimen was ambiguous). We thus infer that a representative number of specimens from Groups 1 and 4 were sequenced in Step 1. Of our remaining museum specimens, 423 specimens (88%) are successfully assigned to groups and 58 (12%) are not (Table 4) but may represent additional diversity.
Step 5: identification of possible species Two additional groups are identified based on COI sequences alone, each representing one MOTU (Groups 12 and 13; Fig. 6A; Table 4). Group 12 contains all three eastern Tropical Pacific (ETP) sequences from Maas et al. (2013;listed as D. vanutrechti). Group 13 contains a single sequence from the Indian Ocean near South Africa and its phylogenetic position is between other Diacavolinia groups and Cavolinia uncinata. No morphological information is available for these groups, because Group 12 is solely represented by sequences from Maas et al. (2013) and Group 13 is represented by a single juvenile specimen with a damaged shell.
Step 6: global biogeography Biogeographic ranges of the identified groups are diverse and varied in size (Fig. 7). In the Atlantic Ocean, two groups occur: Group 1 (N = 276) is present

SyntheSiS
Following our integrative approach, a total of 752 specimens (77.6% of 969 available specimens) is assigned to 13 groups (Figs 1, 8; Table 4). We distinguish two groups in the Atlantic Ocean and 11 groups in the Indo-Pacific. We consider there to be sufficient genetic and morphometric evidence for the validity of the Atlantic Group 1 and Indo-Pacific Groups 3-5 as integrative species, with possible additional diversity in Group 5. The morphological variation of Atlantic Group 1 is large compared to all other groups, including Group 5, but the extensive sampling in the Atlantic Ocean, the homogeneous distribution of specimens within the morphospace and their low levels of genetic variation make it unlikely that there are additional Atlantic species. Although we had only one specimen each for Groups 6 and 7, they likely also represent separate species based on genetic and morphometric data. We identify Atlantic Group 2 and Indo-Pacific Groups 8-11 as morphological species, but they were not sampled in our recent collection and we lack genetic information for these taxa. Hence, we could not link them to possible species of Groups 12 and 13, for which we only had genetic information. However, we consider it unlikely that Group 13 is identical to any of the morphological species analysed here based on its South African locality with none of Groups 8-11 occurring nearby. Of the unassigned individuals (N = 217), 158 specimens provide no evidence to determine their position in morphospace (no suitable photographs available), and neither can their molecular phylogenetic position be determined due to absent data. Therefore, additional diversity may be present. A further 59 specimens for which geometric morphometric measurements are available, of which 50 are from the Indo-Pacific, could not be confidently assigned to a group, and may also represent additional diversity. More molecular data are needed to establish the presence of additional diversity. Figure 8 gives an overview of typical adult shell shapes (in lateral orientation) and sizes for Groups 1-11. Important morphological characteristics for distinguishing between Diacavolinia species are the presence of an outer hump (B in Fig. 3), direction of the shell rostrum (A in Fig. 3), shape of the spines (H in  To our knowledge, this is the first time a global collection of samples, including both recent and museum-type specimens of a marine zooplankton group, are combined into a single dataset for integrative taxonomic purposes. We collected morphological and/or genetic information from 969 specimens belonging to the pteropod genus Diacavolinia for which there was morphological and/ or genetic information available. Our inclusion of all extant museum material also examined by Van der Spoel et al. (1993) allowed for direct comparison with previous taxonomies for this genus. In this way, we more accurately and objectively resolved species boundaries and species distributions than has been possible in prior work.
Our findings suggest the presence of 13 species in our available samples, while 24 taxa were originally described by Van der Spoel et al. (1993). Especially in Figure 8. Overview of revised Diacavolinia taxonomy with example specimens of Groups 1-11 shown in a lateral orientation. Species names as (re)described by Van der Spoel et al. (1993) are listed below each group based on holo-and paratype specimens. Species names are indicated with an asterisk (*) if type specimens were included in this study. Species without type specimens are listed below the group to which the majority of specimens identified as such was assigned in Steps 3 and 4. They are listed only if the species was originally described from the same ocean basin. Photo sizes are standardized to the 0.5-cm scale bar. Groups 12 and 13 are not shown because no morphological data was available.
the Atlantic Ocean, the number of species should be reduced from eight to two, with one species comprising the currently described D. constricta, D. deblainvillei, D. deshayesi, D. longirostis and D. ovalis, and the other comprising D. atlantica, D. limbata f. africana and D. limbata f. limbata (d'Orbigny, 1836) (Fig. 8). New taxonomic descriptions should reflect the larger morphological variation covered by each group and follow rules of the International Commission on Zoological Nomenclature (ICZN). Although the overall morphological diversity in the Atlantic is large, we found little structure in our morphometric data and low levels of genetic diversity, in contrast to the Indo-Pacific. Although we had type specimens for only eight of the 16 original Indo-Pacific species, and no specimens at all for D. robusta, we found evidence for 11 species, comprising 13 of the original taxa. Based on our findings, D. angulata (Souleyet, 1852) (Group 4), D. triangulata (Group 8), D. strangulata (Group 9) and D. mcgowani (Group 10) are confirmed as valid species. Nine original taxa should be merged into two species: one species containing D. elegans, D. vanutrechti f. vanutrechti, D. vanutrechti f. meisenheimeri Van der Spoel et al., 1993 and D. souleyeti Van der Spoel et al., 1993 (Group 3) and one containing D. bandaensis, D. bicornis, D. grayi, D. pacifica and D. striata (Group 5). The latter species may comprise additional taxa based on the presence of two MOTUs and a relatively large shell-shape variation compared to other taxa. Furthermore, three Indo-Pacific groups identified in this study require new taxonomic names (Fig. 8).
Despite the global distribution of our samples, some taxonomic uncertainties remain, especially in the Indo-Pacific, ETP and Red Sea. The sampling coverage in the Atlantic Ocean was higher than in the Indo-Pacific. The validity of three taxa as described by Van der Spoel (1993) remains unclear: D. aspina Van der Spoel et al., 1993, D. flexipes Van der Spoel et al., 1993 and D. robusta, described from the Indian Ocean, Red Sea and Pacific Ocean, respectively (Van der Spoel et al., 1993). Material from D. aspina (N = 1; non-type) could not be unambiguously placed in any group by LDA. Diacavolinia robusta specimens were not available. We also had no material identified as D. flexipes from the Red Sea, its type locality (Table 2). However, Janssen (2007b) suggested that the validity of D. flexipes based on subtle morphological differences is doubtful and rather represents intraspecific variation.
The higher overall species diversity in the Indo-Pacific compared to the Atlantic supports the hypothesis of an Indo-Pacific origin for Diacavolinia as proposed by Van der Spoel et al. (1993). Diacavolinia was already present in the Indo-Pacific in the Pliocene (Piacenzian, 3.6-2.6 Mya), based on fossils from northern Philippine sediments described as Diacavolinia pristina (Janssen, 2007a). It is unknown how long Diacavolinia has been present in the Atlantic Ocean.
barrierS to diSperSal Persistent dispersal barriers may limit range shifts of some taxa in response to changing ocean conditions, while other taxa may be able to adapt and occupy new ecological niches. The most important biogeographic barriers for Diacavolinia, as inferred from distinct species assemblages, were between the Atlantic and Indo-Pacific oceans and the East and Central Pacific. Biogeographic distributions of proposed Diacavolinia species were as follows: Atlantic (two endemic species), warm waters south of South Africa (one endemic species), Western Indian Ocean (four non-endemic species), Red Sea and Gulf of Aden (three non-endemic species), Indo-West Pacific (six species, one endemic), Central Pacific (three species, one endemic), coastal waters of the North-West Pacific (three non-endemic species) and ETP (three endemic species). Species distributions described here are less patchy and disjunct compared to Van der Spoel et al. (1993) and may better reflect ecological and/or habitat preferences of Diacavolinia species. The distribution patterns of the revised Diacavolinia species are congruent with several well-known biogeographic provinces for other holoplankton or benthic species with pelagic larvae, and provide important information on the range of environmental variation experienced by each species (e.g. Kulbicki et al., 2013;Bowen et al., 2016;Iacchei et al., 2016).
The Agulhas Current in the Indian Ocean intermittently forces warm eddies into the Atlantic (Hutchings et al., 2009;Villar et al., 2015). Although distributions of two Diacavolinia taxa extended to waters south of South Africa, we did not find evidence for recent dispersal between the Atlantic and Indian Oceans. Hence, the degree of connection of Diacavolinia in the Atlantic and Indo-Pacific basins appears to have been overestimated by Van der Spoel et al. (1993), who reported that five of the 24 originally described Diacavolinia taxa occurred in both the Atlantic and the Indo-Pacific. Likewise, no evidence of recent dispersal from the Indian Ocean into the Atlantic Ocean was found in pteropods of the genus Cuvierina (Janssen, 2005;Burridge et al., 2015). More rigorous molecular examination of other warm-water pteropods may identify higher numbers of endemic species in the Atlantic and Indo-Pacific Ocean basins than have been described to date. Other pelagic examples for which Atlantic taxa are isolated from the Indo-Pacific include atlantid heteropods (Atlanta selvagensis de Vera &Seapy, 2006 andProtatlanta sculpta Issel, 1911;Wall-Palmer et al., 2016, several copepods (e.g. Hirai et al., 2015;Goetze, 2011) and populations of two-wing flyfish (Lewallen et al., 2016). Conversely, leakage of Indian Ocean water into the Atlantic Ocean via Agulhas eddies appears to have increased over the last decade (Biastoch et al., 2009). Sporadic species dispersal through Agulhas rings has been demonstrated, e.g. for moray eels and glasseye fish (Reece et al., 2010;Gaither et al., 2015). Additionally, Villar et al. (2015) found an overlap in MOTUs of the plankton community within Agulhas rings between the Indian and South Atlantic oceans based on metabarcoding.
We found evidence for endemism of Diacavolinia species in the ETP. Other genetic studies of East Pacific plankton showed that some, but not all taxa demonstrated East Pacific endemism, and it is likely that cryptic diversity is present within what are now considered single species or species complexes. Some endemic cryptic diversity within the ETP was found in the Pleuromamma piseki Farran, 1929-P. gracilis Claus, 1863 copepod species complex (Halbert et al., 2013), and morphological divergence was found in Glaucus marginatus (Reinhardt & Bergh, 1864) sea slugs (Churchill et al., 2014b). Some taxa with pelagic larval stages demonstrate East Pacific isolation, such as the reef-building coral Porites lobata Dana, 1846 (Baums et al., 2012) and the spiny lobster Panulirus penicillatus (Olivier, 1791) (Iacchei et al., 2016), but the echinoderm Echinothrix diadema (Linnaeus, 1758) does not (Lessios et al., 1998). Within another Pacific group of holoplanktonic gastropods (heteropods, Pterotracheoidea), Atlanta californiensis Seapy & Richter, 1993 is restricted to the California Current upwelling ecosystem (Seapy et al., 2003;Wall-Palmer et al., 2018).
Red Sea endemism appears to have multiple causes, with a cold, nutrient-rich barrier separating the Gulf of Aden from the rest of the Arabian Sea, and a narrow strait separating the Red Sea from the Gulf of Aden. Moreover, circulation patterns and environmental gradients may provide additional isolating barriers to dispersal (DiBattista et al., 2016). We found indications of isolation for Red Sea and Gulf of Aden Diacavolinia from the Indian Ocean based on geometric morphometric information, because specimens from the Gulf of Aden (N = 3) and entrance of the Red Sea (N = 8) could not be unambiguously assigned to any of the groups, although our sample sizes are low. Because no genetic information is available, we could not infer whether the Red Sea has exported Diacavolinia biodiversity over time, or if it is an area of ongoing speciation and local endemism due to its peripheral position and unique oceanographic conditions. Fossil records have indicated that the latter is more likely for pelagic taxa, because of the loss of plankton diversity in the central Red Sea due to low sea level stands during Pleistocene glacials, and hypersaline conditions caused by almost complete isolation from the Indian Ocean (Fenton et al., 2000). Conversely, there is also evidence that some taxa increased in abundance due to freshwater dilution in the Gulf of Aqaba, such as observed in Creseis pteropods and siliceous diatoms (Reiss et al., 1980;Fenton et al., 2000;Almogi-Labin et al., 2008). Furthermore, some Red Sea lineages are older than their respective sister lineages in the Indian Ocean, suggesting export of biodiversity from the Red Sea, such as for some reef fishes (DiBattista et al., 2013). Peripheral speciation was observed for the spiny lobster Panulirus penicillatus with a nine-month larval stage, as well as for several reef fish taxa (Liu et al., 2014;Fernandez-Silva et al., 2015;Iacchei et al., 2016).
The high number of Diacavolinia species in the Indo-West Pacific (IWP) reflects the high overall marine diversity in the area (e.g. Renema et al., 2008;Becking et al., 2016). We found that the most common taxa were distributed across a wide range in the coastal waters of the north-west Pacific, central Pacific and Indian Oceans. Present-day waters in the Indo-West Pacific are deep enough to enable diel vertical migration in epi-and upper mesopelagic zones. However, more genetic sampling in this area is needed to further resolve species boundaries, clarify the high genetic and presumably cryptic diversity within Group 5, and explore IWP endemism.
Across all oceans, we observe no obvious equatorial dispersal barriers separating the distributions of extant Diacavolinia species, but our Group 2 (no genetic information available) only occurs in subtropical Atlantic waters north of the equator. We also found no evidence for equatorial genetic breaks in other Diacavolinia groups. Equatorial dispersal barriers may be important drivers of pelagic evolution for other taxa, such as shown for subtropical Atlantic copepods (Goetze et al., , 2017 and Cuvierina (Burridge et al., 2015).

CONCLUSIONS
Combining varied datasets in an integrative taxonomic framework may be suitable for a wide array of morphologically diverse marine taxa and is an important first step to predicting species-specific responses to climate change. Museum collections prove to be an invaluable resource for assessing species boundaries and biogeographic distributions, enabling comparison of modern findings with prior works based on the same specimens. We assessed species boundaries in Diacavolinia pteropods based on rigorous sampling over a period of 104 years and data collection across the Atlantic, Pacific and Indian Oceans, which would have never been possible within a single research project. Our results show that taxonomic revisions of Diacavolinia are needed, and will be reported on in subsequent work. However, not all species boundaries were resolved in this study and some species were only sampled once or not at all in recent years. New specimens could be added to the current framework based on new species descriptions, standardized photographs and information on COI or 28S. It is also remarkable that a morphospecies from the Atlantic Ocean (Group 2, suggested species name: D. atlantica; Fig. 8) was abundant in museum samples but has not been encountered in recent years, despite multiple sampling expeditions in this region. Although it is common to uncover new species in zooplankton when integrative approaches are applied (e.g. Wall-Palmer et al., 2018), in this study, we observed fewer genetic lineages than expected based on prior taxonomic descriptions of morphological species.