Influence of the Great Amazon Reef System and Pleistocene sea-level drops on the phylogeography of Haemulon aurolineatum (Haemulidae)

Sample collections and data gathering

Sixty-five individuals of Haemulon aurolineatum were collected across the GC, with an additional 36 sequences downloaded from GenBank (http://www.ncbi.nlm.nih.gov/genbank/) and BOLDSYSTEMS (http://www.boldsystems.org/). In total, 101 individuals from 14 locations were included, covering the entire GC and representatives from each biogeographic province as defined by Robertson & Cramer (2014), as well as the Brazilian province (Fig. 1; Tables S1, S1.1 & S1.2).

Most samples were obtained with the assistance of fishermen while others were collected using a pole spear. We took a fin clip from each individual and preserved it in absolute ethanol at −70 °C. Tissue samples and voucher specimens were deposited in the fish collection at the Universidad Michoacana de San Nicolás de Hidalgo (CPUM) and the Fish Collection at the Museo de Zoología of the Universidad de Costa Rica. Organism collection was supported and allowed by the following institutions and permits by Mexican Commission of Fisheries and Aquaculture under collection permits (CONAPESCA-PPF/DGOPA-/2013 and CONAPESCA-PPF/DGOPA-262/17); for Costa Rica collection permits R-SINAC-SE-DT-PI-029-2023 and the Comisión Institucional de Biodiversidad of the Universidad de Costa Rica (RESOLUCIÓN No. 377); and for Colombia, Universidad Nacional de Colombia, collection permits 76991-2024 ANLA-Unal-MinAmbiente and 3914-27.

Sequencing, alignment and substitution models

We extracted genomic DNA of the 65 samples collected using the Phenol-Chloroform protocol (Sambrook, Fritsch & Maniatis, 1989). We amplified fragments of one mitochondrial and three nuclear DNA markers: Cytochrome oxidase subunit 1 (coxI), ribosomal protein S7 (S7), rhodopsin (Rho) and myosin (Myh). PCR reactions were performed with a final volume of 12.5 µl (Table S2). The amplicons were purified with Exo-sap enzymes and sent for sequencing to Macrogen Inc., Seoul, South Korea. The sequences obtained for each molecular marker were visualized, edited, and manually aligned using MEGA v. 10.2.2 software (Tamura, Stecher & Kumar, 2021). To determine the best substitution model for each molecular marker, we used the program jModeltest v2.1.5 (Darriba et al., 2012) with the corrected Akaike Information Criterion (AIC). To detect and resolve heterozygous individuals, we used PHASE module from DnaSP v. 6 (Rozas et al., 2017). Analyses were conducted independently for each molecular marker. All sequences were deposited on GeneBank (Accession numbers: coxI: PQ571864 –PQ571924; S7: PQ588478 –PQ588521; Rho: PQ588522 –PQ588559; Myh: PQ588560 –PQ588601). The coxI mitochondrial marker was chosen due to the known variability in other population studies (Torres-Hernández et al., 2022; Bernal-Hernández et al., 2024; Torres-García et al., 2024) and the possibility to find a high number of sequences in the public repositories (e.g., Brazilian samples). The nuclear markers were chosen to try to have a high arrange of mutation rate, and the possibility to detect both, recent and ancestral events, as have been show in previous studies (Gaither et al., 2011; Sandoval-Huerta et al., 2019; Palmerín-Serrano et al., 2021).

Genetic diversity and haplotype networks

We calculated the number of haplotypes (Hn), haplotypic diversity (h), nucleotide diversity (π), and the number of polymorphic sites (S) of mtDNA and nDNA by biogeographic province in ARLEQUIN v. 3.5 (Excoffier & Lischer, 2010). In order to analyze the genealogical relationships between haplotypes and their geographical correspondence, we constructed a statistical parsimony network for each molecular marker using the software PopArt with the Median–Joining algorithm (Bandelt, Forster & Röhl, 1999).

Genetic structure and divergence time estimates

We conducted an analysis of molecular variance (AMOVA) to assess genetic differentiation. The groups were tested as follows: (1) a panmictic population (single genetic group), (2) by GC biogeographic provinces (Northern vs. Central vs. Southern Caribbean) (Robertson & Cramer, 2014) plus the Brazilian province (for coxI only), (3) the Northern Caribbean vs. Central, Southern Caribbean and Brazilian provinces (SCB-Provinces), and (4) between haplogroups obtained from the statistical parsimony haplotype network. Furthermore, we estimated population genetic differentiation by computing pairwise Fst values (1) among biogeographic provinces, (2) the Northern Caribbean compared to Central, Southern Caribbean and Brazilian province) and (3) between haplogroups obtained from the statistical parsimony haplotype network. Both analyses were performed using ARLEQUIN v. 3.5 (Excoffier & Lischer, 2010), with each analysis including 1,000 permutations to estimate significance. We also calculated uncorrected p-Distances (p-D) for all markers based on the same grouping scheme as the Fst values using MEGA v. 10.2.2 software (Tamura, Stecher & Kumar, 2021).

Divergence times were estimated in BEAST v. 1.8.4 software (Drummond & Rambaut, 2007) through a Bayesian inference phylogenetic reconstruction. To calibrate the molecular clock, we used seven species from the closely related genus Anisotremus, and four from Haemulon (Table S3) (Tavera, Acero & Wainwright, 2018). Due to the availability of gene markers in the outgroups, only the mtDNA marker coxI and nDNA marker S7 were included in the analysis. Fossil ages were sourced from the Paleobiology Database (https://paleobiodb.org/). The first calibration point was based on the split between Haemulon and Anisotremus estimated at 23–16 Mya (Million years ago) (mean: 0.2, SD: 0.8 offset: 16; calibration: Haemulidae) (Aguilera & De Aguilera, 2004). The second calibration was in the crown-group of Haemulon in the late Miocene: 11.6–7.24 Mya (mean: 0.9, SD: 0.4 offset: 7.24; calibration: Haemulon spp.) as reported by Gillette (1984). The third calibration point was in the crown-group of Anisotremus, with Pleistocene ages 5.33–3.6 Mya (mean: 0.2, SD: 0.1 offset: 3.6; calibration: Anisotremus spp.), as reported by Fitch & Lavenberg (1983). The analyses were run using an uncorrelated relaxed clock model with a lognormal Yule speciation process. Three independent MCMC runs were performed, each consisting of 10 million generations, sampling every 1,000 generations. Convergence of the chains was assessed by discarding the first 10% of generations as burn-in, using Tracer v 1.7 (Rambaut, Drummond & Suchard, 2018). Estimated parameters from the independent runs were pooled using the LogCombiner module in BEAST, and the maximum clade credibility tree was generated with the TreeAnnotator module. The analyses were executed on the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2010). Finally, the tree was visualized and edited on FigTree v 1.4.2 (Rambaut, 2014).

Isolation by distance

We estimated patterns of isolation by distance (IBD) using coxI gene throughout a Mantel test (Mantel, 1967) with 10,000 permutations in the “vegan” library on R (R Studio Team, 2022; Oksanen et al., 2022). To perform the correlation between genetic and geographic distance, we use linearized Fst values (Fst/1-Fst) (Rousset, 1997) as genetic distances, whereas the geographic distances between localities were obtained using the package “fossil” on R (Vavrek, 2011).

Species delimitation

We conducted two species delimitation analyses using four molecular markers. First, we inferred a species tree using BEAST v1.8.4 (Drummond, Rambaut & Suchard, 2016). For this approach, we designated Haemulon vittatum as the outgroup species (Tavera, Acero & Wainwright, 2018). The analysis employed an uncorrelated relaxed clock model with a normal Yule speciation process. The Markov chain Monte Carlo (MCMC) runs consisted of 30 million generations, sampling every 3,000 generations. To ensure convergence, we discarded the first 10% of generations as burn-in and assessed convergence using Tracer v1.7 (Rambaut, Drummond & Suchard, 2018). Second, we performed a Species Tree and Classification Estimation (STACEY v1.3.1; Jones, 2017) analysis using the STACEY template in BEAST v2.6 (Bouckaert et al., 2019). STACEY infers a “Species or Minimal Clusters” (SMC) tree under the birth-death-collapse tree prior, which does not require a guide tree (Petzold & Hassanin, 2020). This analysis was conducted using a strict clock model and a birth-death-collapse tree prior with the following parameters: CollapseHeight = 0.0001; CollapseWeight = 0.5; BirthDiffRate = 100; RelativeDeathRate = 0.5; OriginHeight = 100. The MCMC runs for this analysis also consisted of 30 million generations, with results saved every 30,000 generations. Convergence and ESS values were evaluated using Tracer v1.7 (Rambaut, Drummond & Suchard, 2018). The resulting files were processed in the Species Delimitation Analyser with a collapse height of 0.0001 and a similarity cut-off of 1.0, following the exclusion of the first 10% of generations as burn-in.

Sequences and genetic diversity

For the mitochondrial gene coxI we obtained 97 sequences with a length of 628 base pairs (bp), exhibiting 43 polymorphic sites and 39 haplotypes. For S7, we obtained 44 sequences with a length of 387 bp, showing 46 polymorphic sites and 34 haplotypes. Regarding Rho, dataset contained 38 sequences with a length of 806 bp, featuring 17 polymorphic sites and 21 haplotypes. The dataset for the Myh gene consisted of 42 sequences with a length of 692 bp, displaying 34 polymorphic sites and 25 haplotypes. Samples from Brazilian province for the nDNA markers were not included in this study. The coxI gene exhibited the highest values of nucleotide diversity and haplotypic diversity in the Central Caribbean region and the lowest in the Brazilian province. The nuclear gene S7 showed the highest values of nucleotide diversity and haplotypic diversity in the Southern Caribbean region, and the lowest in the Northern Caribbean. In Rho, the highest values of nucleotide diversity and haplotypic diversity were found in the Southern Caribbean region, and the lowest in the Northern Caribbean. Lastly, in the Myh the highest haplotypic diversity was found in Southern Caribbean, the highest nucleotide diversity in Northern Caribbean and the lowest in the Central Caribbean (Table 1). We acknowledge that the absence of nDNA markers for the Brazilian province could introduce bias into our results. However, we believe the observed pattern of genetic homogeneity in the coxI Central and Southern Caribbean, as well as in the Brazilian province, is highly likely to persist even if nDNA markers are included in the analyses. This conclusion is supported by our own findings, where similar patterns of variation were observed in both mtDNA and nDNA data across the broader sampling area and provinces. Additionally, since mtDNA is known to exhibit greater variability than nDNA (see Vawter & Brown, 1986; Allio et al., 2017), we expect this pattern to hold. Nevertheless, we strongly recommend including nDNA markers in future studies to enhance the robustness of the findings.

Haplotype networks

The coxI haplotype network shows two haplogroups (Hg1 and Hg2) separated by eight mutation steps. Although Hg1 was composed of most individuals from the Northern Caribbean, also included five individuals from Central and one from Southern Caribbean, most of them showing peripheral position separated by three mutation steps from the nearest Northern Caribbean sample. Hg2 is composed by most of the individuals from the CSB-provinces, as well as two from Florida (St. Lucie County) and one from Yucatán (Sisal), mainly mixed with haplotypes of central province. The nDNA S7 marker shows more variation, with almost each sample forming their own haplotype without mixing, a somewhat geographic segregation is observed for the Northern Caribbean. The nDNA Rho and Myh genes show the haplotypes from the Northern Caribbean mostly segregated from those of the Southern and Central Caribbean samples (Fig. 2).

Genetic structure and divergence time estimates

The results of the AMOVA for a single population found a high and significant ϕ_ST, rejecting the null hypothesis of panmixia. The haplogroup arrangement (Hg1 vs. Hg2) maximized variation among groups for the coxI gene (Table 2). In contrast, for the nuclear gene S7 the highest variation was observed within populations. While Rho found the highest variation among populations within groups. Finally, for Myh, the variation was maximized within populations. Fst values as well as p-D show the same pattern as the haplotype networks in all molecular markers, where the Northern Caribbean is well differentiated in relation to the CSB-provinces (Fst = 0.58 to 0.80 (p < 0.05); p-D = 0.42 to 1.89; Table 3). The split of both, the Northern Caribbean and CSB-Provinces H. aurolineatum lineages was dated ca. 0.8 Mya (95% HPD: 0.35−1.38 Mya) (Fig. 3).

Isolation by distance

Mantel test shows a non-significant correlation between genetic and geographic distances (r = 0.0413; P = 0.3551).

Species delimitation

The species tree supported the separation of the Northern Caribbean from the Central and Southern Caribbean clusters with a posterior probability of one (Fig. 3). The results from STACEY for the global test, where we defined three minimal clusters (Northern vs. Central and Southern Caribbean vs. Outgroup), according to the hypothesis obtained in both phylogenetic reconstructions, showed high support (pp = 0.99).

Divergence of the two genetic groups

All the analyses conducted here reveal two well-supported genetic groups, one mainly distributed in the Northern Caribbean and the other in the CSB-Provinces, with genetic p-D of 1.89% and divergent event calculated to occur at 0.8 Mya. Nevertheless, previous studies in Haemulidae family provided older estimates (Tavera, Acero & Wainwright, 2018), placing the split between the two H. aurolineatum lineages at 2 Mya. The discrepancies between estimations could be associated with the used calibration points. Tavera, Acero & Wainwright (2018) used fossils of the stem lineage of Siganidae and Scatophagidae with an estimation date between 63.9 to 55.8 Mya, a stem lineage of the subfamily Acanthurinae with an estimation date between 57.3 to 50 Mya and a stem lineage of Luvaridae and Acanthuridae with an estimation date between 63.9 to 55.8 Mya. These three calibration points are concentrated in the Eocene-Paleocene epoch, which may bias the time estimates towards a more distant past (see: Magallón, 2010; Molak et al., 2012; Duchêne, Lanfear & Ho, 2014). On the contrary, in the present study we used three Haemulidae fossils, including the stem group of Haemulon and Anisotremus with dates between 23 to 16 Mya, the crown group of Anisotremus with dates between 5.3 to 3.6 Mya and the crown group of Haemulon with dates between 11.6 to 7.2 Mya (Palmerín-Serrano et al., 2021). These calibration points are better distributed over time, covering a window from 23 to 3.6 Mya. This broader time span in the calibration points enhances the accuracy of dating recent divergence events and helps mitigate temporal biases toward older estimates (Parham et al., 2012; Molak et al., 2012; Duchêne, Lanfear & Ho, 2014), as observed in the case of the two distinct genetic populations of H. aurolineatum. Consequently, we consider the dating findings in the present study to be more accurate than those reported by Tavera, Acero & Wainwright (2018).

Since IBD analysis indicates a non-significant correlation between genetic and geographic distance, the results suggest that the divergent event must be influenced by geological, oceanographical or biological-ecological barrier. Other fish species also show the isolation between Northern Caribbean (GOM) and CSB-Provinces samples, but they vary in divergences time, being calculated at 0.19 Mya in Dormitator maculatus (Galván-Quesada et al., 2016); or in genetic p-D, being 1.86% in Awaous banana (McMahan et al., 2021), to 5.3% in Gobiomorus dormitor (Guimarães Costa et al., 2017). In most of these studies, the genetic break was attributed to the effect of the Loop Current acting as a soft barrier. Loop Current is a warm water current that originates at the Yucatán Channel, flows through the strait between the Yucatán Peninsula and Cuba, loops back, and ultimately pass through the strait between Florida and Cuba (Vukovich, 1988; Oey, Lee & Schmitz Jr, 2003). This current could transport larvae from the Caribbean through the Yucatan Channel to the East, preventing its entrance to the GOM (Richards et al., 1993). This hypothesis is also supported by most of the samples (two of three) from the Atlantic side of Florida, which belong to the CSB-Populations haplogroup (Hg2), instead of the Northern Haplogroup (Hg1). Our divergence estimation (800,000 years ago) is in accordance with the dates estimated for the Günz Glaciation occurred between 850,000 to 600,000 years ago (Corrêa, 2021). The drastic climatic oscillations between repeated cycles of glaciation during the Pleistocene (Corrêa, 2021) dropped the sea level between 60 to 120 m lower than today, and coastlines extended horizontally between 10 and 100 kilometers, leading to the fragmentation and alteration of marine habitats (Jackson, 1992; Greenstein, Curran & Pandolfi, 1998; Ludt & Rocha, 2015; Diester-Haass, Billups & Lear, 2018).

The combination of the sea level drop (exposing large areas of the Florida and Yucatán peninsulas), with the dramatic sea surface cooling in the Gulf of Mexico (GOM), the prevailing north-easterly winds over the Caribbean, the strong outflow of the Mississippi River and the topography of the Yucatán Peninsula, restricted the Loop Current entry into the GOM. Instead, the current veered directly towards the Florida Straits and Cuba (Schmidt, Spero & Lea, 2004; Nürenberg et al., 2008; Arellano-Torres, Amezcua-Montiel & Casas-Ortiz, 2023) (Fig. 4). As a result, the GOM became isolated, leading to the genetic separation of two evolutionary lineages in H. aurolineatum. Nevertheless, a comprehensive explanation of genetic differentiation between GOM and CSB-Provinces populations must consider not only their origins but also their maintenance under present conditions. Moreover, if we take into consideration the mixture of haplotypes in the Caribbean of Mexico and Belize as well as Florida, high rates of gene flow over a few generations can eliminate genetic differentiation between populations. Therefore, it appears that the oceanographic barriers to dispersal or selection against migrants (or hybrids) are likely acting to maintain geographic variation currently.

Homogeneity of the Central, Southern and Brazilian provinces

All our analyses across all genes show a pattern of genetic homogeneity in H. aurolineatum among the Central, Southern Caribbean and Brazilian provinces. Similar genetic homogeneity between the Central and Southern Caribbean provinces has been reported in other species, such as Abudefduf saxatilis, Sparisoma viride, Dormitator maculatus, Agonostomus monticola, and Awaous banana (McMahan et al., 2012; McMahan et al., 2021; Galván-Quesada et al., 2016; Piñeros & Gutiérrez-Rodríguez, 2017; Loera-Padilla et al., 2021). However, genetic homogeneity between the Southern-Central Caribbean and Brazilian provinces is an uncommon pattern in reef-dwelling fishes. Typically, this region displays genetic breaks, with the emergence of sister species pairs or populations among shallow-water reef fishes in the WTA, such as Ophioblennius macclurei and Ophioblennius trinitatis (Muss et al., 2001), Acanthurus bahanius (Rocha et al., 2002), Sparisoma frondosum—Sparisoma griseorubra (Robertson et al., 2006), Halichoeres radiatus—Halichoeres brasiliensis, Sparisoma axillare—Sparisoma rubripinne (Rocha et al., 2008a; Rocha et al., 2008b), Haemulon atlanticus, Coryphoptetus thrix, Enneanectes altivelis (Araujo et al., 2022), Centropomus irae (Carvalho-Filho et al., 2019)—Centropomus undecimalis (Malcher et al., 2023), among others. This commonly observed pattern of genetic differentiation is often attributed to the Amazon-Orinoco plume barrier, which assumed its current configuration during the early Pleistocene (∼2.4 Mya). Over time, the permeability of this barrier fluctuated with sea-level changes: it became less permeable during glacial periods, preventing species from crossing through the Great Amazon Reef System as the sea level dropped. In contrast, during interglacial periods, it became more permeable as the sea level rose, leading to the opening of the Great Amazon Reef System and allowing species to cross (Rocha, 2003; Ludt & Rocha, 2015; Araujo et al., 2022). Despite this, H. aurolineatum does not exhibit the typical genetic break seen in species like Acanthurus chirurgus (Rocha et al., 2002), Haemulon plumierii and Citharichthys spilopterus (De Jesus Gama-Maia et al., 2024).

To explain the observed genetic homogeneity of the tomtate grunt, we considered three possible scenarios. First, during interglacial periods, the increased sea levels may have made the Amazon-Orinoco plume barrier more permeable (Rocha, 2003; Araujo et al., 2022), allowing species like H. aurolineatum to cross, resulting in genetic homogeneity, as has been explained for Opisthonema oglinum (Ferreira-Araújo et al., 2024), Selene stepannis and H. aurolineatum (Araujo et al., 2022). Second, juveniles of H. aurolineatum, like other species in the Haemulidae family, may use mangrove and estuarine habitats during their life cycle (Castro-Aguirre, Pérez & Schmitter-Soto, 1999; Bravo, Eslava & González, 2009), potentially facilitating the exchange of individuals between the Greater Caribbean and Brazilian provinces. A similar pattern of connectivity across sandy and mangrove gaps has been observed in Anisotremus interruptus in the Tropical Eastern Pacific (Palmerín-Serrano et al., 2021). However, only two studies suggest that these species inhabit estuaries in their early stages, while most evidence points to H. aurolineatum being strictly marine (Darcy, 1983; Böhlke & Chaplin, 1993; Ornellas & Coutinho, 1998; Hoese & Moore, 1977; Robertson et al., 2023). This may be attributed to the findings of Castro-Aguirre, Pérez & Schmitter-Soto (1999) and Bravo, Eslava & González (2009), which report the presence of juvenile specimens in mangrove ecosystems within coastal lagoons, particularly in estuarine environments characterized by hypersalinity or negative salinity (Cervigón, 1986; Potter et al., 2010; Tweedley et al., 2019; Lasso-Alcalá et al., 2023). Thirdly, a continuous distribution and interchange of individuals through the Amazonian Coast plume may have occurred (including the Amazon-Orinoco plume zone). This region is characterized by large positive estuaries (Cervigón et al., 1992; Potter et al., 2010). This extends from the Orinoco River Delta (Venezuela) through the Guianas (Guyana, Suriname, and French Guiana) to the mouth of the Amazon River, and continues further east along the northern coast of Brazil, spanning throughout the states of Pará, Maranhão, and Piauí, up to the Parnaiba River Delta (see Fig. 1 and Fig. S1). This hypothesis is the most plausible, as data from the Global Biodiversity Information Facility (GBIF, 2024) shows a continuous distribution of the species (see Fig. S1). Notably, while H. aurolineatum is typically recorded in coastal reefs across its distribution range, in the Amazonian Coast region, the species has been recorded far from the coast and at more than 30–40 m deep (Lowe-McConnell, 1969; Collette & Rützler, 1977; Uyeno, Matsuura & Fujii, 1983; Guéguen, 2000; Moura et al., 2016; Pinheiro et al., 2018; Mrceniuk et al., 2021). These records from the Amazonian Coast region are aligned with a large area of reef patches known as the Great Amazon Reef System, which extends from eastern French Guiana to the north-eastern mouth of the Amazon River, off the coasts of Amapá, Pará, and part of Maranhão states in Brazil (Rosemary & McConnell, 1962; Collette & Rützler, 1977; Francini-Filho et al., 2018; Mrceniuk et al., 2021; Carneiro et al., 2022). This mesophotic reef (30–120 m depth), with its carbonate structures, acts as a corridor for genetic connectivity of reef-associated species across a wide range of depths (Cordeiro et al., 2015; Moura et al., 2016; De Mahiques et al., 2019; Banha et al., 2022; Carneiro et al., 2022). This reef habitat bridge appears to enable gene connectivity between H. aurolineatum populations in the Central-Southern Caribbean and Brazilian provinces, avoiding the outflow of the Amazon-Orinoco plume, as has been suggested for Coryphopterus venezuelae and species of the Halichoeres genus (Rocha, 2003; Volk et al., 2021).

Taxonomic implications

Phylogenetic studies of the Haemulidae family have revealed significant divergences within H. aurolineatum, identifying two distinct lineages (Tavera et al., 2012; Tavera, Acero & Wainwright, 2018). Taxonomically, two closely related species have been described: Haemulon aurolineatum, based on type specimens collected in San Domingo, Brazil (Cuvier, 1830), and Haemulon rimator, described from specimens collected in Charleston, South Carolina, and Pensacola, Florida (Jordan & Swain, 1884). Haemulon rimator was later considered a subspecies of H. aurolineatum by Ginsburg (1948) and subsequently synonymized by Courtenay (1961). However, Castro-Aguirre, Pérez & Schmitter-Soto (1999) treated it as a valid species. Most recently, Fricke, Eschmeyer & van der Laan (2024) again regarded H. rimator as a synonym of H. aurolineatum. Our population genetic analyses and species delimitation analyses confirm the presence of two well-differentiated lineages along the WTA, with all posterior probability values in the tree species delimitation analyses conducted of pp = 1−0.99, with p-D of 1.89% and a divergence event estimated to have occurred around 0.8 Mya. One lineage is predominantly distributed in the northern Caribbean, including samples near the known type locality of H. rimator, while the other is found mainly in the CSB-provinces, where H. aurolineatum was originally described (San Domingo, Brazil). According to Castro-Aguirre, Pérez & Schmitter-Soto (1999), these two species differ in the number of anal fin rays and coloration. They also noted differences in distribution, with H. aurolineatum ranging from Massachusetts to Brazil and H. rimator from Massachusetts to the Gulf of Mexico—findings that align closely with the distribution of the two lineages identified in this study. However, our species delimitation analyses must be taken with caution, since these methods are sensitive to sample size (Carstens et al., 2013). Accordingly, we recommend a comprehensive and integrative systematic analysis of H. aurolineatum across its distribution range, incorporating both morphological and genomic data at intra- and inter-regional spatial scales.