Exploring Genetic and Morphological Integrity across Ocean Basins: A Case Study of the Mesopelagic Shrimp Systellaspis debilis (Decapoda: Oplophoridae)

: Plankton communities often consist of cosmopolitan species with an extensive gene ﬂow between populations. Nevertheless, populations of some plankton species are genetically structured, owing to various barriers such as ocean currents, hydrological fronts, and continents. Drivers that could explain the genetic structures of most mesopelagic species remain unknown on an ocean-basin scale, and our study aims to analyze the genetic and morphological differences between populations of a cosmopolitan mesopelagic shrimp, Systellaspis debilis, from the Southern and Northern Atlantic Ocean, and the Southwest Indian Ocean. We analyzed the ITS-1 and COI markers of 75 specimens and assessed the genetic integrity and within-species variability of these genes. We also coded 32 morphological characteristics in 73 specimens, analyzed their variability, and assessed the correlation between morphological and genetic characteristics using a Redundancy analysis and Mantel test. Systellaspis debilis was genetically cohesive across the whole Atlantic and Southwest Indian Oceans, which is possibly a result of an intensive gene ﬂow through ecological barriers, the resistance of species to hydrological gradients, a purifying selection of mitochondrial genes, etc. In contrast, we found signiﬁcant morphological differences between populations from different regions, which mirrors morphological diversiﬁcation and calls for further genomic approaches in order to understand the basis of these variations and uncover potential local adaptations.


Introduction
Oceans cover most of the Earth's surface area and habitat volume; the vast deeppelagic habitat between the sunlit layers (upper 200 m) and the seafloor is the largest and least-understood environment on our planet [1,2].This habitat contains the mesopelagic (from 200 m to ca.1000 m depth) and the deeper bathy and abyssopelagic zones.Our limited knowledge of these ecosystems is increasingly problematic as they may be vulnerable to global issues such as climate warming, deoxygenation, acidification, commercial fishing, seabed mining, and other threats with an unknown potential for feedback to the climate system [3,4].Albeit greatly underexplored, the mesopelagic zone provides a better chance for ecosystem analyses than deeper layers that require even more time-and cost-consuming efforts.Recent analyses, based mainly on an expert opinion on the distributional patterns of pelagic fauna relative to environmental proxies, allowed a global biogeographic classification of the mesopelagic zone [2].The same authors declared that "work remains to be done to produce a comprehensive and robust mesopelagic biogeography," and this work should be based on numerous empirical observations of the factors driving biodiversity of individual species within the mesopelagic zone.Attention should be paid to zooplankton, which are the key element in the mesopelagic zone because they are the basic trophic link primary producers with larger predators, and abundant enough to be representatively sampled (e.g., [5,6]).
Studies on genetic diversity are focus on finer and promising tools for a deeper insight into pelagic biogeography; this tool, however, has so far been applied to a limited number of zooplankton species and showed that the patterns of genetic structuring of populations are species-specific [17,18,23].In other words, we need much more research on individual species for a proper understanding of drivers that explain the true (genetic and morphological) biodiversity of the mesopelagic zone.
In this paper, we make the next, among many, step in this direction, and focus on the population structure of a cosmopolitan species, Systellaspis debilis (A.Milne-Edwards, 1881), that makes a significant contribution to mesopelagic ecosystems [24] and is the fourth most common pelagic shrimp in the Atlantic Ocean [25].In contrast to previous studies on mesoplankton, this is a macroplanktonic decapod, a group still unexplored in this context despite their prominent role in the mesopelagic zone (40% of the total mesopelagic plankton biomass [26]).Systellaspis debilis occurs in many mesopelagic biogeographical provinces (sensu Sutton et al. [2]) in the Atlantic, Indian, and Pacific Oceans between 63 • N to 58 • S [27].Such an extensive range of species always raises questions about genetic homogeneity and population structure.In this paper, we describe and analyze the genetic and morphological diversity of S. debilis in order to assess the degree of isolation between populations from various basins.Due to the high requirement of material to analyze (undamaged specimens for morphological analysis, "fresh" alcohol-fixed individuals for genetic analyses), our studies were restricted to the Northern Atlantic, Southern Atlantic, and Southwestern Indian Ocean.
We tested the hypothesis that populations of S. debilis are genetically and morphologically distinct in these three ocean basins, and analyzed the accordance of their geographic distribution with the proposed scheme of mesopelagic zonation [2].In order to test our assumption, we assessed the distribution of genetic and morphological variability in S. debilis populations across the Atlantic and in the Southwest Indian Oceans.

Material
The material was collected in the Atlantic Ocean and the Southwestern Indian Ocean during cruises from 2013 to 2020 (Figure 1, Table S1) with a Bogorov-Rass plankton net (mouth area 1 m 2 , 500 µm mesh size) and an Isaacs-Kidd midwater trawl (mouth area 5.5 m 2 , mesh size 5 mm).A total of 75 specimens of S. debilis were identified using the key of Lunina et al. [28], fixed in 96% ethanol just after retrieval, and stored at −20 • C in the laboratory for further analysis.The COI gene was successfully sequenced in all of the 75 specimens in our collections (46 specimens from the North Atlantic, 5 from the South Atlantic, and 24 from the Indian Ocean); an additional 31 sequences were mined from GenBank and added to the dataset.and HCOI 2198 (TAAACTTCAGGGTGARDAAAAAATCA) [29], or decapod-specific primers COL6 (5′-ACAAATCATAAAGATATYGG-3′) and COH6 (5′-TADACTTCDGGRTGDRDAAARAAYCA-3′) in cases where the former failed.The primers ITS1FW (5'-CACACCGCCCGTCGCTACTA-3') and ITS3R (5′-TCGACSCACGAGCCRAGTGATC-3′) [30] were used to amplify the ITS1 gene.PCR reactions were made in a reaction volume of 20 µL, containing 2.4 µL of the Encyclo Plus PCR kit (Eurogen, Russia), 0.2 µL of each primer, 1.6 µL of DNA template, 15.3 µL Mil-liQ water, and 0.3 µL of 50 X Encyclo polymerase (Eurogen, Russia).The PCR cycling profiles and annealing temperatures are listed in Table S2.The PCR products were purified and sequenced with the same primer sets on an ABI Prism 3500 xl genetic analyzer in the Resource Center Development of Molecular and Cellular Technologies of Saint Petersburg State University.Forward and reverse COI and ITS1 sequences were assembled in Geneious ® 7.1.3.and manually treated for ambiguities and heterozygotes (in the case of ITS1).Additionally, COI sequences were checked for stop codons using Geneious ® 7.1.3software.All sequences were deposited in the NCBI GenBank database [31] (Table S1; accession numbers: OR398994-OR399068 and OR415900-OR415922).

Sequence Alignment and Phylogenetic Analyses
In addition to our material, we used all available COI sequences of S. debilis and the most closely related species, S. liui (no.KT946751), deposited in GenBank (https://www.ncbi.nlm.nih.gov/genbank/,accessed on 1 September 2023).Two species of the superfamily Oplophoroidea, S.curvispina (no.KP076159) and Acanthephyra quadrispinosa (no.KP076178), were chosen as outgroups to root the tree.Multiple alignments of

DNA Extraction, Amplification, and Sequencing
DNA was extracted either from the fifth pair of the pleopods or from the pleonic muscle tissue using the IG-Spin™ DNA Prep 200 kit for DNA extraction following the manufacturer's protocol.The extracted DNA was used as a matrix for the amplification of the mitochondrial cytochrome c oxidase subunit I gene fragment I (COI), and the nuclear gene of the first internal transcribed spacer (ITS1).PCR amplification of the COI gene fragment was accomplished with the universal primers LCOI 1490 (GGTCAACAAAT-CATAAAGATATTGG) and HCOI 2198 (TAAACTTCAGGGTGARDAAAAAATCA) [29], or decapod-specific primers COL6 (5 -ACAAATCATAAAGATATYGG-3 ) and COH6 (5 -TADACTTCDGGRTGDRDAAARAAYCA-3 ) in cases where the former failed.The primers ITS1FW (5'-CACACCGCCCGTCGCTACTA-3') and ITS3R (5 -TCGACSCACGAGCCRAG TGATC-3 ) [30] were used to amplify the ITS1 gene.PCR reactions were made in a reaction volume of 20 µL, containing 2.4 µL of the Encyclo Plus PCR kit (Eurogen, Russia), 0.2 µL of each primer, 1.6 µL of DNA template, 15.3 µL MilliQ water, and 0.3 µL of 50 X Encyclo polymerase (Eurogen, Russia).The PCR cycling profiles and annealing temperatures are listed in Table S2.The PCR products were purified and sequenced with the same primer sets on an ABI Prism 3500 xl genetic analyzer in the Resource Center Development of Molecular and Cellular Technologies of Saint Petersburg State University.Forward and reverse COI and ITS1 sequences were assembled in Geneious ® 7.1.3.and manually treated for ambiguities and heterozygotes (in the case of ITS1).Additionally, COI sequences were checked for stop codons using Geneious ® 7.1.3software.All sequences were deposited in the NCBI GenBank database [31] (Table S1; accession numbers: OR398994-OR399068 and OR415900-OR415922).

Sequence Alignment and Phylogenetic Analyses
In addition to our material, we used all available COI sequences of S. debilis and the most closely related species, S. liui (no.KT946751), deposited in GenBank (https: //www.ncbi.nlm.nih.gov/genbank/,accessed on 1 September 2023).Two species of the superfamily Oplophoroidea, S.curvispina (no.KP076159) and Acanthephyra quadrispinosa (no.KP076178), were chosen as outgroups to root the tree.Multiple alignments of all sequences were made in Geneious ® 7.1.3using the MUSCLE algorithm [32] (25 repeats).The final alignment for the COI fragment was 539 bp and included 109 sequences, and for the ITS1 fragment, 23 sequences of 328 bp.In the case of the ITS1 gene, the sequences were not found in public sources, so only newly generated sequences were analyzed.
The phylogenetic reconstruction of the COI gene by the Maximum Likelihood (ML) estimator was conducted with the RAxML (ver.7.2.8 [33]), using the GTR + G nucleotide substitution model for each codon position.Statistical support was assessed using the bootstrap method involving 1000 replicates.Bootstrap values greater than 70% were considered statistically significant.Before the Bayesian analysis was conducted on the COI dataset, the most appropriate nucleotide substitution models and partitioning scheme were selected for each codon, using the Akaike Information Criterion (AICc) in the PartitionFinder2 software [34]).As a result, the nucleotide substitution patterns were as follows: GTR + I + G for the first codon, GTR + I for the second, and GTR + G for the third.A Bayesian analysis was performed using MrBayes 3.3 software [35].Two parallel calculations of 10,000,000 generations, with tree selection every 1000 generations, were performed, and the first 25% of trees were excluded from the calculation of the posterior probabilities.Bayesian posterior probabilities greater than 0.95 were considered statistically significant.

Morphological Analysis
In order to assess the within-species morphological variability of S. debilis, we selected and coded the 32 most variable characteristics linked to carapace (5 characteristics), pleon (7 characteristics), antenna (1 characteristics), telson (2 characteristics), and pereopods (17 characteristics) (Table 1, Figure 2).The carapace length was measured from the posterior margin of the eye orbit to the dorsal posterior end of the carapace; the carapace height was measured at the highest point.All measurements are presented in Table S3.We coded morphological characteristics in 73 specimens ranging from 3.5 mm to 15.8 mm in carapace length: 43 females, 26 males, and 4 juveniles.
all sequences were made in Geneious ® 7.1.3using the MUSCLE algorithm [32] (25 repeats).The final alignment for the COI fragment was 539 bp and included 109 sequences, and for the ITS1 fragment, 23 sequences of 328 bp.In the case of the ITS1 gene, the sequences were not found in public sources, so only newly generated sequences were analyzed.
The phylogenetic reconstruction of the COI gene by the Maximum Likelihood (ML) estimator was conducted with the RAxML (ver.7.2.8 [33]), using the GTR + G nucleotide substitution model for each codon position.Statistical support was assessed using the bootstrap method involving 1000 replicates.Bootstrap values greater than 70% were considered statistically significant.Before the Bayesian analysis was conducted on the COI dataset, the most appropriate nucleotide substitution models and partitioning scheme were selected for each codon, using the Akaike Information Criterion (AICc) in the PartitionFinder2 software [34]).As a result, the nucleotide substitution patterns were as follows: GTR + I+G for the first codon, GTR + I for the second, and GTR + G for the third.A Bayesian analysis was performed using MrBayes 3.3 software [35].Two parallel calculations of 10,000,000 generations, with tree selection every 1000 generations, were performed, and the first 25% of trees were excluded from the calculation of the posterior probabilities.Bayesian posterior probabilities greater than 0.95 were considered statistically significant.

Morphological Analysis
In order to assess the within-species morphological variability of S. debilis, we selected and coded the 32 most variable characteristics linked to carapace (5 characteristics), pleon (7 characteristics), antenna (1 characteristics), telson (2 characteristics), and pereopods (17 characteristics) (Table 1, Figure 2).The carapace length was measured from the posterior margin of the eye orbit to the dorsal posterior end of the carapace; the carapace height was measured at the highest point.All measurements are presented in Table S3.We coded morphological characteristics in 73 specimens ranging from 3.5 mm to 15.8 mm in carapace length: 43 females, 26 males, and 4 juveniles.
Statistical analyses of the morphological data and comparisons of the morphological and genetic parameters were run using R 4.0.5 [40].Missing characteristics (1.2% of the database) were replaced with their mean values characteristics [41].The juveniles (carapace lengths < 5 mm) were removed from the morphological analysis as the proportions of this species greatly change during ontogenesis.
In order to remove the influence of individual size, we used carapace length and carapace height as predictors in a Redundancy analysis (RDA: [42]), and the other 26 characteristics as dependent variables in the matrix.The analysis was conducted using the RDA function from the "vegan" package [42], yielding both constrained and unconstrained axes.
The canonical axes were influenced by the size of the individuals, while the unconstrained axes provided insight into the structure of the residual matrix from regression models, allowing us to examine the relationship between morphological characteristics without the impact of size.Therefore, we excluded the canonical axes (RDA1 and RDA2) from further analysis and focused on the two most informative unconstrained axes, PCA1 and PCA2, which facilitated a more accurate analysis of morphological characteristics without interference from the influence of individual size.
In order to assess the correlation between the morphological features and genetic characteristics, we used the Mantel test [41] and created two distance matrices.The first matrix included Euclidean distances between individuals in PCA1 and PCA2 space, whereas the second one included square roots of pairwise genetic distances between the sequences of the COI gene.Genetic distances were calculated using the dist.alignment function from the "seqinr" package [43].The mantel correlation between the two matrices was calculated using the Mantel function from the "vegan" package [42].The statistical significance of the test was assessed using the permutation method (9999 permutations).The results of the statistical analyses were visualized using the package "ggplot2" [44].

Genetic Variability and Spatial Structure
The phylogenetic reconstruction retrieved two supported clades, the most abundant, Clade 1 (1/76, Bayesian posterior probabilities/ML bootstrap), comprised 96% of the COI sequences (Figure 3A).This clade included all specimens from the North and South Atlantic (sixty-eight and six, respectively) and twenty-seven specimens from the Indian Ocean.Clade 2, sister of Clade 1, did not gain bootstrap support (0.98/57) and encompassed four specimens of S. debilis collected off the north coast of Madagascar, and one specimen of S. liui (KT946751) from the western Pacific.
Specimens from Clade 1 showed a moderate haplotype diversity (H d ) of 0.547 ± 0.059 (range: 0.611-1.000)and a low nucleotide diversity (π) of 0.0016 ± 0.000 (range: 0.0020-0.0056)across all three regions (Table 2).In the COI minimum-spanning network, 21 unique haplotypes were observed across 102 specimens, with 68 of these representing a shared central haplotype across all three regions (Figure 3B).Specimens from Clade 2 (including S. liui) had higher values of H d (1.000 ± 0.126) and π (0.0122 ± 0.003), and unique haplotypes separated by 29 substitutions from the Clade 1 haplogroup (Figure 3B).The Tajima's D neutrality test resulted in a rejection of the neutral model for Clade 1 overall (D = −2.338,p < 0.001) and the North Atlantic population (D= −1.913, p < 0.05), which is typical of a recently expanded population.
The ITS1 gene marker was analyzed in the specimens from Clade 1, as the Clade 2 representatives were absent in our collection.We randomly sorted five to ten specimens from each region and successfully sequenced ten specimens from the North Atlantic, five from the South Atlantic, and eight from the Indian Ocean.Genetic diversity was very low among the sequences: nineteen out of twenty-three were identical, and the others differed in one to twelve substitutions.
The phylogenetic reconstruction retrieved two supported clades, the most abundant, Clade 1 (1/76, Bayesian posterior probabilities/ ML bootstrap), comprised 96% of the COI sequences (Figure 3A).This clade included all specimens from the North and South Atlantic (sixty-eight and six, respectively) and twenty-seven specimens from the Indian Ocean.Clade 2, sister of Clade 1, did not gain bootstrap support (0.98/57) and encompassed four specimens of S. debilis collected off the north coast of Madagascar, and one specimen of S. liui (KT946751) from the western Pacific.

Morphological Variability
Four characteristics (numerous lateral spines arranged in two or more rows on the telson; the presence of the medial dorsal groove on the scaphocerite; number of the movable spines on the ischium of the third pereopod (anterior row of spines); and number of the movable spines on the carpus of the third pereopod (anterior row of spines) with no variance) were removed from the analysis.The remaining 28 characteristics examined in 73 specimens were used in further statistical analyses.
The RDA model, with carapace length and height as the predictors, was statistically significant (F = 7.3655, p perm = 0.0001, N perm = 9999).The two canonical axes described 18% of the total variance, which suggested that 82% of the morphological variance was not related to body size, and 38% of the residual variability was determined by PCA1 and PCA2 (Table 3).The RDA model showed a positive correlation between body size and some morphological characteristics (Figure 4A), such as number of movable spines in the posterior row of the ischium of the fifth pereopods, and number of serrations on the lateral margin of the fourth and fifth abdominal somites.We excluded the influence of size and further analyzed the residual RDA variability, i.e., the variability of the non-canonical axes, such as Principal Components (PC) (Figure 4B).The RDA model showed a positive correlation between body size and some morphological characteristics (Figure 4A), such as number of movable spines in the posterior row of the ischium of the fifth pereopods, and number of serrations on the lateral margin of the fourth and fifth abdominal somites.We excluded the influence of size and further analyzed the residual RDA variability, i.e., the variability of the non-canonical axes, such as Principal Components (PC) (Figure 4B).  1.
Morphological variability (PC1) was significantly linked to the sampling region (ANOVA: F = 5.306, p = 0.0073) (Figure 4B).Specifically, specimens from the North Atlantic region exhibited lower PC1 values compared to those from the South Atlantic and the Indian Ocean.However, we did not observe any significant relationship between PC2 and the sampling location (ANOVA: F = 0.01, p = 0.99).
Some morphological characteristics showed significant dependence on location (Figure 5).Numbers of spines in the posterior row on the merus of the third and in the anterior row of the merus of the fourth pereopod were significantly higher in the Indian Ocean than in the Atlantic.Similarly, individuals from the South Atlantic had a higher average number of spines in the posterior row on the merus of the fourth pereopod than individuals from other geographic areas.Furthermore, the number of lateral serrations on the pleon on the left side of the fourth abdominal segment of the South Atlantic shrimp was only slightly higher than that of individuals from the Indian Ocean.The lowest average number of spines on the third and fourth pereopod and the number of teeth on the left side of the fourth pleonic somite were observed in the North Atlantic.1.
Morphological variability (PC1) was significantly linked to the sampling region (ANOVA: F = 5.306, p = 0.0073) (Figure 4B).Specifically, specimens from the North Atlantic region exhibited lower PC1 values compared to those from the South Atlantic and the Indian Ocean.However, we did not observe any significant relationship between PC2 and the sampling location (ANOVA: F = 0.01, p = 0.99).
Some morphological characteristics showed significant dependence on location (Figure 5).Numbers of spines in the posterior row on the merus of the third and in the anterior row of the merus of the fourth pereopod were significantly higher in the Indian Ocean than in the Atlantic.Similarly, individuals from the South Atlantic had a higher average number of spines in the posterior row on the merus of the fourth pereopod than individuals from other geographic areas.Furthermore, the number of lateral serrations on the pleon on the left side of the fourth abdominal segment of the South Atlantic shrimp was only slightly higher than that of individuals from the Indian Ocean.The lowest average number of spines on the third and fourth pereopod and the number of teeth on the left side of the fourth pleonic somite were observed in the North Atlantic.

Population Structure of Systellaspis debilis (Clade 1)
In the Atlantic and Southwest Indian Ocean, S. debilis encompasses two divergent and reciprocally monophyletic mitochondrial clades with different geographic distributions.The Atlantic harbors only representatives of Clade 1; since the type locality of S. debilis is the Bahama Channel (North Atlantic), we consider this clade as S. debilis.This clade is genetically homogenous throughout the Atlantic and the Southwest Indian Oceans.The genetic similarity of the COI and ITS1 genes could be the result of several scenarios.
Firstly, there is an intensive gene flow through ecological barriers which are usually impede gene flow between populations of mesoplankton (animals much smaller than the species considered here) [12,13,15,16,[20][21][22].These mesoplankton species, even having haplotypes with panoceanic distribution, show significant variations of haplotype frequencies between oceans or ecoregions [8,12,20,45,46].In contrast to most mesoplankton, S. debilis is a macroplankton species undertaking intensive diurnal vertical migrations through vertical abiotic gradients [47], which makes the species resistant to horizontal gradients of the frontal zones.In addition, a long life cycle (5 to 8 years for oplophoroid shrimps) [48] provides a better opportunity for individual transfer through geographic regions with oceanic currents, and thus also contributes to high levels of gene flow between distant regions.
Another possible scenario suggests that barriers to the gene flow currently exist but were established not long ago, and not enough time has passed to provide genetic differentiation between populations.In this case, purifying selection is effective in eliminating even slightly disadvantageous mutations and maintaining genetic homogeneity in each of the distant populations [49,50].In fact, the purifying selection was thought to be a constraint on genetic diversity and differentiation between two distant populations of this species in the North-West Atlantic [11].Subtle or no genetic differentiation at the global (circumtropical) scale was also reported for some large pelagic fishes and was explained by the large effective size of their populations and/or high capacity for dispersion, which can obscure signals of spatial genetic differentiation [51-53].In order to analyze the relationship between the genetic and morphological characteristics of individuals, we ran the Mantel test that assessed the similarity of the two distance matrixes (genetic distance matrix and spatial distance matrix of PC1 and PC2).The results of this test showed that there is a statistically significant similarity between the two matrixes (r = 0.1791, p = 0.003, 9999 permutations).

Population Structure of Systellaspis debilis (Clade 1)
In the Atlantic and Southwest Indian Ocean, S. debilis encompasses two divergent and reciprocally monophyletic mitochondrial clades with different geographic distributions.The Atlantic harbors only representatives of Clade 1; since the type locality of S. debilis is the Bahama Channel (North Atlantic), we consider this clade as S. debilis.This clade is genetically homogenous throughout the Atlantic and the Southwest Indian Oceans.The genetic similarity of the COI and ITS1 genes could be the result of several scenarios.
Firstly, there is an intensive gene flow through ecological barriers which are usually impede gene flow between populations of mesoplankton (animals much smaller than the species considered here) [12,13,15,16,[20][21][22].These mesoplankton species, even having haplotypes with panoceanic distribution, show significant variations of haplotype frequencies between oceans or ecoregions [8,12,20,45,46].In contrast to most mesoplankton, S. debilis is a macroplankton species undertaking intensive diurnal vertical migrations through vertical abiotic gradients [47], which makes the species resistant to horizontal gradients of the frontal zones.In addition, a long life cycle (5 to 8 years for oplophoroid shrimps) [48] provides a better opportunity for individual transfer through geographic regions with oceanic currents, and thus also contributes to high levels of gene flow between distant regions.
Another possible scenario suggests that barriers to the gene flow currently exist but were established not long ago, and not enough time has passed to provide genetic differ-entiation between populations.In this case, purifying selection is effective in eliminating even slightly disadvantageous mutations and maintaining genetic homogeneity in each of the distant populations [49,50].In fact, the purifying selection was thought to be a constraint on genetic diversity and differentiation between two distant populations of this species in the North-West Atlantic [11].Subtle or no genetic differentiation at the global (circumtropical) scale was also reported for some large pelagic fishes and was explained by the large effective size of their populations and/or high capacity for dispersion, which can obscure signals of spatial genetic differentiation [51][52][53].
The second scenario is supported by the star-like structure of the haplotype network, the lack of transversion mutations, and the negative and significant Tajima's D values.In fact, the low haplotype diversity of S. debilis is unusual for a globally distributed zooplankton species; a similar effect was found only in a few species over a much more limited distribution: the northern krill Meganyctiphanes norvegica [54] and the neritic chaetognath Sagitta setosa [55].Despite their large population size, marine pelagic species may be susceptible to population crashes with measurable effects on their genetic makeup.Bottlenecks resulting from range contractions during the Pleistocene were proposed to have occurred in two copepod and one chaetognath species in the North Atlantic, which displayed lower levels of genetic variation than expected from their estimated population sizes [56,57].Systellaspis debilis may have emerged relatively recently and spread over a huge area due to the possession of some evolutionary/ecological advantage.As there are no fossil records for S. debilis or related species, no correlation can now be found with any specific event in the past.
In the Indian Ocean, most S. debilis specimens were collected between 20 • S and 34 • S within the Southern Indian Ocean and Agulhas Current mesopelagic ecoregions, according to Sutton et al. [2].Only one specimen from GenBank was collected at ~13 • S, northwest of Madagascar, in the Mid-Indian Ocean ecoregion.Surprisingly, the same site harbored four genetically different specimens, which, along with S. liui from the West Pacific, comprised Clade 2. As no specimens of Clade 2 were found in the Atlantic Ocean or south of 20 • S in the Indian Ocean, we suggest that the geographic boundary between both clades occurs between 13 • S and 20 • S in the Indian Ocean.This is consistent with the boundary between Mid-Indian and Southern Indian zones [2].
The genetic break between the clades can be caused by the Agulhas current, which originates at ~27 • S [2] and transports water from the southwest tropical Indian Ocean to the Southern Atlantic, thereby stimulating the dispersal of Clade 1 representatives as it was shown earlier for some marine species [58].Conversely, Clade 2 appears to be confined to the northern Indian Ocean, where the Gyre system prevails [2].As the tropical Indian Ocean likely harbors representatives of Clade 2, the Mozambique Channel acts as a meeting place for both clades.

Morphological Variability of Systellaspis debilis (Clade 1)
Our analyses showed that specimens from the North Atlantic had lower average PC1 values than those from the South Atlantic and Indian Oceans (Figure 4B) and suggested significant morphological differences between shrimp populations in these regions.In particular, the number of spines in the posterior row on the merus of the third pereopod and anterior row of the merus of the fourth pereopod were significantly higher in individuals collected from the Indian Ocean.Individuals from the South Atlantic, on the other hand, had a higher average number of spines in the posterior row on the merus of the fourth pereopod than individuals from other geographic groups.Finally, the lowest average of both the number of spines on the third and fourth pereopod, and the number of teeth on the left side of the fourth segment, were observed in the North Atlantic group, highlighting the differences in morphology across geographic regions.
Alain Crosnier [59] was the first to find variations in the rostrum-carapace length of S. debilis specimens from different geographic locations, including Northern and Southern Madagascar, the Northern and Southern Atlantic, and the Philippines.Our results indicated that the length of the rostrum to the length of the carapace (Lr/Lc) ratio of specimens from the Northern Atlantic was less homogeneous compared to those from the Indian Ocean, but significant variations between the basins that could potentially be correlated with haplotypes were not found.
The observed morphological differences between populations might be influenced by a combination of genetic and environmental factors.The Mantel test suggests that there is a significant similarity between genetic and spatial distance matrices, indicating that the observed morphological variability is likely to be driven, at least in part, by genetic differences between populations.The results suggest a link between the genetic makeup of individuals and their morphology; the exact mechanisms driving this correlation remain unclear and warrant further investigation.
The correlation between genetic and morphological traits does not necessarily imply a direct causal relationship.Other factors, such as environmental conditions and developmental plasticity, could also play a role in shaping the observed morphological variability [60].Future studies that incorporate environmental and developmental factors will be necessary to fully understand the complexity of the relationship between genetic and morphological traits in this species of shrimp.

The status of Systellaspis liui and Related Specimens (Clade 2)
Clade 2 encompasses four specimens deposited in GenBank as 'Systellaspis debilis' [61] and one specimen of S. liui [62].Sha and Wang [62] described Systellaspis liui based on a single female specimen from the Western Pacific (Philippine Sea) and suggested that four specimens of 'S.debilis' from Aznar-Cormano et al. [61] and S. liui are synonyms.According to the original description, genetic and morphological differences between S. debilis and S. liui were sufficient to suggest the new species, but the validity of S. liui was considered to be controversial [28].In fact, morphological variations in S. debilis from various locations [59,62] makes defining S. liui as a distinct species difficult.
Sha and Wang [62] proposed five morphological characteristics to distinguish S. liui from S. debilis.Our morphological analysis showed that four of them were present in 26-100% of observed S. debilis: a medial dorsal groove on the scaphocerite; a carina on the dorsal margin of the third abdominal somite; movable spines on the pereopods; and three teeth on the posterior margin of the fifth abdominal somite.The only morphological characteristic not found in our specimens of S. debilis was the presence of additional spines on the telson, a characteristic not common for the family Oplophoridae and likely attributed to an abnormal specimen [28].
In this study, we compared interspecific and intraspecific K2P distances for all species of the genus Systellaspis, including S. debilis Clade 1, S. debilis Clade 2, and S. liui.The distance between Clade 2 (with and without S. liui) and Clade 1 was 7.1%, which was the lowest observed pairwise distance (ranging from 8.0% to 32.5%) (Table 4).A comparable divergence (8.0%) was observed between Systellaspis braueri and Systellaspis paucispinosa, which also differ mainly in the spination of the telson.
Overall, the observed divergences in COI sequences within and between S. debilis Clade 1 and S. debilis Clade 2 suggest that S. liui is similar to the "S.debilis" of Aznar-Cormano et al. [61], and represents a separate mitochondrial clade.However, further investigations with additional material for morphological studies and the analysis of additional nuclear genes are necessary to clarify the taxonomic status of S. liui and the previously mentioned 'S.debilis' of Aznar-Cormano et al. [61].

Conclusions
Our data indicates that S. debilis is a genetically cohesive species throughout its distribution range in the whole Atlantic and the Southwest Indian Ocean.Populations of S. debilis are genetically homogenous in three geographically distant ocean basins separated by oceanographic fronts, which is unusual for plankton species studied thus far.In contrast to genetic homogeneity, statistically significant morphological differences were found, and populations from the North Atlantic, South Atlantic, and Southwest Indian Oceans differ in the spination of pereopods and the serration of pleonic somites.Scenarios to explain the observed phenomenon include intensive gene flow through ecological barriers owing to resistance to horizontal oceanographic gradients, and long life cycle and/or purifying selection of mitochondrial genes.In both cases, morphological variation between regions may be a result of phenotypic plasticity or have a genetic (not mitochondrial-linked) basis.The use of genomic approaches will clarify this question and unveil finer detail about population structure and possible local adaptations of the species, as was shown for some other pelagic organisms [53,71].We encourage marine biologists to further study the population structure of mesopelagic shrimps that are the key component of deep-sea communities and a target for possible commercial exploitation.
Systellaspis debilis is distinct from both S. liui and the five specimens mentioned as 'S.debilis' in GenBank, which creates a separate clade distributed in the West Indian Ocean and the West Pacific.Both clades are parapatric, and a geographic boundary occurs between 13 • S and 20 • S in the Indian Ocean.The taxonomic status of the 'S.liui' clade (e.g., species, subspecies, and species complex) needs further clarification through additional material for morphological studies and additional nuclear genes analyses.

Figure 1 .
Figure 1.Sampling locations of S. debilis in the Atlantic and Indian Oceans, along with their basinscale grouping.The symbols on the map indicate the type of data used (references are on the legend).

2. 2 .
DNA Extraction, Amplification, and Sequencing DNA was extracted either from the fifth pair of the pleopods or from the pleonic muscle tissue using the IG-Spin™ DNA Prep 200 kit for DNA extraction following the manufacturer's protocol.The extracted DNA was used as a matrix for the amplification of the mitochondrial cytochrome c oxidase subunit I gene fragment I (COI), and the nuclear gene of the first internal transcribed spacer (ITS1).PCR amplification of the COI gene fragment was accomplished with the universal primers LCOI 1490 (GGTCAACAAATCATAAAGATATTGG)

Figure 1 .
Figure 1.Sampling locations of S. debilis in the Atlantic and Indian Oceans, along with their basinscale grouping.The symbols on the map indicate the type of data used (references are on the legend).

Figure 3 .Figure 3 .
Figure 3. (A) Bayesian consensus phylogram of S. debilis based on mitochondrial cytochrome c oxidase I (COI) gene fragment (539bp).The horizontal scale bar marks the number of expected substitutions per site.Statistical support indicated as Bayesian posterior probabilities (left) and Maximum Likelihood bootstrap values for 1000 pseudoreplicates (right).(B) Minimum-spanning networks of S. debilis COI gene fragment.The size of the filled circles represents the number of indi-Figure 3. (A) Bayesian consensus phylogram of S. debilis based on mitochondrial cytochrome c oxidase I (COI) gene fragment (539 bp).The horizontal scale bar marks the number of expected substitutions per site.Statistical support indicated as Bayesian posterior probabilities (left) and Maximum Likelihood bootstrap values for 1000 pseudoreplicates (right).(B) Minimum-spanning networks of S. debilis COI gene fragment.The size of the filled circles represents the number of individuals with each haplotype, with the smallest circles representing one individual with that haplotype, color represents sampling regions.Hatch marks on the branches represent the number of mutational steps.

Figure 4 .
Figure 4.The RDA results.(A).Ordination of morphological characteristics in the space of canonical axes RDA1 and RDA2.(B).Variability of individuals of S. debilis from the Indian Ocean, Northern, and Southern Atlantic along the non-canonical PC1 axis.See characteristic coding in Table1.

Figure 4 .
Figure 4.The RDA results.(A).Ordination of morphological characteristics in the space of canonical axes RDA1 and RDA2.(B).Variability of individuals of S. debilis from the Indian Ocean, Northern, and Southern Atlantic along the non-canonical PC1 axis.See characteristic coding in Table1.

Diversity 2023, 15 , 1008 10 of 17 Figure 5 .
Figure 5. Variation in morphological characteristics in populations of S. debilis from the Indian Ocean (21 specimens), Northern Atlantic (43), and Southern Atlantic (5): (A) Number of spines in the posterior row on the merus of the fourth pereopods; (B) Number of spines in the posterior row on the merus of the fourth pereopods, (B).Number of spines in the posterior row on the merus of the third pereopods, (C).Number of spines in the posterior row on the merus of the fourth pereopods, (D).Serrations on the lateral margin of the fourth abdominal somite, right side.

Figure 5 .
Figure 5. Variation in morphological characteristics in populations of S. debilis from the Indian Ocean (21 specimens), Northern Atlantic (43), and Southern Atlantic (5): (A) Number of spines in the posterior row on the merus of the fourth pereopods; (B) Number of spines in the posterior row on the merus of the fourth pereopods, (B).Number of spines in the posterior row on the merus of the third pereopods, (C).Number of spines in the posterior row on the merus of the fourth pereopods, (D).Serrations on the lateral margin of the fourth abdominal somite, right side.

Table 1 .
Morphological characteristics and their minimum, maximum, and average values with standard deviation (SD) for 73 S. debilis specimens collected in three geographical regions (North Atlantic, 44; South Atlantic, 5: Indian Ocean, 24 specimens).In qualitative characteristics "−" indicates absence, and "+" indicates presence of the morphological character.

Table 3 .
Partitioning of variance in morphological characteristics of S. debilis based on RDA results.The constrained axes correspond with the body size, and the unconstrained describe the characteristics that do not correlate with the size.