Next Article in Journal
The Sum of Two Halves May Be Different from the Whole—Effects of Splitting Sequencing Samples Across Lanes
Next Article in Special Issue
The Complete Mitochondrial Genome of Dendrogale murina (Tupaiidae) and Phylogeny of Scandentia
Previous Article in Journal
Characterizing Macrophages Diversity in COVID-19 Patients Using Deep Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genomic Analyses Implicate the Amazon–Orinoco Plume as the Driver of Cryptic Speciation in a Swimming Crab

1
Department of Biology, Institute of Environment, Florida International University (FIU), Miami, FL 33199, USA
2
Laboratory of Bioecology and Systematics of Crustaceans (LBSC), Faculty of Philosophy, Sciences and Letters at Ribeirão Preto (FFCLRP), University of São Paulo (USP), Ribeirão Preto 14040-901, Brazil
3
Department of Invertebrate Zoology, National Museum of Natural History-Smithsonian, Washington, WA 20013-7012, USA
4
College of Fisheries and Ocean Sciences, University of Alaska Fairbanks, Fairbanks, AK 99775, USA
5
Auke Bay Laboratories, Alaska Fisheries Science Center, NOAA National Marine Fisheries Service, Juneau, AK 99801, USA
*
Author to whom correspondence should be addressed.
Genes 2022, 13(12), 2263; https://doi.org/10.3390/genes13122263
Submission received: 14 October 2022 / Revised: 25 November 2022 / Accepted: 25 November 2022 / Published: 1 December 2022
(This article belongs to the Special Issue Phylogenomics and Molecular Evolution)

Abstract

:
The Amazon–Orinoco plume (AOP) is the world’s largest freshwater and sediment discharge into the ocean. Previous studies limited to mtDNA suggest that the swimming crab Callinectes ornatus Ordway, 1863 exists as two distinct genetic clusters separated by the AOP. However, questions concerning migration, diversification time, and species delimitation are unresolved. Densely sampling markers across the genome (SNPs) could elucidate the evolutionary processes within this species. Here, we combined mtDNA data and ddRAD-seq to explore the diversification patterns and processes within the swimming crab C. ornatus. We show great genetic differentiation between groups on the north and south sides of the plume but also signs of hybridization. Demographic modeling indicates the divergence between groups starting around 8 Mya following the AOP’s formation. After a period of isolation, we detect two incidences of secondary contact with stronger migration in concordance with the North Brazil Current flow. Our results suggest speciation with gene flow explained by the interplay among the AOP, oceanographic currents, and long larval dispersal. This work represents the first investigation employing ddRAD-seq in a marine invertebrate species with distribution encompassing the north and south Atlantic and sheds light on the role of the AOP in the diversification of a marine species.

1. Introduction

Understanding marine species diversification remains challenging for ecologists and evolutionary biologists [1]. The interplay among complex life cycles (including multiple larval stages), pelagic larval duration (PLD), and oceanographic processes results in convoluted diversity patterns in marine systems (e.g., [2,3]). Moreover, processes of diversification in the marine system can be easily overlooked because many marine barriers are not conspicuous (e.g., mountains in a terrestrial environment) [4]. Most marine barriers are “soft”, that is, their permeability depends on the species, representing physical isolation for some species or environmental gradients for others [4]. Therefore, it is hard to predict the effects of marine barriers on diversification at the intra- or interspecific level. There is evidence for low or high genetic structure, allopatric and sympatric speciation, speciation with gene flow, and cryptic speciation scenarios [3,5,6,7,8,9,10,11]. However, the mechanisms responsible for each response are usually unknown and local/taxa-dependent. By investigating the effects of soft barriers on marine taxa, advancements can be made toward understanding marine diversification.
One of the most notable marine barriers is the Amazon–Orinoco plume (AOP) [12]. Located in the tropical western Atlantic, the AOP is delimited by the Amazon and Orinoco river mouths in the northern part of South America and represents the world’s largest freshwater and sediment discharge into the ocean [12]. The AOP’s water outflow is 6300 km3/year and extends 400 km off the coast and 30 m in depth [13]. This represents an environmental barrier known to change water salinity (36 to 32 psu) and turbidity (silicate from 1 to 9 μmol/kg) [13,14]. The AOP originated approximately 8–10 million years ago (Mya) following the Andes uplift and continued to develop by intensifying river outflow and sedimentation until 2–5 Mya when it was fully established [15]. In some species, results from few molecular markers show evidence for genetic structure and sister species separated by the AOP; however, some other specific marine groups are not affected by this barrier (annelids—[16]; crustacean—[17]; fish—[18,19,20]; manatees—[21]; mussel—[22]). The factors driving each scenario are unclear, but we can hypothesize that low tolerance to freshwater and/or turbid water may be driving diversification patterns. For instance, there is evidence for a widespread western Atlantic coral species showing genetic groups separated by the AOP [23] and evidence of coral larvae being affected by salinity changes [24]. Moreover, the comparison of two swimming crab species with known distinct physiological salinity tolerances shows that they are affected in different ways by the AOP [25]. Unfortunately, salinity tolerance data are not available for many species, and the effect of the AOP upon diversification is assumed a posteriori.
When investigating genetic structure in marine organisms, we must also consider the effect of pelagic larval duration (PLD). Many marine species show indirect development with larval stages during their life cycle. The larvae are released in the water and have the potential to be transported by ocean currents and reach distant locations [26]. There is a positive association between the PLD and the distance traveled by the larvae [27]. Consequently, longer PLD would theoretically result in a lack of population differentiation (e.g., [8]). However, across marine species, the association between PLD and population differentiation has been shown to range from weak to null [28]. Indeed, many species show some level of genetic structure across their distribution [3,6].
Marine populations or species separated by the AOP represent promising candidates for investigations on the interplay between PLD and marine barriers because they can help us to understand if one factor can overpower the other. Studies analyzing species occurring on both sides of the AOP have only used a few markers (Sanger sequencing)—usually mitochondrial DNA (mtDNA) (e.g., [16,17,18,19,20,22]). There are examples of mtDNA not capturing structure found by loci across the nuclear genome (single nucleotide polymorphisms—SNPs; [29]), as well as cases in which mtDNA delineates several cryptic species that are not supported by SNP data [30], an event called mitonuclear discordance [31]. Thus, studies employing genomic methods examining more loci across the nuclear genome (single nucleotide polymorphisms—SNPs) are warranted (but see [29] and [32]) because genetic structure and speciation can be overlooked in species occurring on both sides of the AOP [33,34]. Some topics that demand further investigation are the divergence time between genetic groups separated by the AOP estimated by hundreds of loci, comparisons between results coming from mtDNA (Sanger sequencing) and genomic approaches, and differentiation across the genome.
The swimming crab Callinectes ornatus Ordway, 1863 represents the ideal model organism to investigate the effects of a soft barrier (AOP) on the diversification of marine species in the tropical western Atlantic. The species is widespread along the western Atlantic, occurring in coastal waters (up to 75 m depth) from south Brazil to North Carolina, USA [35,36]. Despite nine larval stages and a PLD of 50 days [37,38,39], the species shows a strong genetic structure, composed of two separated groups: one north and one south of the AOP [25]. Unlike many species for which physiological data are not available, there is evidence showing that C. ornatus is sensitive to salinity changes [40]. Considering this a priori data showing that C. ornatus is not tolerant to low salinities, the authors hypothesized that salinity tolerance plays a significant role in shaping these two groups [25]. Unfortunately, due to the limited resolution of the data (mtDNA—cytochrome oxidase I (COI) and 16S ribosomal RNA (16S rRNA)), questions of demographic and diversification processes remain to be addressed. For instance, there is no clear resolution if two groups represent different species, and the results might represent only the mtDNA evolutionary history and not the species tree [31]. Therefore, using more loci can improve our characterization of these lineages, shed light on the role of the AOP upon western Atlantic species, and determine if C. ornatus represents a case of cryptic speciation.
Here, we combined mtDNA (COI gene) and ddRAD-seq to explore the genetic structure and demographic history of C. ornatus, focusing on the role of the AOP in the diversification within this species. We quantified genetic differentiation and investigated demographic processes between groups on both sides of the AOP (north and south tropical western Atlantic). Considering the scenario shown by mtDNA [25], we hypothesize that we will find evidence of speciation driven by the AOP when examining SNPs from across the genome and evidence that shared mtDNA haplotypes are the result of ancient polymorphisms. To the best of our knowledge, our work represents the first investigation employing ddRAD-seq in a tropical marine invertebrate species with distribution encompassing both sides of the AOP. We provide empirical evidence of a dynamic demographic history indicating that the AOP has led to diversification in a species displaying long PLD but low salinity tolerance.

2. Materials and Methods

2.1. Sampling and DNA Extraction for mtDNA and ddRAD-Seq

We used the same cytochrome c oxidase subunit I (COI) dataset for the mtDNA analyses as in [25]. We expanded it with sequences available on GenBank, plus de novo sequences generated and submitted to GenBank from regions not included in previous analyses (northern Brazil, French Guiana (region within the plume), North Carolina (USA), n = 12) (Table S1). For the ddRAD-seq analyses, we obtained 63 individuals of C. ornatus (Figure 1) from the following collections: Crustacean Collection of the Department of Biology (CCDB)—Faculty of Philosophy, Sciences, and Letters at Ribeirão Preto (FFCLRP) of the University of São Paulo (USP); the Invertebrate Zoology Collection—Florida Museum of Natural History (FLMNH) of the University of Florida; and the Florida International University Crustacean Collection (FICC). Many of them were used for both mtDNA and ddRAD analyses, which enabled us to compare the results from both types of markers. Our sampling covers the species’ described range (Figure 2A and Table S1).
Genomic DNA was extracted from muscle tissues using the salt extraction method [41,42] or with the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany), following the protocol provided by the manufacturer.

2.2. COI Amplification and Analyses

COI sequences were amplified using the primer COL6b (5′-ACAAATCATAAAGATATYGG-3′)/COH6 (5′-TADACTTCDGGRTGDCCAAARAAYCA-3′) [43], following the PCR cycle: initial denaturing for 5min at 94–95 °C; annealing for 35–40 cycles: 45 s at 95 °C, 45 s at 42–48 °C, 1 min at 72 °C; final extension of 5 min at 72 °C. PCR products were purified using a SureClean Plus kit (Bioline Reagents, UK) following the protocol provided by the manufacturer, and sequenced in an ABI 3730 xl DNA Analyzer (Applied Biosystems, Waltham, MA, USA) following Applied Biosystems protocols. Primer removal and quality checks were performed in Geneious Prime 2020.2.4 (https://www.geneious.com). Pseudogenes were identified by translating the consensus sequences and checking for indels and stop codons, removing them when present [44]. We aligned the sequences using MAFFT v.7 [45] and prepared a haplotype file in DnaSP v.6 [46], which was used to access the relationship among haplotypes using a statistical parsimony network through the TCS method [47] implemented in PopArt v.1.7 [48].

2.3. ddRAD-Seq Library Preparation and Data Processing

Double digest RADseq libraries were prepared according to the ddRADseq method [49]. Briefly, after enzyme trials to determine the best enzyme set, DNA from all individuals was digested with a combination of NlaIII and NotI (New England Biolabs, Ipswich, MA, USA). Following digestion, custom barcoded adapters were ligated to the fragments and pooled into twelve sublibraries. Each sublibrary was size selected (250–300 bp) on a PippinPrep with a 2% Agarose Gel Cassette (SageScience, Beverly, MA, USA). Size-selected fragments were then amplified via PCR with Phusion Hi-Fidelity Polymerase (Thermo Scientific, Waltham, MA, USA), which also incorporated indices (i7) and Illumina adapters into the fragments, allowing for the pooling of sublibraries into the final library. A washing step using AMPure XP beads (Beckman Coulter, Brea, CA, USA) was performed after digestion, adapter ligation, and PCR. Sequencing was done at the Genewiz Facility in South Plainfield, New Jersey, by an Illumina HiSeq 4000 (paired-end [PE] 150).
Raw sequence files were demultiplexed (--inline-index), cleaned (-c), quality-filtered (q), and rescued (-r), allowing two mismatches using the process_radtags program in STACKS v.2.3d [50]. Short reads were indexed and aligned to the blue crab (Callinectes sapidus) reference genome (Csap_IMET_V1, [51]) using Burrow–Wheller Aligner (BWA-MEM) v. 0.7.17 [52,53,54] with default settings. Generated SAM files were converted to BAM files using SAMtools v.1.9 [55]. Using BAM files as input, we ran the ref_map.pl wrapper program on STACKS v.2.3.d to call SNPs. We applied an iterative filtering strategy to minimize the effects of missing data and low-quality sequencing when calling SNPs. Using STACKS v.2.3d, we first ran the populations module to detect genotype call rate per locus (i.e., the proportion of individuals a locus is called in a given population; -r) using no threshold. Then, we used the function --missing-indv on VCFtools v. 0.1.17 [56] to assess missing data per individual and excluded individuals with >95% missing data. A second round excluded loci with a call rate <25% (-r 0.25) and individuals with >85% missing data. A third round was used to retrieve all loci with a genotype call rate of at least 50% (-r 0.5). This strategy is effective to exclude low-quality loci and individuals [57]. We set a minor allele frequency of 5% (--min-maf 0.05) because low-frequency alleles can affect population structure inferences [58]. The population module was also used to filter one SNP per rad-loci (--write-random-snp).

2.4. Linkage-Disequilibrium and Outlier Detection

Because most of the analyses assume that SNPs are independent, SNPs potentially under linkage disequilibrium (LD) were detected by calculating the squared correlation coefficient (r2) between SNP pairs using the function --geno-r2 in VCFtools v. 0.1.17. Only pairwise SNP comparisons within the same chromosome and with more than half of our final sample size were considered. SNPs showing r2 > 0.8 were excluded from the final dataset. Two different approaches were implemented to identify outlier SNPs (i.e., non-neutral SNPs, possibly under selection), BayeScan v.2.1 [59] and PCAdapt [60], as the use of multiple methods has been recommended to reduce type 1 errors [61,62,63]. BayeScan is a method to identify putative adaptive SNPs based on different allele frequencies among populations, and we performed it using default settings (prior odds to 10, iterations to 5000, and burn-in to 50,000). Outlier SNPs were identified at a q-value (i.e., false discovery rate) of 0.01. PCAdapt implements a hierarchical method (not assuming groups a priori) based on principal component analysis that identifies SNPs excessively related to population structure, probably due to selection. We ran PCAdapt exploring twenty PCs (K = 20) to select the optimal K following Cattell’s rule to retain the best K PC value [60]. These PCs detect the SNPs most associated with population structure. We detected putative non-neutral SNPs based on a q-value of 0.01. An SNP was considered non-neutral when both analyses indicated that the SNP is potentially non-neutral.

2.5. Genetic Diversity, Population Structure, and Genetic Differentiation

Summary genetic diversity statistics were calculated in GENODIVE v3.0 [64]. We employed the Bayesian program STRUCTURE v2.3.4 [65] to test for population structure within the data. Eight K-values were tested (K = 8) 10 times each under the admixture model. Following a burn-in of 10,000 generations, 100,000 Markov Chain Monte Carlo generations were run. The most optimal K value was estimated with the posterior probability for each model using the ΔK method [65,66], as implemented in STRUCTURE HARVESTER v0.6.94 [67]. We calculated corrected pairwise FST in GENODIVE v3.0 between the groups in the north and south of the AOP. We also accessed genetic differentiation across the genome by calculating kernel-smoothed FST considering a sliding window of 900 Kbp in total length using the population module in STACKS v.2.3d.

2.6. Demographic History: Divergence Time, Migration, and Effective Population Size

We performed demographic modeling analyses to find the best scenario explaining the two genetic groups—north and south, separated by the AOP (see results). Simulated and observed two-dimensional joint site frequency spectra (2D-JSFS) were compared using the Genetic Algorithm for Demographic Model Analysis (GADMA, [68]). Loci were considered not linked. Populations were projected down to smaller sample sizes (projection parameter set to north: 9 and south: 44) for 2D-JSFS estimation. Ordinary differential equations implemented in the moments engine [69] were used to simulate 2D-JSFS within GADMA. The Genetic Algorithm (GA) was used for global optimization and Powell’s conjugate method was used for local optimization during the model search in the parameter space. Optimizations were run for 50 repeats for each model. We set effective population size changes and migration to “true”, no demographic changes before the first split were allowed, and two demographic changes after the first split were allowed (initial structure: [1,1]; final structure [1,2]). All other parameters were set to default. Models were compared using Akaike’s Information Criterion (AIC). Parameters are given in genetic units in relation to reference effective population size (Nref). Nref was determined using θ/θ0, in which θ is the population-scaled mutation rate given in the output of the GADMA run and θ0 = 4 × µ (mutation rate per base per generation) × L (length of the sequence that was used to build the data). Because there is no µ estimate for C. ornatus or swimming crabs, we used the value calculated for Drosophila melanogaster [70] and Daphnia pulex [71]: 4 × 10−9. The average sequence length of the ddRAD-loci used as L = 188.41 bp. Callinectes ornatus generation time was considered as one year [72].

3. Results

3.1. mtDNA (COI)

We used 66 sequences (570 bp), including 12 new sequences not included in [25], representing North Carolina, northern Brazil, and French Guiana (within the AOP). The haplotype network depicts a clear split between the individuals from the north and south of the AOP, but one individual from the north group (FLMNH 32103) falls within the south network. Individuals from the AOP region show the most common haplotype, shared by the south group (Figure 2B).

3.2. ddRAD-Seq

After process_radtags (demultiplexing, cleaning, and quality-filtering), we retained ~570 million reads across 63 individuals (69% retained reads). After the ref_map pipeline, the unfiltered set of markers consisted of 91,098 loci, composed of 15,501,392 sites, covering approximately 2% of the genome. After the filtering steps, 10 individuals were excluded (16% of our total) and we kept 53 individuals (Table S1). The ref_map pipeline returned 1870 loci (mean = 188.41 bp) and 912 SNPs (one SNP per rad-locus).

3.3. Linkage-Disequilibrium and Outlier Detection

Thirteen SNPs showed signals of LD and were excluded from the dataset. Although BayeScan and PCAdapt detected potentially non-neutral SNPs, the analyses did not converge to the same SNP set. Therefore, we decided not to exclude any SNPs and assume differences in allele frequency and genetic differentiation in some loci due to neutral processes (e.g., genetic drift). Our final dataset for downstream analyses was composed of 899 SNPs.

3.4. Genetic Diversity, Population Structure, and Genetic Differentiation

Summary statistics of genetic diversity for the north and south groups are similar and are shown in Table 1. STRUCTURE results show K = 2 as the most probable number of genetic clusters within our dataset. The two groups represent a group on the north and one on the south of the AOP (Figure 2C). Some individuals are totally assigned to the north or south group, but the STRUCTURE analysis also detected admixed individuals (FLMNH 34910, FLMNH 32103, FLMNH 3982, FLMNH 1409b, CCDB 6105, CCDB 6130a, CCDB 6130f, CCDB 1537b, CCDB 353). Admixed individuals are not restricted to a specific geographical location. When we combine the results coming from ddRAD data and mtDNA haplotypes, we also see evidence of mitonuclear discordance. For instance, in the STRUCTURE plot, an individual from Saint Martin (FLMNH 32103) has the south mtDNA, but it has the genetic signatures of both the south and north. In contrast, individuals from Florida (FLMNH 26242, FLMNH 1476) are assigned to the north group based on mtDNA but are assigned to the south cluster under the ddRAD analysis (Figure 2B and 2C). The pairwise FST indicates high divergence between genetic groups detected by STRUCTURE (FST = 0.267). FST across the genome indicates moderate to high differentiation in all chromosomes (Figure 2D).

3.5. Divergence Time, Migration, and Effective Population Size

GADMA indicates that the most probable demographic scenario is allopatric divergence followed by secondary contact at different times (log-likelihood = −260.37; AIC score = 548.74) (Figure 3). The split between the north and south groups occurred at approximately ~8 Mya (Genetic Units = 5.09). The north group shows a linear population size increase until ~150 thousand years ago (Kya) (Genetic Units = 0.09) when it shows a decrease. The south group shows a constant population size until 150 Kya when it shows signals of population expansion. After a period of no contact between the two genetic groups, GADMA detected migration from north to south at approximately ~4 Mya and bidirectional migration starting at ~150 Kya. Migration from south to north is higher than from north to south, indicating asymmetrical migration.

4. Discussion

Our results indicate allopatric divergence following the emergence of the Amazon–Orinoco plume (AOP) between C. ornatus on the north and south sides of the AOP. However, we also found evidence for secondary contact, the first occurring during the Miocene (~4 Mya) and characterized by migration from the south to the north, and the second occurring during the Pleistocene and characterized by asymmetrical migration, which causes admixture between both lineages. These results challenge the validity of C. ornatus as a single species and suggest a cryptic speciation scenario with gene flow. Employing mtDNA and ddRAD-seq, we provide an example of how dynamic the diversification process can be in the marine realm and stress the importance of the AOP for the diversification of marine species within this region.

4.1. The Role of the Amazon–Orinoco Plume (AOP) in the Diversification of Callinectes ornatus

We found high genetic differentiation between groups separated by the AOP (north and south) in both types of markers (mtDNA and SNPs). Here, we added to [25] COI dataset sequences from three other locations: the north coast of Brazil, the northern coast of South America (within the plume), and the North Carolina samples. As expected, the North Carolina mtDNA haplotypes fall within the north network, but sequences from within the plume are represented by the most common haplotype in the south network. This pattern can be explained by the North Brazil Current (NBC—[73,74]) transporting individuals in the northward direction (south to north—see Figure 2A). Uniparental inheritance, smaller effective population size, and high mutation rate in comparison to nuDNA [75,76] in combination with the AOP acting as a barrier can explain the deep split we see between north and south despite continuous distribution. When looking at ddRAD markers, pairwise FST and STRUCTURE also find two genetic groups (north and south), showing high genetic differentiation (FST = 0.26). However, the STRUCTURE plot also shows admixture between north and south, indicating that complete differentiation has not occurred, as expected in cases of complete speciation (e.g., [77,78]). We find individuals being assigned to the opposite cluster, despite their geographic location and individuals showing signs of mixture and mitonuclear discordance. Admixed individuals are found in Florida and Saint Martin and also in the northeastern and southeastern coast of Brazil; mitonuclear discordance was only found in individuals on the north side of the AOP.
Considering that there are sister species separated by the AOP in a broad range of taxa (e.g., [22]—mussel, [79]—fish), our scenario could be explained by complete isolation but not complete divergence due to large effective population sizes reducing genetic drift [80]. We modeled the most probable demographic scenario to gain insight into the mechanisms that may lead to our results. However, our most probable scenario indicates divergence followed by secondary contact at different times. The estimated divergence time between the north and south groups is approximately 8 Mya, but then we detect a weak signal of migration from south to north around 4 Mya and a following phase of asymmetric migration. The geological formation of the AOP is characterized by different events, with the initial formation at 8–10 Mya [15]. A meta-analysis investigating the divergence time of reef-fish species along the Western Atlantic shows that this is the moment in which many species started to diverge, reinforcing the role of the AOP in the diversification of marine species from this region [20]. It is still unclear which traits play a role in how species respond to the barrier, but data on fish show that larger species seem to be less affected by the barrier [20]. This is probably because size seems to be associated with the ability to cross soft barriers, long-distance migration, and salinity and turbid water resistance [81]. Based on the species distribution data of all major animal groups on the north and south sides of the AOP, low salinity is evoked as one of the major factors contributing to the hampering dispersal [82]. In our case, salinity tolerance can be the trait responsible for species diverging at this initial phase of the AOP formation considering available data on C. ornatus’s sensitivity to salinity changes. A comparative phylogeographic study investigating C. ornatus and C. danae, two swimming crabs with the same dispersal potential but different salinity tolerance, indicates that the tolerant species (C. danae) is not affected by the plume [25]. Therefore, the onset of the AOP might be the kickstart of the divergence between C. ornatus groups on the north and south sides of the plume.
GADMA detected low migration rates from south to north around 4 Mya, which can be explained by the closure of the Isthmus of Panama, leading to the formation of the North Brazil Current (NBC) [83]. However, this migration ceased, probably because of the progression of the Andean uplift, which resulted in increased freshwater and sediment outflow around the same time [15]. The isolation between the two groups onwards can also be explained by the sea level changes throughout the Pleistocene, causing decreases in coastal habitats and closing deep reefs below the plume that could be used as steppingstones [84]. However, our demographic model also indicates another signal of migration around 150 Kya, characterized by an asymmetric migration (south to north > north to south). Secondary contact after a period of isolation can also be explained by the changing cycles during the Pleistocene [84]. For instance, although coastal areas decreased during sea level falls, they also increased during the next moment of sea level rise [84]. Our demographic model shows population expansion of the north around 150 Kya, which might have increased the chances of contact between the north and south. This secondary contact during this time can also be explained by changes that occurred during the Pleistocene, such as changes in marine currents’ dynamic [15,85,86,87]. The NBC retroflection (i.e., the event when the NBC flows southwards and west) is a common phenomenon from the Pleistocene until the present, which might intensify the transport of larvae out of the AOP influence and make the contact between south and north more frequent [88]. Therefore, our results indicate that the AOP represents a barrier for C. ornatus, but variations in the plume over time might open opportunities for larvae to be transported between regions, especially from south to north, because of the NBC flowing northwards. Even with two moments of secondary contact (4 Mya and present), the allopatric periods between north and south might have been enough to accumulate divergence caused by genetic drift. This result is corroborated by the fact that we found regions showing high FST across the whole genome and not concentrated in specific regions, indicating that genome-wide divergence has been accumulating for a long time [89]. At the same time, a long PLD coupled with the plume and current dynamics might result in occasional gene flow between north and south groups. The diversification process in marine species should be seen not only considering potential barriers but also the likelihood of larvae overcoming these barriers over time. Therefore, multiple scenarios are possible, specifically in marine systems [90].

4.2. Speciation with Gene Flow in Callinectes ornatus

A traditional paradigm postulates marine environments as open systems because of the lack of barriers and potential broad dispersal due to an extended larval stage of several marine species [91]. However, empirical studies reveal diversification occurring without complete isolation and cryptic speciation as a common phenomenon at sea [5,92]. In fact, diversification is characterized by not only strict isolation but also isolation with migration, ancient migration, secondary contact, or periodic connectivity [90]. The interplay between semi-permeable barriers and dispersal potential results in speciation with gene flow as a likely scenario in marine species [11], which also seems to be the case for C. ornatus.
Even isolated for millions of years, lineages can interbreed after secondary contact if no reproduction isolation has evolved [89,93]. As we discussed, larvae might occasionally travel across the AOP and promote hybridization between groups. Potential hybrids (i.e., individuals showing admixture) were found in Florida, Saint Martin, and Brazil, indicating no hybrid zone or unique region of contact between the genetic groups. Higher migration in the northward direction might explain the occurrence of mitonuclear discordance only on the north side of the AOP. The system might be characterized as a population–species continuum throughout the species distribution [94,95]. This is in accordance with the fact that adults are not migratory, but larvae can reach further regions. Based on genetic and distributional data, previous studies showed that the connection between both sides of the AOP is always in the northward direction [20,82]. Using more regions across the genome and modeling the demographic history, our results expand this current vision and provide evidence of asymmetric migration. The occurrence of admixtured individuals and mitonuclear discordance on both sides of the plume indicate that hybridization has been occurring since the south and north came in contact again, reinforcing connectivity in both directions.
Speciation is often a continuous process that eventually results in complete reproductive isolation [96]. In marine systems, we are finding speciation with gene flow increasingly common [11,97]. Under a general model of speciation with gene flow, the first phase represents positive selection on a few genes, while most of the genome shows low differentiation [89]. Later, divergence hitchhiking creates “islands of differentiation,” and, eventually, genome hitchhiking leads to great differentiation across the whole genome [89,98]. In our case, genomic differentiation across the genome shows multiple highly divergent regions and can be considered an early/intermediate moment in the speciation process [89]. We could not detect any regions under selection, which can be explained by small sample sizes or the proportion of the genome sampled, but it can also represent that genetic drift plays a role in accumulating differences between genomes over a long period of no contact. Therefore, not only does the mtDNA show levels of differentiation compatible with different species (>4%, [25]), but the whole genome seems to show the same pattern, indicating undocumented speciation within C. ornatus.

5. Conclusions

Many questions remain concerning the diversification patterns and processes of marine species [1,97]. Our work contributes to the field by employing mtDNA and ddRAD-seq to investigate the interplay between PLD and a marine barrier (AOP). Moreover, our work expands current knowledge by having an arthropod as a model since three-quarters of marine diversification studies are on Chordata, Mollusca, and Cnidaria [90]. We find evidence that the swimming crab C. ornatus has gone through a divergence period likely driven by the emergence of the AOP, but the two lineages mixed after secondary contact caused by larvae transported by marine currents. Low salinity tolerance is likely the trait responsible for differentiation, and PLD is responsible for secondary contact. A thorough morphological investigation has failed to reveal a diagnostic character that allows species-level differentiation ([99], personal communication). Here, we take a conservative approach to not nominating a new species due to the lack of operational criteria [100]. However, our findings provide evidence for a case of speciation with gene flow, revealing the presence of two separately evolving lineages [100].

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13122263/s1, Table S1: Individuals of Callinectes ornatus analyzed in COI and ddRAD analyses. Vouch. Coll.: Voucher Collection ID. COI and GB COI: indicates if the specimen was used in COI analyses with the accession number. ddRAD: indicates all individuals used (total) and included in the analyses after filtering steps (final dataset). Locality indicates where they were sampled. Latitude and longitude are shown as indicated on the original Voucher Collection ID tag or GB; otherwise, we indicated them as not available (n/a). CCDB = Crustacean Collection of Department of Biology, FFCLRP/USP, Brazil; FLMNH = Florida Museum of Natural; HBG = Florida International Crustacean Collection; ULLZ = University of Louisiana at Lafayette Zoological Collection; MNHN = Muséum National D’Histoire Naturelle; USNM = National Museum of Natural History.

Author Contributions

Conceptualization, P.A.P. and F.L.M.; methodology, P.A.P., H.B.-G. and L.E.T.; formal analysis, P.A.P., H.B.-G. and L.E.T.; investigation, P.A.P.; resources, F.L.M. and H.B.-G.; data curation, P.A.P.; writing—original draft preparation, P.A.P.; writing—review and editing, P.A.P., H.B.-G., L.E.T. and F.L.M.; supervision, F.L.M. and H.B.-G.; visualization, P.A.P.; data curation, P.A.P. and F.L.M.; project administration, F.L.M. and H.B.-G.; funding acquisition, F.L.M. and H.B.-G. All authors have read and agreed to the published version of the manuscript.

Funding

This research resulted from the PhD dissertation of Pedro A. Peres, supported by the São Paulo Research Foundation—FAPESP, grant number 2017/12376-6 and 2017/12376-6. All financial support was provided by grants from FAPESP Temático Biota, 2010/50188-8, and Biota INTERCRUSTA, 2018/13685-5; Coleções Científicas, 2009/54931-0; PROTAX, 2016/50376-5, Conselho Nacional de Desenvolvimento Científico e Tecnológico—CNPq, 301359/2007-5; 473050/2007-2; 302748/2010-5; 471011/2011-8; PQ 302253/2019-0; and Coordenação de Aperfeiçoamento de Pessoal—CAPES—Código de Financiamento 001, Ciências do Mar II—2005/2014–23038.004308/201414, granted to Fernando L. Mantelatto. The work completed in the Bracken-Grissom lab was partially supported by the National Science Foundation Division of Biology, grant number DEB-1856667, awarded to Heather Bracken-Grissom at Florida International University.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Demultiplexed sequence reads (forward and reverse) for every specimen are deposited in SRA (BioProject PRJNA901448). The COI haplotype data generated and used in this study are deposited in NCBI Nucleotide Database (Table S1).

Acknowledgments

We are grateful to Gustav Paulay and John Slapcinsky (FLMNH), Cleverson Santos (MPEG), and all the colleagues of the Laboratory of Bioecology and Crustacean Systematics (LBSC) for their help with sampling, donation, and lending the specimens used herein. We are also grateful to Pedro A. Peres’ committee members, Fernando Franco, Fosca Pedini Pereira Leite, Mariana Terossi, Rafael Robles, and Sónia Cristina da Silva Andrade, for suggested improvements to the manuscript. All analyses were run on the Florida International University High-Performance Computing Cluster (HPCC). The collections of species conducted in this study complied with the current applicable state and federal laws of Brazil (permanent license to FLM for collection of Zoological Material No. 11777-2 MMA/IBAMA/ SISBIO). This is contribution #1511 from the Institute of Environment at Florida International University.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bowen, B.W.; Shanker, K.; Yasuda, N.; Celia, M.; Malay, M.C.M.D.; von der Heyden, S.; Paulay, G.; Rocha, L.A.; Selkoe, K.A.; Barber, P.H.; et al. Phylogeography unplugged: Comparative surveys in the genomic era. Bull. Mar. Sci. 2014, 90, 13–46. [Google Scholar] [CrossRef] [Green Version]
  2. Eble, J.A.; Toonen, R.J.; Bowen, B.W. Endemism and dispersal: Comparative phylogeography of three surgeon fishes across the Hawaiian Archipelago. Mar. Biol. 2009, 156, 689–698. [Google Scholar] [CrossRef]
  3. Kelly, R.P.; Palumbi, S.R. Genetic structure among 50 species of the northeastern Pacific rocky intertidal community. PLoS ONE 2010, 5, e8594. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Briggs, J.C.; Bowen, B.W. Marine shelf habitat: Biogeography and evolution. J. Biogeogr. 2013, 40, 1023–1035. [Google Scholar] [CrossRef]
  5. Knowlton, N. Sibling species in the sea. Annu. Rev. Ecol. Evol. Syst. 1993, 24, 189–216. [Google Scholar] [CrossRef]
  6. Pelc, R.A.; Warner, R.R.; Gaines, S.D. Geographical patterns of genetic structure in marine species with contrasting life histories. J. Biogeogr. 2009, 36, 1881–1890. [Google Scholar] [CrossRef]
  7. Gaither, M.R.; Toonen, R.J.; Robertson, D.R.; Planes, S.; Bowen, B.W. Genetic evaluation of marine biogeographical barriers: Perspectives from two widespread Indo-Pacific snappers (Lutjanus kasmira and Lutjanus fulvus). J. Biogeogr. 2009, 37, 133–147. [Google Scholar] [CrossRef]
  8. Reece, J.S.; Bowen, B.W.; Smith, D.G.; Larson, A. Comparative phylogeography of four Indo-Pacific moray eel species (Muraenidae) reveals comparable ocean-wide genetic connectivity despite five-fold differences in available adult habitat. Mar. Ecol. Prog. Ser. 2011, 437, 269–277. [Google Scholar] [CrossRef] [Green Version]
  9. Leasi, F.; Andrade, S.C.D.S.; Norenburg, J. At least some meiofaunal species are not everywhere. Indication of geographic, ecological and geological barriers affecting the dispersion of species of Ototyphlonemertes (Nemertea, Hoplonemertea). Mol. Ecol. 2016, 25, 1381–1397. [Google Scholar] [CrossRef]
  10. Álvarez-Campos, P.; Giribet, G.; Riesgo, A. The Syllis gracilis species complex: A molecular approach to a difficult taxonomic problem (Annelida, Syllidae). Mol. Phylogenet. Evol. 2017, 109, 138–150. [Google Scholar] [CrossRef]
  11. Potkamp, G.; Fransen, C.H. Speciation with gene flow in marine systems. Contrib. Zool. 2019, 88, 1–40. [Google Scholar] [CrossRef] [Green Version]
  12. Curtin, T.B. Physical observations in the plume region of the Amazon River during peak discharge—II. Water masses. Cont. Shelf. Res. 1986, 6, 53–71. [Google Scholar] [CrossRef]
  13. Pailler, K.; Bourlès, B.; Gouriou, Y. The barrier layer in the western tropical Atlantic Ocean. Geophys. Res. Lett. 1999, 26, 2069–2072. [Google Scholar] [CrossRef]
  14. Wright, L.D.; Nittrouer, C.A. Dispersal of river sediments in coastal seas: Six contrasting cases. Estuaries 1995, 18, 494–508. [Google Scholar] [CrossRef]
  15. Hoorn, C.; Wesselingh, F.P.; Ter Steege, H.; Bermudez, M.A.; Mora, A.; Sevink, J.; Sanmartín, I.; Sanchez-Meseguer, A.; Anderson, C.L.; Figueiredo, J.P.; et al. Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity. Science 2010, 330, 927–931. [Google Scholar] [CrossRef] [Green Version]
  16. Nunes, F.L.; Van Wormhoudt, A.; Faroni-Perez, L.; Fournier, J. Phylogeography of the reef-building polychaetes of the genus Phragmatopoma in the western Atlantic Region. J. Biogeogr. 2017, 44, 1612–1625. [Google Scholar] [CrossRef] [Green Version]
  17. Tourinho, J.L.; Solé-Cava, A.M.; Lazoski, C. Cryptic species within the commercially most important lobster in the tropical Atlantic, the spiny lobster Panulirus argus. Mar. Biol. 2012, 159, 1897–1906. [Google Scholar] [CrossRef]
  18. Joyeux, J.C.; Floeter, S.R.; Ferreira, C.E.L.; Gasparini, J.L. Biogeography of tropical reef fishes: The South Atlantic puzzle. J. Biogeogr. 2008, 28, 831–841. [Google Scholar] [CrossRef]
  19. Silva, D.; Martins, K.; Oliveira, J.; da Silva, R.; Sampaio, I.; Schneider, H.; Gomes, G. Genetic differentiation in populations of lane snapper (Lutjanus synagris—Lutjanidae) from Western Atlantic as revealed by multilocus analysis. Fish. Res. 2018, 198, 138–149. [Google Scholar] [CrossRef]
  20. Araujo, G.S.; Rocha, L.A.; Lastrucci, N.S.; Luiz, O.J.; Di Dario, F.; Floeter, S.R. The Amazon-Orinoco Barrier as a driver of reef-fish speciation in the Western Atlantic through time. J. Biogeogr. 2022, 49, 1407–1419. [Google Scholar] [CrossRef]
  21. Luna, F.D.O.; Beaver, C.E.; Nourisson, C.; Bonde, R.K.; Attademo, F.L.; Miranda, A.V.; Torres-Florez, J.P.; De Sousa, G.P.; Passavante, J.Z.; Hunter, M.E. Genetic connectivity of the West Indian manatee in the southern range and limited evidence of hybridization with Amazonian manatees. Front. Mar. Sci. 2021, 7, 1089. [Google Scholar]
  22. Trovant, B.; Basso, N.G.; Orensanz, J.M.; Lessa, E.P.; D’Incao, F.; Ruzzante, D.E. Scorched mussels (Brachidontes spp., Bivalvia: Mytilidae) from the tropical and warm-temperate southwestern Atlantic: The role of the Amazon River in their speciation. Ecol. Evol. 2016, 6, 1778–1798. [Google Scholar] [CrossRef] [Green Version]
  23. Souza, J.N.; Nunes, F.L.; Zilberberg, C.; Sanchez, J.A.; Migotto, A.E.; Hoeksema, B.W.; Serrano, X.M.; Baker, A.C.; Lindner, A. Contrasting patterns of connectivity among endemic and widespread fire coral species (Millepora spp.) in the tropical Southwestern Atlantic. Coral. Reefs 2017, 36, 701–716. [Google Scholar] [CrossRef]
  24. Vermeij, M.J.A.; Fogarty, N.D.; Miller, M.W. Pelagic conditions affect larval behavior, survival, and settlement patterns in the Caribbean coral Montastraea faveolata. Mar. Ecol. Progr. Ser. 2006, 310, 119–128. [Google Scholar]
  25. Peres, P.A.; Mantelatto, F.L. Salinity tolerance explains the contrasting phylogeographic patterns of two swimming crabs species along the tropical western Atlantic. Evol. Ecol. 2020, 34, 589–609. [Google Scholar] [CrossRef]
  26. Hedgecock, D. Is gene flow from pelagic larval dispersal important in the adaptation and evolution of marine invertebrates? Bull. Mar. 1986, 39, 550–564. [Google Scholar]
  27. Shanks, A.L. Pelagic larval duration and dispersal distance revisited. Biol. Bull. 2009, 216, 373–385. [Google Scholar] [CrossRef] [Green Version]
  28. Weersing, K.; Toonen, R.J. Population genetics, larval dispersal, and connectivity in marine systems. Mar. Ecol. Prog. Ser. 2009, 393, 1–12. [Google Scholar] [CrossRef] [Green Version]
  29. Pedraza-Marrón, C.D.R.; Silva, R.; Deeds, J.; Van Belleghem, S.M.; Mastretta-Yanes, A.; Domínguez-Domínguez, O.; Rivero-Vega, R.A.; Lutackas, L.; Murie, D.; Parkyn, D.; et al. Genomics overrules mitochondrial DNA, siding with morphology on a controversial case of species delimitation. Proc. R. Soc. B 2019, 286, 20182924. [Google Scholar] [CrossRef] [PubMed]
  30. Hinojosa, J.C.; Koubínová, D.; Szenteczki, M.A.; Pitteloud, C.; Dincă, V.; Alvarez, N.; Vila, R. A mirage of cryptic species: Genomics uncover striking mitonuclear discordance in the butterfly Thymelicus sylvestris. Mol. Ecol. 2019, 28, 3857–3868. [Google Scholar] [CrossRef]
  31. Toews, D.P.; Brelsford, A. The biogeography of mitochondrial and nuclear discordance in animals. Mol. Ecol. 2012, 21, 3907–3930. [Google Scholar] [CrossRef]
  32. Titus, B.M.; Blischak, P.D.; Daly, M. Genomic signatures of sympatric speciation with historical and contemporary gene flow in a tropical anthozoan (Hexacorallia: Actiniaria). Mol. Ecol. 2019, 28, 3572–3586. [Google Scholar] [CrossRef] [PubMed]
  33. Seehausen, O.; Butlin, R.K.; Keller, I.; Wagner, C.E.; Boughman, J.W.; Hohenlohe, P.A.; Peichel, C.L.; Saetre, G.P.; Bank, C.; Brännström, Å.; et al. Genomics and the origin of species. Nat. Rev. Gen. 2014, 15, 76–192. [Google Scholar] [CrossRef] [PubMed]
  34. Kulmuni, J.; Butlin, R.K.; Lucek, K.; Savolainen, V.; Westram, A.M. Towards the completion of speciation: The evolution of reproductive isolation beyond the first barriers. Philos. Trans. R. Soc. Lond. B Biol. Sci. 2020, 375, 20190528. [Google Scholar] [CrossRef] [PubMed]
  35. Norse, E.A. An experimental gradient analysis: Hyposalinity as an “upstress” distributional determinant for Caribbean portunid crabs. Biol. Bull. 1978, 155, 586–598. [Google Scholar] [CrossRef]
  36. Mantelatto, F.L.; Tamburus, A.F.; Magalhães, T.; Buranelli, R.C.; Terossi, M.; Negri, M.; Castilho, A.L.; Costa, R.C.; Zara, F.J. Checklist of decapod crustaceans from the coast of the São Paulo state (Brazil) supported by integrative molecular and morphological data: III. Infraorder Brachyura Latreille, 1802. Zootaxa 2020, 4872, 1–108. [Google Scholar] [CrossRef] [PubMed]
  37. Costlow, J.D., Jr.; Bookhout, C.G. The larval development of Callinectes sapidus Rathbun reared in the laboratory. Biol. Bull. 1959, 116, 373–396. [Google Scholar] [CrossRef]
  38. Bookhout, C.G.; Costlow Jr, J.D. Larval development of Callinectes similis reared in the laboratory. Bull. Mar. Sci. 1977, 27, 704–728. [Google Scholar]
  39. Mantelatto, F.L.; Reigada, A.L.; Gatti, A.C.; Cuesta, J.A. Morphology of the first zoeal stages of five species of the portunid genus Callinectes (Decapoda, Brachyura) hatched at the laboratory. An. Acad. Bras. Cienc. 2014, 86, 755–768. [Google Scholar] [CrossRef] [Green Version]
  40. Garçon, D.P.; Masui, D.C.; Mantelatto, F.L.; McNamara, J.C.; Furriel, R.P.M.; Leone, F.A. K+ and NH4+ modulate gill (Na+, K+)-ATPase activity in the blue crab, Callinectes ornatus: Fine tuning of ammonia excretion. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 2006, 147, 145–155. [Google Scholar] [CrossRef]
  41. Miller, S.; Dykes, D.; Polesky, H. A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Res. 1988, 16, 1215. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Mantelatto, F.L.; Robles, R.; Biagi, R.; Felder, D.L. Molecular analysis of the taxonomic and distributional status for the hermit crab genera Loxopagurus Forest, 1964 and Isocheles Stimpson, 1858 (Decapoda, Anomura, Diogenidae). Zoosystema 2006, 28, 495–506. [Google Scholar]
  43. Schubart, C.D.; Huber, M.G.J. Genetic comparisons of German populations of the stone crayfish, Austropotamobius torrentium (Crustacea: Astacidae). Bull. Fr. Pêche Piscic. 2006, 380–381, 1019–1028. [Google Scholar] [CrossRef] [Green Version]
  44. Song, H.; Buhay, J.E.; Whiting, M.F.; Crandall, K.A. Many species in one: DNA barcoding overestimates the number of species when nuclear mitochondrial pseudogenes are coamplified. Proc. Natl. Acad. Sci. USA 2008, 105, 13486–13491. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Katoh, K.; Standley, D.M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 2013, 30, 772–780. [Google Scholar] [CrossRef] [PubMed]
  46. Rozas, J.; Ferrer-Mata, A.; Sánchez-DelBarrio, J.C.; Guirao-Rico, S.; Librado, P.; Ramos-Onsins, S.E.; Sánchez-Gracia, A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol. Biol. Evol. 2017, 34, 3299–3302. [Google Scholar] [CrossRef]
  47. Clement, M.; Posada, D.C.K.A.; Crandall, K.A. TCS: A computer program to estimate gene genealogies. Mol. Ecol. 2000, 9, 1657–1659. [Google Scholar] [CrossRef] [Green Version]
  48. Leigh, J.W.; Bryant, D. popart: Full-feature software for haplotype network construction. Methods Ecol. Evol. 2015, 6, 1110–1116. [Google Scholar] [CrossRef]
  49. Peterson, B.K.; Weber, J.N.; Kay, E.H.; Fisher, H.S.; Hoekstra, H.E. Double digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS ONE 2012, 7, e37135. [Google Scholar] [CrossRef] [Green Version]
  50. Rochette, N.C.; Rivera-Colón, A.G.; Catchen, J.M. Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics. Mol. Ecol. 2019, 28, 4737–4754. [Google Scholar] [CrossRef]
  51. Bachvaroff, T.R.; McDonald, R.C.; Plough, L.V.; Chung, J.S. Chromosome-level genome assembly of the blue crab, Callinectes sapidus. G3 2021, 11, jkab212. [Google Scholar] [CrossRef] [PubMed]
  52. Li, H.; Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 2009, 25, 1754–1760. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Li, H.; Durbin, R. Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics 2010, 26, 589–595. [Google Scholar] [CrossRef] [Green Version]
  54. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv 2013, arXiv:1303.3997. [Google Scholar]
  55. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R. The sequence alignment/map format and SAMtools. Bioinformatics 2009, 25, 2078–2079. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Danecek, P.; Auton, A.; Abecasis, G.; Albers, C.A.; Banks, E.; DePristo, M.A.; Handsaker, R.E.; Lunter, G.; Marth, G.T.; Sherry, S.T.; et al. The variant call format and VCFtools. Bioinformatics 2011, 27, 2156–2158. [Google Scholar] [CrossRef]
  57. O’Leary, S.J.; Puritz, J.B.; Willis, S.C.; Hollenbeck, C.M.; Portnoy, D.S. These aren’t the loci you’re looking for: Principles of effective SNP filtering for molecular ecologists. Mol. Ecol. 2018, 27, 3193–3206. [Google Scholar] [CrossRef] [Green Version]
  58. Linck, E.; Battey, C.J. Minor allele frequency thresholds strongly affect population structure inference with genomic data sets. Mol. Ecol. Resour. 2019, 19, 639–647. [Google Scholar] [CrossRef]
  59. Foll, M.; Gaggiotti, O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics 2008, 180, 977–993. [Google Scholar] [CrossRef] [Green Version]
  60. Luu, K.; Bazin, E.; Blum, M.G. pcadapt: An R package to perform genome scans for selection based on principal component analysis. Mol. Ecol. Resour. 2017, 17, 67–77. [Google Scholar] [CrossRef]
  61. Narum, S.R.; Hess, J.E. Comparison of FST outlier tests for SNP loci under selection. Mol. Ecol. Resour. 2011, 11, 184–194. [Google Scholar] [CrossRef] [PubMed]
  62. Villemereuil, P.; Frichot, É.; Bazin, É.; François, O.; Gaggiotti, O.E. Genome scan methods against more complex models: When and how much should we trust them? Mol. Ecol. 2014, 23, 2006–2019. [Google Scholar] [CrossRef] [Green Version]
  63. François, O.; Martins, H.; Caye, K.; Schoville, S.D. Controlling false discoveries in genome scans for selection. Mol. Ecol. 2016, 25, 454–469. [Google Scholar] [CrossRef] [Green Version]
  64. Meirmans, P.G. Genodive version 3.0: Easy-to-use software for the analysis of genetic data of diploids and polyploids. Mol. Ecol. Resour. 2020, 20, 1126–1131. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Pritchard, J.K.; Stephens, M.; Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 2000, 155, 945–959. [Google Scholar] [CrossRef]
  66. Evanno, G.; Regnaut, S.; Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 2005, 14, 2611–2620. [Google Scholar] [CrossRef]
  67. Earl, D.A. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 2012, 4, 359–361. [Google Scholar] [CrossRef]
  68. Noskova, E.; Ulyantsev, V.; Koepfli, K.P.; O’Brien, S.J.; Dobrynin, P. GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data. GigaScience 2020, 9, giaa005. [Google Scholar] [CrossRef]
  69. Jouganous, J.; Long, W.; Ragsdale, A.P.; Gravel, S. Inferring the joint demographic history of multiple populations: Beyond the diffusion approximation. Genetics 2017, 206, 1549–1567. [Google Scholar] [CrossRef] [Green Version]
  70. Schrider, D.R.; Houle, D.; Lynch, M.; Hahn, M.W. Rates and genomic consequences of spontaneous mutational events in Drosophila melanogaster. Genetics 2013, 194, 937–954. [Google Scholar] [CrossRef] [Green Version]
  71. Keith, N.; Tucker, A.E.; Jackson, C.E.; Sung, W.; Lledó, J.I.L.; Schrider, D.R.; Schaack, S.; Dudycha, J.L.; Ackerman, M.; Younge, A.J.; et al. High mutational rates of large-scale duplication and deletion in Daphnia pulex. Genome Res. 2015, 26, 60–69. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Mantelatto, F.L.M.; Fransozo, A. Reproductive biology and moulting cycle of the crab Callinectes ornatus (Decapoda, Portunidae) from the Ubatuba region, São Paulo, Brazil. Crustaceana 1999, 72, 63–76. [Google Scholar] [CrossRef]
  73. Silveira, I.C.A.D.; Schmidt, A.C.K.; Campos, E.J.D.; Godoi, S.S.D.; Ikeda, Y. A corrente do Brasil ao largo da costa leste brasileira. Braz. J. Oceanogr. 2000, 48, 171–183. [Google Scholar] [CrossRef]
  74. Lumpkin, R.; Johnson, G.C. Global ocean surface velocities from drifters: Mean, variance, El Niño–Southern Oscillation response, and seasonal cycle. J. Geophys. Res. Oceans 2013, 118, 2992–3006. [Google Scholar] [CrossRef]
  75. Avise, J.C. Phylogeography: The History and Formation of Species; Harvard University Press: Cambridge, MA, USA, 2000. [Google Scholar]
  76. Allio, R.; Donega, S.; Galtier, N.; Nabholz, B. Large variation in the ratio of mitochondrial to nuclear mutation rate across animals: Implications for genetic diversity and the use of mitochondrial DNA as a molecular marker. Mol. Biol. Evol. 2017, 34, 2762–2772. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Hughes, L.C.; Cardoso, Y.P.; Sommer, J.A.; Cifuentes, R.; Cuello, M.; Somoza, G.M.; González-Castro, M.; Malabarba, L.R.; Cussac, V.; Habit, E.M.; et al. Biogeography, habitat transitions and hybridization in a radiation of South American silverside fishes revealed by mitochondrial and genomic RAD data. Mol. Ecol. 2020, 29, 738–751. [Google Scholar] [CrossRef]
  78. Pertierra, L.R.; Segovia, N.I.; Noll, D.; Martinez, P.A.; Pliscoff, P.; Barbosa, A.; Aragón, P.; Raya Rey, A.; Pistorius, P.; Trathan, P.; et al. Cryptic speciation in gentoo penguins is driven by geographic isolation and regional marine conditions: Unforeseen vulnerabilities to global change. Divers. Distrib. 2020, 26, 958–975. [Google Scholar] [CrossRef]
  79. Rocha, L.A. Patterns of distribution and processes of speciation in Brazilian reef fishes. J. Biogeogr. 2003, 30, 1161–1171. [Google Scholar] [CrossRef] [Green Version]
  80. Allendorf, F.W.; Hohenlohe, P.A.; Luikart, G. Genomics and the future of conservation genetics. Nat. Rev. Genet. 2010, 11, 697–709. [Google Scholar] [CrossRef]
  81. Luiz, O.J.; Madin, J.S.; Robertson, D.R.; Rocha, L.A.; Wirtz, P.; Floeter, S.R. Ecological traits influencing range expansion across large oceanic dispersal barriers: Insights from tropical Atlantic reef fishes. Proc. R. Soc. B Biol. Sci. 2011, 279, 1033–1040. [Google Scholar] [CrossRef]
  82. Tosetto, E.G.; Bertrand, A.; Neumann-Leitão, S.; Nogueira Júnior, M. The Amazon River plume, a barrier to animal dispersal in the Western Tropical Atlantic. Sci. Rep. 2022, 12, 537. [Google Scholar] [CrossRef]
  83. Heinrich, S.; Zonneveld, K.A.F. Influence of the Amazon River development and construction of the Central American Seaway on Middle/Late Miocene oceanic conditions at the Ceara Rise. Palaeogeogr. Palaeoclimatol. Palaeoecol. 2013, 386, 599–606. [Google Scholar] [CrossRef]
  84. Ludt, W.B.; Rocha, L.A. Shifting seas: The impacts of Pleistocene sea-level fluctuations on the evolution of tropical marine taxa. J. Biogeogr. 2014, 42, 25–38. [Google Scholar] [CrossRef] [Green Version]
  85. Milliman, J.D.; Summerhayes, C.P.; Barretto, H.T. Quaternary sedimentation on the Amazon continental margin: A model. GSA Bull. 1975, 86, 610–614. [Google Scholar] [CrossRef]
  86. Maslin, M.; Knutz, P.C.; Ramsay, T. 2006. Millennial-scale sea-level control on avulsion events on the Amazon Fan. Quat. Sci. Rev. 2006, 25, 3338–3345. [Google Scholar] [CrossRef]
  87. Figueiredo, J.J.J.P.; Hoorn, C.; Van der Ven, P.; Soares, E. Late Miocene onset of the Amazon River and the Amazon deep-sea fan: Evidence from the Foz do Amazonas Basin. Geology 2009, 37, 619–622. [Google Scholar] [CrossRef]
  88. Wilson, K.E.; Maslin, M.A.; Burns, S.J. Evidence for a prolonged retroflection of the North Brazil Current during glacial stages. Palaeogeogr. Palaeoclim. Palaeoecol. 2011, 301, 86–96. [Google Scholar] [CrossRef] [Green Version]
  89. Feder, J.L.; Egan, S.P.; Nosil, P. The genomics of speciation-with-gene-flow. Trends Genet. 2012, 28, 342–350. [Google Scholar] [CrossRef]
  90. De Jode, A.; Le Moan, A.; Johannesson, K.; Faria, R.; Stankowski, S.; Westram, A.; Butlin, R.; Rafajlović, M.; Fraïsse, C. Ten years of demographic modelling of divergence and speciation in the sea. Evol. Appl. 2022. [Google Scholar] [CrossRef]
  91. Palumbi, S.R. Marine speciation on a small planet. Trends Ecol. Evol. 1992, 7, 114–118. [Google Scholar] [CrossRef]
  92. Bierne, N.; Bonhomme, F.; David, P. Habitat preference and the marine-speciation paradox. Proc. R. Soc. B Biol. Sci. 2003, 270, 1399–1406. [Google Scholar] [CrossRef] [Green Version]
  93. Nosil, P.; Feder, J.L. Genomic divergence during speciation: Causes and consequences. Philos. Trans. R. Soc. Lond. B Biol. Sci. 2012, 367, 332–342. [Google Scholar] [CrossRef]
  94. Losos, J.B.; Glor, R.E. Phylogenetic comparative methods and the geography of speciation. Trends Ecol. Evol. 2003, 18, 220–227. [Google Scholar] [CrossRef]
  95. Edwards, S.V.; Potter, S.; Schmitt, C.J.; Bragg, J.G.; Moritz, C. Reticulation, divergence, and the phylogeography–phylogenetics continuum. Proc. Natl. Acad. Sci. USA 2016, 113, 8025–8032. [Google Scholar] [CrossRef] [Green Version]
  96. Abbott, R.; Albach, D.; Ansell, S.; Arntzen, J.W.; Baird, S.J.; Bierne, N.; Boughman, J.; Brelsford, A.; Buerkle, C.A.; Buggs, R.; et al. Hybridization and speciation. J. Evol. Biol. 2013, 26, 229–246. [Google Scholar] [CrossRef] [Green Version]
  97. Miglietta, M.P.; Faucci, A.; Santini, F. Speciation in the sea: Overview of the symposium and discussion of future directions. Integr. Comp. Biol. 2011, 51, 449–455. [Google Scholar] [CrossRef] [PubMed]
  98. Martin, S.H.; Dasmahapatra, K.K.; Nadeau, N.J.; Salazar, C.; Walters, J.R.; Simpson, F.; Blaxter, M.; Manica, A.; Mallet, J.; Jiggins, C.D. Genome-wide evidence for speciation with gene flow in Heliconius butterflies. Genome Res. 2013, 23, 1817–1828. [Google Scholar] [CrossRef] [Green Version]
  99. Williams, A.B. The swimming crabs of the genus Callinectes (Decapoda: Portunidae). Fish. Bull. 1974, 72, 685–798. [Google Scholar]
  100. De Queiroz, K. Species concepts and species delimitation. Syst. Biol. 2007, 56, 879–886. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Adult male of Callinectes ornatus Ordway, 1863 (ULLZ 15990). Photo by DL Felder.
Figure 1. Adult male of Callinectes ornatus Ordway, 1863 (ULLZ 15990). Photo by DL Felder.
Genes 13 02263 g001
Figure 2. (A) Map of the western Atlantic showing Callinectes ornatus sampling localities. Circles represent the localities used for ddRAD-seq analyses; triangles represent COI new sequences added to [25]’s dataset. Yellow: north group; blue: south group. Blue lines represent the Amazon River and Orinoco River mouth, resulting in the Amazon–Orinoco Plume. Arrows represent oceanic currents; (B) haplotype network result for C. ornatus. The size of the network circles is proportional to the haplotype frequency. Crossed lines indicate the number of substitutions between haplotypes. Colors represent the region from where the haplotype was sampled. North group: yellow; south group: blue. Individuals from within the AOP region are represented in a separate section in the central haplotype of the south group. Twenty substitutions separate the two networks (C) STRUCTURE plot resulting from the analysis of C. ornatus individuals. Each vertical bar represents an individual. Different colors represent different genetic clusters (K). The voucher numbers of each individual and geographic location are indicated under the plot. (D) Smooth FST across the genome by calculating kernel-smoothed FST considering a sliding window of 900 Kbp. Each circle represents a genome region.
Figure 2. (A) Map of the western Atlantic showing Callinectes ornatus sampling localities. Circles represent the localities used for ddRAD-seq analyses; triangles represent COI new sequences added to [25]’s dataset. Yellow: north group; blue: south group. Blue lines represent the Amazon River and Orinoco River mouth, resulting in the Amazon–Orinoco Plume. Arrows represent oceanic currents; (B) haplotype network result for C. ornatus. The size of the network circles is proportional to the haplotype frequency. Crossed lines indicate the number of substitutions between haplotypes. Colors represent the region from where the haplotype was sampled. North group: yellow; south group: blue. Individuals from within the AOP region are represented in a separate section in the central haplotype of the south group. Twenty substitutions separate the two networks (C) STRUCTURE plot resulting from the analysis of C. ornatus individuals. Each vertical bar represents an individual. Different colors represent different genetic clusters (K). The voucher numbers of each individual and geographic location are indicated under the plot. (D) Smooth FST across the genome by calculating kernel-smoothed FST considering a sliding window of 900 Kbp. Each circle represents a genome region.
Genes 13 02263 g002
Figure 3. GADMA model selected (left) and model validation (right). On the graph depicting the demographic model, arrows represent migration direction, the arrow’s thickness represents migration intensity, and blue figures represent genetic groups over time. Graphs depicting model validation indicate a good fit between the simulated demographic model and our data and indicate residuals near zero.
Figure 3. GADMA model selected (left) and model validation (right). On the graph depicting the demographic model, arrows represent migration direction, the arrow’s thickness represents migration intensity, and blue figures represent genetic groups over time. Graphs depicting model validation indicate a good fit between the simulated demographic model and our data and indicate residuals near zero.
Genes 13 02263 g003
Table 1. Summary statistics of Callinectes ornatus genetic groups from the north and south sides of the Amazon–Orinoco Plume (AOP). Values are based on the analysis of 899 SNPs from across the genome. Geographic location: indicates the number of individuals collected on the north or south side of the plume and included in the analyses after all filtering steps. N: number of individuals after all filtering steps. HE: expected heterozygosity. HO: observed heterozygosity. GIS: inbreeding coefficient.
Table 1. Summary statistics of Callinectes ornatus genetic groups from the north and south sides of the Amazon–Orinoco Plume (AOP). Values are based on the analysis of 899 SNPs from across the genome. Geographic location: indicates the number of individuals collected on the north or south side of the plume and included in the analyses after all filtering steps. N: number of individuals after all filtering steps. HE: expected heterozygosity. HO: observed heterozygosity. GIS: inbreeding coefficient.
Geographic LocationNHEHOGIS
North90.2820.1480.282
South440.2360.1650.236
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Peres, P.A.; Bracken-Grissom, H.; Timm, L.E.; Mantelatto, F.L. Genomic Analyses Implicate the Amazon–Orinoco Plume as the Driver of Cryptic Speciation in a Swimming Crab. Genes 2022, 13, 2263. https://doi.org/10.3390/genes13122263

AMA Style

Peres PA, Bracken-Grissom H, Timm LE, Mantelatto FL. Genomic Analyses Implicate the Amazon–Orinoco Plume as the Driver of Cryptic Speciation in a Swimming Crab. Genes. 2022; 13(12):2263. https://doi.org/10.3390/genes13122263

Chicago/Turabian Style

Peres, Pedro A., Heather Bracken-Grissom, Laura E. Timm, and Fernando L. Mantelatto. 2022. "Genomic Analyses Implicate the Amazon–Orinoco Plume as the Driver of Cryptic Speciation in a Swimming Crab" Genes 13, no. 12: 2263. https://doi.org/10.3390/genes13122263

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop