Comparative Phylogeography of Three Marine Species with Different PLD Modes Reveals Two Genetic Breaks Across the Southern Caribbean Sea

1. Introduction

The ocean is considered a continuous environment and an open system. This implies that marine species exhibit genetic connectivity throughout their distribution due to the currents influencing genetic flow between populations [1]. However, most studies have shown that marine species exhibit a level of population genetic subdivision in response to historical, geological, or ecological-oceanographic factors, although in a segregating manner in taxa with contrasting life histories [1,2,3,4,5,6,7,8,9,10,11]. Combining these factors with biological (e.g., reproductive biology, functional traits, and morphological adaptations) and ecological information allows us to understand marine organisms' phylogeographic and genetic structure [2,5,6,9,12,13,14,15,16]. Among the biological traits, the dispersal potential of marine species is proposed as one of the attributes correlated with species phylogeographic and genetic structure [2,6,17], and the dispersal potential is interpreted as the pelagic larval duration (PLD). Several studies argue that species with a prolonged pelagic larval lifespan exhibit low substructuring or panmixia [11,13,18,19,20,21,22], while those species that exhibit a PLD of a few days or lack pelagic larvae [18,19] tend to present phylogeographic breaks or population genetic substructuring in the presence of contemporary or historical barriers [6,10,13]. These barriers involve hydrological processes such as marine currents and gradients in the physicochemical properties of seawater due to continental river discharges and upwelling zones. Other factors include the absence of specific habitats and variation in shoreline geomorphology, including large distances that limit the dispersal of adult or larval stage organisms [5,6,7,9,11,23,24].

The marine sector of Colombia is in the southern of the Caribbean Sea [25,26] and presents different characteristics on its coasts that allow the evaluation of different biogeographical hypotheses [17,26,27,28]. This marine sector has been characterized by various historical events, such as changes in the geomorphology of the coastline due to variations in sea level during the last glaciation and the uplift and movement of Caribbean basin mountain systems, such as the Sierra Nevada de Santa Marta [27]. This movement changed the continental shelf, the direction of marine currents, and the location of the mouth of the Magdalena River [29], which may have affected the phylogeographic structure of the marine-coastal organisms. Indeed, various phylogeographic studies in the Caribbean region used different genetic markers for some species with different larval dispersal times and suggested the existence of a break in the genetic connectivity of populations on both sides of the Caribbean, principally among the central Caribbean and West Indies locations [17,30,31,32,33,34,35,36]. This genetic break is placed between Venezuela and Colombia and may be caused by the Magdalena River plume (MRP). However, several studies have rejected this hypothesis because a sampling design needs to be proposed. For example, some studies have not considered locations in Colombia [32,33,35,36,37] and those that do have only sampled close to one side of the river [34,38]. In other cases, the studies sampled both sides but only considered species with pelagic larvae of more than 12 days in duration and concluded that the Magdalena River plume is a permeable barrier to gene flow [i.e., coral, shrimp, mollusk, sea urchin, and fish species: [21,34,39,40,41,42,43,44,45,46,47]. Therefore, future investigations are necessary to evaluate the MRP effect on other marine species.

A second barrier has been believed to be located in 74 – 71°W, including Santa Marta and the La Guajira Peninsula [17,27,28], which is attributed to the permanent upwelling that transports upwelled water toward the center of the Caribbean Sea and that would limit larval dispersal toward the southwestern Caribbean coast of Colombia [48]. This putative barrier possibly affects the genetic and phylogeographic structure of species associated with coral reefs and shallow rocky bottoms, whose marine ecosystems are absent for more than 300 km of coastline between Cabo de la Vela (La Guajira) and Tayrona National Natural Park, where upwelling occurs (Figure 1).

Mollusks and fish are the most diverse and abundant species groups in the Colombian Caribbean [28,49], and serve as the biological models for investigating how the Magdalena River plume (MRP) and the combination of the absence of shallow rocky bottoms (ARB) with permanent upwelling in La Guajira (PUG) affect the phylogeography of their populations. For this purpose, we selected three marine species with different PLDs: the reef fish Acanthemblemaria rivasi (with PLD < 22 days) [50,51], the exposed rocky shore snails Cittarium pica (with PLD < 6 days) [52,53] and Nerita tessellata (with PLD > 60 days) [54,55,56]. With this design, we followed aspect III of the multispecies approach proposed by [4], which suggests comparative phylogeography studies to investigate whether the same factors that cause spatial patterns of lineages can affect multiple codistributed taxa [4,57]. Therefore, we evaluated the following questions: Are the Magdalena River plume (MRP) and the combination of the absence of a rocky bottom and permanent upwelling in La Guajira (ARS+PUG) physical barriers influencing the phylogeographic patterns of marine species with varying dispersal potential across the southern Caribbean Sea? What is the population model of each species? Furthermore, what is the level of congruence between the phylogeographic patterns of the three species?

2. Materials and Methods

2.1. Sampling

For each species, 95 specimens were collected across five marine sectors (1. Cabo de la Vela (12° 27' 29.3" N, 71° 39' 59.6" W - 12° 12' 58.0" N, 72° 10' 41.1" W); 2. Santa Marta-Tayrona National Natural Park (11° 18' 51.6" N, 74° 11' 39.1" W - 11° 20' 17.7" N, 74° 03' 14.2" W); 3. Cartagena-PNNN Corales del Rosario y San Bernardo (10° 22' 29.3" N, 75° 34' 16.9" W - 10° 15' 22.2" N, 75° 37' 04.5" W - 10° 10' 31.2" N, 75° 46' 07.1" W); 4. Isla Fuerte (9° 23' 15.4" N, 76° 10' 23.6" W), and 5. Capurganá-Sapzurro (8° 38' 30.6" N, 77° 20' 45.0" W - 8° 38' 40.4" N, 77° 20' 14.6" W - 8° 39' 40.6" N, 77° 21' 40.3" W; Figure 1a; Table S1). The sampling was performed by scuba diving (Acanthemblemaria rivasi) and hand collection (Cittarium pica and Nerita tessellata). Foot (mollusks) and fin tissue samples were fixed with 96% ethanol. The final number of processed samples per location for each species is presented in Table S1.

2.2. Laboratory Procedures

To obtain the SNP loci, double digest restriction-site associated DNA (ddRADseq) was developed using the library preparation protocol of [58]. The library preparation and sequencing of ddRADseq for 95 individuals were performed using the enzymes EcoRI and HpyCh4IV for A. rivasi, EcoRI and NlaIII for C. pica, and EcoRI and HpyCh4IV for N. tessellata. The laboratory procedures were developed by the Australian Genome Research Facility. The libraries were size selected, targeting fragments 280–375 bp in size with Blue Pippin (Sage Science), and PCR amplified with indexed primers. Finally, the single-ended fragments were sequenced on the Illumina NextSeq500 with 150 cycles in high output mode [59].

2.3. De Novo Assembly Reads

The processes to demultiplex reads by barcode and to remove reads with poor sequence quality or uncalled bases were developed in the program process_radtags in STACKS v2.0 [60,61]. Prior to the assembly of the RAD sequences, the optimal values of parameters m (minimum number of raw reads required to form a stack of a putative allele), M (number of mismatches allowed between putative alleles to merge them into a putative locus), and n (number of mismatches allowed between putative loci during the construction of the catalog) were obtained following the recommendations of [62]. First, the RADProc pipeline [63] was executed by selecting the 20 individuals (4 per location) with the highest number of reads. This program uses the same parameters defined in STACKS (m, M, and n) for de novo locus formation and catalog building. The output of the results allowed us to graph the number of loci, the number of polymorphic loci, and the number of SNPs of each combination of the parameter values [62], with the highest numbers observed when M = 6, m = 3, and n = 2 for A. rivasi, M = 5, m = 3, and n = 2 for C. pica, and M = 4, m = 3, and n = 2 for N. tessellata. These values were established to assemble the rad loci with the denovo_map program in STACKS v.2.0 [61].

Without a reference genome, RAD-seq loci were assembled de novo using the denovo_map program in STACKS [61]. The final filtered genotype file was created by applying the following filters: a minimum percentage of individuals in the populations required to process a locus (-R = 0.8) equivalent to 80% and a minor allele frequency (MAF) (--min_maf = 0.02). The negative effects of missing data were minimized when calculating frequency-based genetic distance parameters. SNPs were also screened for allele coverage, with any SNPs displaying a local and global minor allele frequency (MAF) threshold of less than 1% removed from the dataset. In cases where multiple SNPs were found within the same read, only one locus was retained (--write_single_snp) in order to avoid statistical bias from the physical linkage. Individuals who exceeded the missing data percentage (less than 5%) were removed from the final file. After the final filtering, 85, 65 and 51 individuals of A. rivasi, C. pica and N. tessellata, respectively, remained. The VCF file was also converted into other formats using PGDSpider v2.1 software [64].

2.4. Analysis of Data

2.4.1. Genetic and Phylogeographic Structure

Genetic diversity indices were estimated by the number of alleles, effective number of alleles, observed heterozygosity (Ho), heterozygosity within populations (Hs) and inbreeding coefficient (G_IS), which were generated with GenoDive [65]. As the study aims to determine levels of variation at the intrallocality level upon biogeographic barriers, the results of these parameters are presented in Table S1.

We determined the spatial genetic structure of each species to identify the probable number of populations (K) using the Bayesian approach with the program STRUCTURE v.2.3.3 [66]. Each run included 200,000 burn-in iterations, followed by 800,000 MCMC iterations used to calculate a probable K between one and n + 1 populations in each case, with seven runs for each simulated K value. Admixture and allele frequency correlated models were assumed here. We determined the number of most likely populations with the methods proposed by [67,68,69] using the STRUCTURE SELECTOR program [70], which includes the CLUMPAK program [71], to combine and visualize the STRUCTURE results. Individuals within populations were ordered according to sampling localities. A principal component analysis (PCA) was also performed for each species with the dudi.pca function of the ade4 package [72] in R.

We performed estimates of Ф_ST to determine genetic differentiation between pairs of localities with STACKS [61] using the "populations" module with the "--fstats" indicator. In addition, we performed an analysis of molecular variance (AMOVA) to determine the effects of putative barriers on the phylogeographic pattern of each species. This evaluated several levels of grouping according to the localities on each side of the putative barriers: 1. the effect of the absence of shallow rocky bottoms for more than 300 kilometers of coastline between Cabo de La Vela and Santa Marta, as well as the permanent upwelling of La Guajira (ARB+PUG); 2. the effect of the Magdalena River plume; 3. a third analysis was performed for the species that presented substructuring according to the results of the Bayesian analysis of STRUCTURE (only for C. pica and A. rivasi). AMOVAs were performed with the poppr.amova function of the poppr package in R [73].

The phylogeographic pattern analysis was performed from phylogenetic analyses based on consensus sequences for all concatenated individuals. The files were generated using the STACKS population program with the --phylip_var function based on 477,072 SNPs for A. rivasi, 367,676 SNPs for C. pica, and 229,422 SNPs for N. tessellata. The generated matrices were used to perform maximum likelihood (ML) analyses with the program IQTREE 1.3.10 [74], with two independent runs, including model selection with SNP verification bias correction, 1000 ultrafast bootstrap replicates, and 10000 iterations. The SNP-based sequence evolution model was selected by IQTREE following the Bayesian Information Criterion, with TVM+ASC+G4 for A. rivasi, TVMe+ASC+G4 for C. pica, and TVMe+ASC+G4 for N. tessellata.

2.4.2. Population Model

We determined the population model of each species (open, isolated, or hierarchical populations) with the Mantel test to correlate the genetic (Ф_ST) and geographic (km) distance matrices using the mantel function of the R package vegan 2.5.7 [75]. Two geographic distances were calculated from Google Earth: one directly measured the distance between locations (linear), and the other measured along the coastline. We developed the analyses without any previous manipulation of data [76]. Additionally, a Mantel correlogram analysis was also performed, in which geographic distances were divided into class-marked submatrices [76]. We did this to determine false-positives of isolation by distance when the populations were significantly understructured because of biogeographic barriers [77].

2.4.3. Phylogeographic Concordance between Species

The phylogeographic concordance analysis among species was performed by comparing the topology of the dendrograms constructed with the pairwise Ф_ST matrices calculated between localities. In this sense, the Kendall procedure [78] was used to compare tree topology according to tip category names using a Concordance Factor (CF). The CF measure can estimate values between 0 and 1. It is 1 when the concordance is complete. For this measure, comparisons were made between A. rivasi and C. pica, A. rivasi and N. tessellata, and C. pica and N. tessellata. The procedure was performed with the treeConcordance function of the R package treeSpace [79], entering the matrix in the Newick format of each dendrogram. Another estimator employed was the congruence between distance matrices test (CADM) [80], which was used to determine topological and genetic congruence. This test evaluates the hypothesis of complete incongruence between species trees, corresponding to phylogenies with different topologies or branch lengths [80]. Kendall's W statistic was used to determine the level of congruence between species, which has values between 0 and 1, with 0 indicating no congruence and 1 being complete topological or genetic congruence of the trees. The topological and genetic congruence tests used the Ф_ST matrices calculated between pairs of locations.

The CADM test was run for each congruence type and between pairs of species, with the CADM.global function of the R package Ape [81]. Only for the topological congruence test was the cophenetic.phylo function of the R package Ape [81] previously executed. The significance value was estimated from 99,999 permutations. As a complement, the Mantel test was performed between the triangular matrices of the paired Ф_ST values of the codistributed species with the mantel function of the R vegan 2.5.7 package [75]. Congruence between the matrices is indicated by r values close to 1 and p <0.05. The significance value was estimated from 99,999 permutations.

3. Results

We obtained 10,122,277 reads for 85 individuals of Acanthemblemaria rivasi, 8,764,529 reads for 65 individuals of Cittarium pica, and 3,385,649 reads for 51 individuals of Nerita tessellata. The average number of usable polymorphic SNP loci per sample varied: 66,325 were determined for A. rivasi, 55,151 for C. pica, and 21,333 for N. tessellata. A summary of genetic diversity indices is presented in Table S1.

3.1. Genetic and Phylogeographic Structure

Acanthemblemaria rivasi exhibited considerable genetic differentiation between samples from Cabo de la Vela and the other localities (Ф_ST > 0.290 values; p < 0.05), with greater magnitude observed for the most distant location (Capurganá; Ф_ST = 0.345; p < 0.05). Samples from Santa Marta also showed high genetic differentiation compared to the rest of the localities (Ф_ST > 0.262; p < 0.05), while those collected in Cartagena, Isla Fuerte, and Capurganá had moderate differentiation among them (Ф_ST < 0.097; Table 2). The results of the Bayesian analysis provided variable results among the estimates of the probable number of populations. For example, the ∆K estimate provided a K = 2. According to the bar graph, the first population is distributed in Cabo de la Vela and the second in Santa Marta, Cartagena, Isla Fuerte, and Capurganá (Figure 1b). MedMedK and MedMeanK indicated a K = 3, with the first population located at Cabo de la Vela again, the second in Santa Marta, and the third across Cartagena, Isla Fuerte, and Capurganá. For the LnP(K), MaxMedK, and MaxMeanK estimates, K was 5. These estimates proposed the same first three populations as with K=3, but with the fourth (higher proportion) and fifth populations (in a minimum proportion) codistributed mainly between Cartagena and Isla Fuerte (Figure S1a). Principal component analysis (PCA) and the network tree showed three clusters, which are congruent with the K = 3 of the Bayesian analysis and the levels of differentiation based on Ф_ST (Figure 2a).

In the case of C. pica, the samples from Cabo de la Vela were also differentiated from the rest of the localities. Although this genetic differentiation was moderate (Ф_ST > 0.160; p < 0.05), the highest value was observed for the neighboring locality of Santa Marta (Ф_ST = 0.160; p < 0.05). The localities of Cartagena, Isla Fuerte, and Capurganá presented a lower level of differentiation between their samples (Ф_ST < 0.013; p < 0.05; Table 2). Concerning Bayesian analysis, the LnP(K), MedMedK, MedMeanK, MaxMedK, and MaxMeanK estimates calculated a K = 3, with one population located in Cabo de la Vela, the second in Santa Marta, with a minimal proportion in Cartagena (< 5%), and the third codistributed across Cartagena, Isla Fuerte, and Capurganá (Figure 1c). The ∆K estimated K = 5, the same as K = 3, but with the fourth and fifth in minimal proportions (< 3%) and codistributed between Isla Fuerte and Capurganá (Figure S1b). The PCA is congruent with the proposed K = 3 and the Ф_ST results (Figure 2b).

Nerita tessellata presented slight genetic differentiation among all localities, with an average value of Ф_ST = 0.006, which was the highest, but of subtle differentiation, in the comparison of those of Cabo de la Vela with the rest of the localities (Ф_ST < 0.009; Table 2). The MedMedK, MedMeanK, MaxMedK, and MaxMeanK estimates calculated a codistributed K = 1 at all localities (Figure 1d). For the cases of ∆K and LnP(K), K values equal to 3 and 2, respectively, were proposed. The two clusters are distributed from Cabo de la Vela to Capurganá, one with a higher proportion than the other. The third cluster proposed by LnP(K) is distributed in a minimal proportion in all localities except Cabo de la Vela (Figure S1c). However, the PCA and network tree agree with the K = 1 proposal (Figure 2c).

3.2. Identification of Phylogeographic Breaks

The genetic differentiation analysis by AMOVA presented a level of phylogeographic structuring in all species (Ф_ST > 0.076, p < 0.05), except in N. tessellata (Ф_ST = 0.001 - 0.002, p > 0.05), which did not present substructuring (Table 3). For A. rivasi, the leading causes of the phylogeographic pattern were two putative barriers: 1. the combination of the absence of shallow rocky bottoms and permanent upwelling (ARB+PUG) located between Cabo de la Vela and Santa Marta (AMOVA: Ф_CT = 0.406, p < 0.05) and 2. the Magdalena River plume (MRP; Ф_CT = 0.420, p < 0.05). These two barriers delimit the three populations that were identified for A. rivasi. Under this last scenario, a third AMOVA was performed assuming K = 3, determining a high genetic differentiation (Ф_CT = 0.495, p < 0.05; Table 3). The maximum likelihood (ML) tree for A. rivasi was concordant with the AMOVA when K = 3 was assumed, confirming the effects of the two putative barriers (Figure 3a). The three clades are observed with 100% bootstrap support. The same as in the Bayesian analysis when K = 3, one clade is constituted by the Cabo de la Vela samples, the second by those from Santa Marta (100% bootstrap), and the third by Cartagena + Isla Fuerte + Capurganá samples (100%).

Concerning C. pica, when the ARB+PUG barrier was assumed, substantial genetic differentiation was observed between samples from both sides (Ф_CT = 0.224, p < 0.05) compared to the MRP barrier, which showed moderate genetic differentiation (Ф_CT = 0.076, p < 0.05). Furthermore, when the analysis was performed assuming K = 3, the differentiation was significant (Ф_CT = 0.174, p < 0.05; Table 3). However, the phylogenetic analysis was concordant with the ARB+PUG barrier (100% bootstrap). Only two clades are configured, one on each side of the barrier (clade 1. Cabo de la Vela and clade 2. Santa Marta + Cartagena + Isla Fuerte + Capurganá). Interestingly, clade 2 shows how the Santa Marta samples are more closely related but separated from the rest of the localities, forming a subclade with 53% bootstrap support (Figure 3b). With respect to N. tessellata, no evidence of phylogeographic breaks was observed (Figure 3c).

3.3. Population Model

Only A. rivasi showed a significant correlation between the genetic distance and linear and coastline distances (r_{m coastline} = 0.623, p = 0.0253; r_{m linear} = 0.774, p = 0.0118; Figure 4a), while the analysis of C. pica (r_{m coastline} = 0.441, p = 0.2059; r_{m linear} = 0.616, p = 0.0571; Figure 4b) and N. tessellata (r_{m coastline} = 0.422, p = 0.227; r_{m linear} = 0.580, p = 0.0831; Figure 4c) was not significant. However, the Mantel correlogram indicated that A. rivasi and C. pica fit the hierarchical population model (Figure S2).

3.4. Phylogeographic Concordance between Species

A concordance factor of 0.6 was determined when comparing the topologies of the dendrograms of A. rivasi-N. tessellata and C. pica-N. tessellata and of 0.5 for those of A. rivasi-C. pica. The global CADM test showed significant topological (Wt = 0.847, X² = 22.86, p = 0.012) and genetic (Wg = 0.852, X² = 23, p = 0.011) differences among the three species. However, in pairwise comparisons, only A. rivasi-C. pica showed significant topological (Wt = 0.937, X² = 16.88, p = 0.046) and genetic (Wg = 0.915, X² = 16.42, p = 0.048) correlations, as well as a significant correlation between genetic distance matrices (r_m = 0.7615, p = 0.001). N. tessellata only presented genetic congruence with A. rivasi (Wg = 0.915, X² = 16.47, p = 0.043; Figure S3a-c).

4. Discussion

4.1. Genetic and Phylogeographic Structure

Genetic structuring analyses based on paired Ф_ST, Bayesian analysis, PCA, AMOVA, and ML trees determined the most likely number of clusters for the three marine species sampled along the southern Caribbean Sea (Colombia sector). However, only N. tessellata presented one population (K = 1), and its high dispersal potential reflected in a larval lifespan of more than 60 days [54,55,56] may be the primary explanation for the absence of genetic structuring. The overall Ф_ST value was low and coincides with those determined for other marine species associated with the southern Caribbean ecoregion, which were estimated with various molecular markers (mitochondrial genes and SNPs, and microsatellite loci), mainly in the fishes Stegastes partitus [39], Lutjanus synagris [40], Caranx hippos [42] and Micropogonias furnieri [82], the rock boring urchin Echinometra lucunter lucunter [44], the southern white shrimp Litopenaeus schmitti [83], and the Caribbean sharpnose shark Rhizoprionodon porosus [84]. Although some Ф_ST values were significant, the levels of genetic differentiation between localities were too low to prove any phylogeographic break (F_ST < 0.05), and the high dispersal potential of these species may be a reason for the lack of structure (PLD > 10 days). Other studies that extended the sampling coverage to other regions of the Great Caribbean documented that species with high PLD do not show a phylogeographic pattern [i.e., Echinolittorina ziczac and Cenchritis muricatus, [41]; Panulirus argus, [85]], and some species were described to have a panmictic population pattern [Sparisoma viride,[21]; Clibanarius tricolor, [22]; Sparisoma aurofrenatum, [46]]. Overall, the timing of the larval period of N. tessellata defined the phylogeographic structure with a single population. It also exhibits high genetic connectivity between localities due to the functioning of the Caribbean current and the Panama-Colombia countercurrent, which could generate larval dispersal along the southern Caribbean in the presence of the barriers evaluated (Figure 1a).

Concerning A. rivasi and C. pica, pairwise Ф_ST, Bayesian analysis, and PCA allowed the identification of three populations (K = 3), which were confirmed with AMOVA. The delimitation of the three populations is influenced by the putative barriers evaluated. For example, based on the Ф_CT statistic of the AMOVA, it was possible to measure the level of genetic differentiation between the samples classified on either side of the barriers. For both species, the first population was in Cabo de la Vela and the second was in Santa Marta. The populations are separated by the barrier defined as a combined effect of the absence of the rocky bottom at a distance greater than 300 km from the coastline (ARB) and the permanent upwelling in La Guajira (hereafter ARB+PUG). This barrier significantly affects A. rivasi (Ф_CT = 0.406) and C. pica (Ф_CT = 0.224), but in a minimal proportion in N. tessellata (Ф_CT = 0.005). Finally, the analysis determined a third population between the localities of Cartagena, Isla Fuerte, and Capurganá, separated from the second by the barrier of the Magdalena River plume (MRP). This barrier significantly affects A. rivasi (Ф_CT = 0.420) and, to a lesser degree, C. pica (Ф_CT = 0.076): for the latter, it is a barrier permeable to gene flow. Cittarium pica presents a PLD of less than six days [52,53]. This time may be sufficient for larvae produced in Cartagena to be transported to Santa Marta via the Caribbean Counter Current (CCC). This current, which originates from the Panama-Colombia Gyre [48,86], transports surface waters over the continental shelf to the north and northeast to the La Guajira front when the northeasterly trade winds weaken [48,86,87]. The CCC may also be responsible for maintaining high genetic connectivity between Capurganá and Cartagena, where the third populations of C. pica and A. rivasi are found. It is assumed that the larvae of both species produced in the reef systems of the Cartagena (Tierrabomba, Barú Island, Rosario Islands, and San Bernardo archipelago), Isla Fuerte, and Capurganá sectors are transported by CCC [87]. This current has been proposed as the oceanographic factor that facilitates the exchange of genetic information between populations of several marine species that are not significantly understructured in the southwestern Caribbean [i.e., S. partitus, [39]; L. synagris, [40]; L. schmitti, [83]; R. porosus, [88]; C. hippos, [42]; Orbicella faveolata, [38]; P. notialis, [45]; Mugil liza, [89]; Acropora palmata and A. cervicornis, [47]; Micropogonias furnieri, [82]].

The phylogenetic analysis confirmed the previous phylogeographic proposal for A. rivasi and C. pica. For C. pica, only two clades were configured in the ML tree, one on each side of the ARB + PUG barrier. Nevertheless, clade 2 shows how the samples from Santa Marta are more closely related but separated from the rest of the localities and formed a subclade within it (53% support). This ML tree suggests that C. pica shows evidence of the effect of a single biogeographic barrier (ARB+PUG). In contrast, for A. rivasi, the simultaneous impact of the two barriers evaluated (ARB+PUG and MRP) is observed. These results are interesting because the action of two barriers was identified for a reef species such as A. rivasi less than 400 km from the coastline. Interestingly, species of the genus Acanthemblemaria have a proposed PLD of 22 to 25 days [51]. If this is true for A. rivasi, under this dispersal scenario, it would be affected by only one of the two barriers, as in the case of C. pica. However, this fish could have a PLD of fewer than ten days, although this must be investigated. Another reason for explaining the phylogeographic pattern of A. rivasi is that Acanthemblemaria fishes are characterized by parental care of eggs by males, low fecundity, and inhabiting invertebrate skeletal orifices embedded in rocks or corals in reef areas with high wave energy [50,51,90]. In addition, these species show larval retention near coral reefs [91,92]. The reproductive and ecological characteristics mentioned above could generate high larval retention and high biological recruitment in the local populations of Cabo de la Vela and Santa Marta, which, added to the effects of ARB+PUG and MRP barriers, define the high level of genetic structuring observed.

The combination of physical, ecological, and biological aspects suggests a vicariant event for A. rivasi because a new species distributed from Santa Marta to Venezuela was recently proposed [A. aceroi, [93]]. However, the high level of genetic differentiation observed among the three clades (Cabo de la Vela, Santa Marta, and Cartagena+Isla Fuerte+Capurganá; AMOVA: Ф_CT = 0.495, p < 0.05) suggests a taxonomic revision, which should include the mitochondrial molecular markers, meristic characters, and morphology to determine the existence of a possible additional species in the A. rivasi complex.

4.2. Identification of Phylogeographic Breaks

The phylogeographic analysis of the three marine species determined the action of two putative barriers across the southern Caribbean (Colombia sector). For the first time, evidence is provided of a phylogeographic break caused by the Magdalena River Plume (MRP), mainly for A. rivasi. The Magdalena River annually delivers 142 x 106 tons year^-1 of sediment to the Caribbean Sea. The MRP extends to 6.5 km offshore and is characterized by a high turbidity concentration (178.6 ± 78.7 mg L^-1), comparable to estimates for the Amazon and Yangtze rivers [94]. Furthermore, MRP manages to dilute salinity from 10.8 ± 3.4 at the mouth to 28.4 ± 0.4 in the adjacent marine sector [95]. These physicochemical conditions and their influence on the Caribbean Sea suggest that the PMR must act as a biogeographic barrier for marine species that are limited in their dispersal, mainly those whose pelagic larvae are unable to survive when attempting to cross this area. Perhaps they cannot tolerate the decrease in salinity, higher temperature (≥ 2 °C) than the marine waters off Santa Marta-Cabo de la Vela, and the high turbidity of the waters, which are possible factors that regulate the success of biological dispersal. The MRP likely operates as a filter for A. rivasi larvae that attempt to cross it through the action of the CCC. However, the possibility that populations formed on either side of the plume have adapted to the different environmental conditions is not excluded. For example, those larvae that cross the MRP may be selected against because the marine waters on the western side are warmer and less salty than those on the eastern side [96].

The MRP was proposed as a barrier across the southern Caribbean [34]. However, many studies were inconclusive because some of them did not consider localities in Colombia [32,33,35,36,37], and when it was possible, they only sampled one side of the MRP [34,38]. The hypothesis of the MRP as a barrier was rejected when studies included both sides of MRP and selected species with pelagic larvae enduring longer than 12 days (S. partitus, L. synagris, L. schmitti, C. hippos, M. incilis, P. notialis, M. liza, M. furnieri, A. cervicornis, and A. palmata). Nevertheless, other world scenarios describe marine species with phylogeographic breaks caused by large river plumes. The Amazon River plume, for example, has been the subject of biogeographic studies on reef fishes, demonstrating the separation of the fish fauna of the Caribbean Sea from Brazil [12]. Recent reviews also document it as a barrier to dispersing species from several taxonomic groups, generating a spatial pattern in marine species diversity between the Greater Caribbean and Brazil [97]. From a phylogeographic point of view, some species are not affected by the Amazon River plume; for example, the reef fishes Chaetodon striatus [98] and Abudefduf saxatilis [99], showed genetic flow between localities in the Caribbean Sea and northern Brazil. The current from northern Brazil flowing toward the Caribbean, passing through the Amazon River plume platform, promotes the genetic connectivity of A. saxatilis between both sectors [99]. In contrast, the Amazon River plum has been shown to act as a barrier for the reef fishes Chromis multilineata and those of the genus Halichoeres, demonstrating levels of genetic differentiation between samples from the Greater Caribbean and localities from Brazil and the Central Atlantic islands [100,101]. This was also observed in the striped snapper L. synagris when comparing samples from Colombia with Brazil [24]. However, another phylogeographic break caused by a river plume has been documented. The limpet Cellana toreuma is distributed along the coasts of China, where the Yangtze River plume causes a phylogeographic break for this species [102].

The other finding was the effect caused by the absence of the rocky coastline for more than 300 km between Cabo de la Vela and the TNNP sector, which operates for A. rivasi and C. pica. The rocky coastline constitutes the specific habitat for both species, which heterogeneous distributions in the southern Caribbean [44,103]. However, the absence of the rocky coastline in this 300 km sector coincides with two upwelling areas. One is almost permanent in La Guajira, with strong upwelling during December-May and a weak effect during June-August when the temperature reaches a minimum of 24 °C (~72°W, Figure 1A; see review in [104]). The other upwelling develops between Santa Marta and TNNP and is usually seasonal (December to March), with a minimum of 24 °C [48,86,104]. In both cases, surface currents transport upwelled water offshore [48,86,104], which is then added to water transported by the Caribbean Current toward the Central Caribbean [48,104] (Figure 1a). This oceanographic feature may be responsible for transporting larvae offshore [104] and regulating the genetic exchange between the Cabo de la Vela and Santa Marta localities. Therefore, it could be a combined effect which generate the phylogeographic break in the two species. This barrier had not been documented with scientific support [44,46] and is interesting because the two species evaluated contrast in their larval life histories, a sign of the generalized effect of the barrier on the communities of marine organisms dependent on the rocky bottoms across the southern Caribbean [Aspect III, [4]]. A similar situation occurs in Norway with the rocky shore fish Symphodus melops [105] and southeastern Australia with the barnacle Catomerus polymerus, and the limpet Cellana tramoserica, [5]], where these marine species associated with the rocky shore exhibits a break ingenetic connectivity due to a biogeographic barrier formed by sandy coastlines.

Some investigations have explained how upwelling zones affect the genetic and phylogeographic structures of marine species. For example, upwelling events at Cape Blanco (Oregon) and Cape Mendocino (California) affect the genetic structure of the barnacle Balanus glandula [106] and five species in the rocky intertidal community [19], respectively. Along the southeastern Pacific coast, the gastropod Crepipatella dilatata [107] and the beach isopod Excirolana hirsuticauda [108] exhibit a break at 32 S, a transition area characterized by upwelling. In addition, the fish Pomatomus saltatrix exhibits this phenomenon in the upwelling area of the Benguela Current [109], and the fish Sebastes thompsoni ocurrs between two sectors of the East Sea, which may be related to current patterns such as eddies and upwelling [110].

The above demonstrates the importance of further investigating the effects of the MRP and ARB+PUG barriers on other species inhabiting the rocky shores and reefs of the Southern Caribbean. Ideally, future research should use a multispecies approach that includes species from different taxonomic groups to test whether the effects of the barrier are widespread throughout the marine community associated with this ecosystem [4,111]. In addition, research should also be based on the multilocus approach [4] to evaluate whether barriers have generalized genomic effects on marine species.

4.3. Population Model

One of the assumptions in population genetics is that genetically structured species fit the isolation by distance (IBD) model, in which samples from nearby locations are less genetically different than those farther apart [112]. However, highly structured species may exhibit false-positives when testing the IBD model for spatial autocorrelation between genetic and geographic distances [76,77]. Some species experienced abrupt separations in the past that are evidenced by high levels of genetic divergence in the present. Thus, Mantel test analyses and Mantel correlograms performed on A. rivasi and C. pica determined that none of them fit the IBD model. In contrast, these analyses demonstrated that they fit the hierarchical population model [77]. This model coincides with those species that show an abrupt genetic change in the face of a biogeographic barrier in a geographical area [1,77,113]. For example, the regression plots could show a linear relationship between geographical and genetic distances, although in some cases not to a significant degree, as was observed in C. pica.

All of this may occur due to the genetic relationships between subpopulations within each clade that bias the test by suggesting a fit to the IBD model. For A. rivasi and C. pica, it was possible to detect the effects of barriers on the relationship between the two distances (Figure 4). In the correlation plots of these species, three relationship groups are observed. One is found toward the upper part of the graph, where the pairs of data from Cabo de la Vela are related to the rest of the localities. Another group is toward the lower part of the graph, where the relationships of the localities of Cartagena, Isla Fuerte, and Capurganá are observed, and another is in the center, including the relationships of the Santa Marta data with those of Cartagena, Isla Fuerte, and Capurganá. For N. tessellata, the Mantel test analyses showed a low correlation between geographic and genetic distances, which was not significant. The slope values of the linear regressions between the two variables tended to be zero. These results, those of the AMOVA, and the low value of Ф_ST suggest that N. tessellata fits the model of open or panmictic populations across the southern Caribbean [1].

Species that do not fit the IBD population model have been investigated across the geographical area studied, principally those species with high PLD values and low levels of genetic differentiation that were estimated from microsatellite loci, as in the cases of S. partitus [39], L. synagris [40], M. incilis [43], E. lucunter lucunter [44], Mugil liza [89], and C. hippos with mitochondrial genes [42]. No effect of geographic distances on genetic distances suggests a process of genetic homogenization among the studied localities. However, the literature review identified that the shrimp P. notialis is the only species that fits the IBD model across the southern Caribbean [45]. For this species, a K of three populations was determined, with the samples from La Guajira being the most genetically different from the rest of the localities in Colombia.

4.4. Phylogeographic Concordance between Species

Phylogeographic concordance factors allow us to quantify the degree of phylogeographic congruence among species to determine which species share a phylogeographic pattern. In the past, phylogeographic studies compared phylogenies among multiple taxa qualitatively and interpreted similarity as a shared response to a historical event [4,114,115]. However, most work concludes that not all codistributed species exhibit identical phylogenies due to incongruence in tree topology and clade divergence times. This argument is posited because each species responds differentially to the factors responsible for the phylogeographic break and to the biological characteristics of the species involved [5,57,116].

In this study, only two species presented phylogeographic congruence. Acanthemblemaria rivasi and C. pica coincided in presenting phylogeographic breaks to the putative barriers evaluated but with a differential response. For example, A. rivasi showed two phylogeographic breaks caused by the ARB+UPG and MRP barriers, while C. pica was only affected by ARB+UPG. This explains why the Kendall concordance factor [79] presented a value of 0.5, indicating some degree of incongruence with respect to the distribution of the topology categories of their dendrograms. Nevertheless, these species showed a high level of genetic and topological congruence as determined by the method of [80]. The biological and genetic attributes of both species are likely responsible for the level of congruence observed: low dispersal potential attributed to a PLD of few (6 days in C. pica) to several days (22 days in A. rivasi); presenting egg parental care (A. rivasi) and low fecundity (A. rivasi); high Ф_ST values (> 0.08); a K = 3 (A. rivasi, C. pica); and two (C. pica) to three clades (A. rivasi) delineated by the phylogenetic analyses. N. tessellata is characterized as having high dispersal potential with a PLD greater than 60 days, high fecundity, no parental care, and is genetically unstructured (Ф_ST < 0.01; K = 1; 1 clade).

4.4. Conservation Aspects

For A. rivasi and C. pica, three populations were identified: one in Cabo de la Vela, another in Santa Marta, and the third codistributed across Cartagena, Isla Fuerte, and Capurganá. This spatial pattern coincides with the zoogeographic subareas proposed by [28] as a result of a study of the biogeography of marine gastropods in the southern Caribbean. This author proposed five subareas, three of which are distributed in the marine sector of Colombia. The first is located from the mouth of the Magdalena River to Costa Rica (Isthmian), which contains the Colombian localities of Cartagena, Isla Fuerte, and Capurganá. The second includes the marine sector from the eastern side of the mouth of the Magdalena River to the eastern side of Tayrona Natural National Park (Samarian), where the third subarea begins and reaches the Paraguaná Peninsula in Venezuela (Goajira). This subarea includes the Cabo de la Vela location. Overall, these three subareas are characterized by differences in the size of the continental shelf, in the predominant types of bottoms, in the conditions of quietude, transparency, temperature, and salinity of the water masses, and in the different types of habitats that they contain [see details in [28]].

The physical, chemical and geological conditions of the three subareas exert substantial effects on the coastal marine species of the Caribbean of Colombia, which allowed us to observe a typical pattern of genetic differentiation between the samples collected in the Goajira subarea and the other subareas [i.e., C. mapale, [27]; E. lucunter, [44]; P. notialis, [45]; R. porosus, [88]; Sciades proops, and Melongena melongena [82]]. Thus, the marine sector of La Guajira in Colombia should be considered a particular management area for the conservation of marine species, mainly for those that are exploited by fishing activities, as in the cases of C. pica, C. mapale, S. proops, P. notialis and R. porosus or those that are part of fragile marine ecosystems, such as coral reefs inhabited by A. rivasi. For each species, La Guajira's population should be considered a genetic management unit (GMU) to prioritize implementing conservation and fisheries management measures [117]. Two marine protected areas are established in La Guajira (Bahía Portete – Kaurrele National Natural Park and Los Flamencos Sanctuary), but they are not sufficient for conservation.

On the other hand, additional C. pica and A. rivasi populations were identified in the Samarian subarea, which should also be treated as a second GMU. This sector has the advantage of including marine protected areas such as Tayrona Park and Salamanca Island Road Park, where the main factors affecting species and ecosystem conservation, such as fishing and tourism activities, are regulated. The population associated with the Isthmian subarea is the largest and should be considered the third GMU for both species. Although the marine protected areas of Corales del Rosario and San Bernardo islands and Acandí are established there, the challenge for the fishery and environmental authorities is to improve the regulation of multiple human activities that affect the conservation of species and ecosystems outside these protected areas.

Finally, additional studies should be conducted to investigate the phylogeographic patterns of each species assuming all areas of its distribution. The results will serve as a basis for implementing a multinational approach to developing conservation strategies in which the countries that host the species will be committed to implementation.

5. Conclusions

This study demonstrated for the first time the phylogeographic break caused by MRP for A. rivasi and the combined effect caused by the absence of the rocky bottom along more than 300 km of coastline and the permanent upwelling in La Guajira (ARB+PUG), which operates on A. rivasi and C. pica. Three populations (K = 3) were identified for A. rivasi and C. pica, while N. tessellata presented one population and exhibited the panmixia model (K = 1). Only A. rivasi and C. pica showed phylogeographic congruence. This could be due to the oceanographic and environmental conditions of the Colombian Caribbean generating discordance over species with short larval phases, specifically in the biological traits of these species (e.g., PLD). On the other hand, species with long larval phases do not show evidence of phylogeographic breaks.