Evaluating the boundaries of marine biogeographic regions of the Southwestern Atlantic using halacarid mites (Halacaridae), meiobenthic organisms with a low dispersal potential

Abstract Aim We evaluated traditional biogeographic boundaries of coastal marine regions in Southwestern Atlantic using DNA sequence data from common, rocky‐shore inhabiting, marine mites of the genera Agauopsis and Rhombognathus, family Halacaridae. Methods We investigated geographic population genetic structure using CO1 gene sequences, estimated divergence times using a multigene dataset and absolute time‐calibrated molecular clock analyses, and performed environmental niche modeling (ENM) of common marine mite species. Results Agauopsis legionium has a shallow history (2.01 Ma) with four geographically differentiated groups. Two of them corresponded to the traditional Amazonian and Northeastern ecoregions, but the boundary between the two other groups was inferred at the Abrolhos Plateau, not Cabo Frio. Rhombognathus levigatoides s. lat. was represented by two cryptic species that diverged 7.22 (multilocus data) or 10.01 Ma (CO1‐only analyses), with their boundary, again at the Abrolhos Plateau. ENM showed that A. legionium has suitable habitats scattered along the coast, while the two R. levigatoides cryptic species differ considerably in their niches, especially in parameters related to upwelling. This indicates that genetic isolation associated with the Abrolhos Plateau occurred in both lineages, but for the R. levigatoides species complex, ecological niche specialization was also an important factor. Main conclusions Our study suggests that the major biogeographic boundary in the Southwestern Atlantic lies not at Cabo Frio but at the Abrolhos Plateau. There two biogeographically relevant factors meet (a) changes in current directions (which limit dispersal) and (b) abrupt changes in environmental parameters associated with the South Atlantic Central Waters (SACW) upwelling (offering distinct ecological niches). We suggest that our result represents a general biogeographic pattern because a barrier at the Abrolhos Plateau was found previously for the fish genus Macrodon (phylogeographic data), prosobranch mollusks, ascidians, and reef fishes (community‐level data).


| INTRODUC TI ON
Concordance of phylogeographic (population level) and biogeographic (above species level) patterns provides a useful framework for investigating and comparing patterns of evolution of faunas around the globe (Dawson, 2001). Such concordance is expected since faunal provinces are based upon endemism (Briggs, 1974), and the same factors (e.g., barriers to dispersal, divergent natural selection) that fuel speciation also lead to population genetic structuring, making both processes interchangeable in the short term (Nosil, Harmon, & Seehausen, 2009;Sukumaran & Knowles, 2017).
As stated by the founders of phylogeography, "macroevolution is ineluctably an extrapolation of microevolution" (Avise et al., 1987).
Counterintuitively, when associated with natural selection, strong genetic population structuring could be found even in species with long-term pelagic larvae, or even completely pelagic species. In these species, local adaptation ("local" often spans over hundreds of kilometers in the marine realm) could lead to speciation despite high immigration rates of poorly fit genotypes outside of the selective environment (Palumbi, 1992;Sanford & Kelly, 2011).
In the marine environment, breaking points in genetic structuring or community composition are commonly associated with: (a) hydrodynamics, such as currents, affecting recruitment and/or availability of planktonic dispersal stages (e.g., Francisco, Mori, Alves, Tambarussi, & Souza, 2018); (b) gradients in ecological factors that can impose a limit on geographic distributions or lead to local adaptation. However, they may be difficult to disentangle because changes in superficial circulation often create dispersal barriers where disparate water masses meet (Gaylord & Gaines, 2000;Pelc, Warner, & Gaines, 2009). Spalding et al. (2007)  Temperate Southwestern Atlantic, with its southern limit at Váldez Peninsula. These provinces were further subdivided into five ecoregions (areas of relatively homogeneous species composition, clearly distinct from adjacent systems) as shown on Figure 1a. This division of coastal and shelf areas of the Brazilian coast greatly influenced subsequent studies on marine conservation and biogeography (e.g., Bissoli & Bernardino, 2018;Huguenin et al., 2018;Mantelatto, Cruz, & Creed, 2018). Importantly, this biogeographic classification adopted an older F I G U R E 1 (a) Map of Brazilian coast showing 36 sampling localities, traditional boundaries of ecoregions (Spalding et al., 2007). Abrolhos Plateau and Vitória-Trindade ridge are highlighted in green, Cabo Frio indicated. (b) Southwestern Atlantic currents. Yellow arrows represent superficial warm, saline and oligotrophic currents, NBC: North Brazil Current; South Equatorial Current (SEC), with three branches: Southern (SSEC), Central (CSEC), and Equatorial (ESEC). Black arrows represent colder, less saline, and more nutritive currents: the South Atlantic Central Water (SACW, 100-500 m depth) and Malvinas (Falklands) current (MC) (a) (b) view that the main biogeographic break in the region is Cabo Frio, which is the southest limit for coral reefs, and is in rough agreement with the 23°C isotherm data (Briggs, 1974(Briggs, , 1995Palacios, 1982 Similarly, coastal fishes of the genus Macrodon form two almost parapatric sister species (Santos, Hrbek, Farias, Schneider, & Sampaio, 2006), separated at the Abrolhos Plateau region (approx. latitudes 15-20°S):

| Brazilian coastal biogeography and phylogeography
king weakfish M. ancylodon and Southern king weakfish M. atricauda (Carvalho-Filho, Santos, & Sampaio, 2010). Niche differentiation due to the SACW upwelling is the most plausible explanation for the presence of these two species, because, conceivably, these fishes can easily migrate over the Abrolhos Plateau region. In summary, overall data from nearly all available studies suggest that the Southern limit for the tropical fauna should be northward of the traditional boundary at Cabo Frio.

| Main currents and coastal waters
Marine biogeographic boundaries usually occur where major oceanic currents interact. In the ocean, colliding currents may create flowinducing boundaries affecting dispersal and distributional patterns near the shore (Gaylord & Gaines, 2000). Currents also determine the characteristics of coastal waters (CW), that is, the outcome of interactions among oceanic water masses and continental drainage  (Spalding et al., 2007) are given; population color-coding follows that on Figure 2, it indicated the Cabo Frio and the Abrolhos Plateau (Silveira, Miranda, & Brown, 1994;Silveira, Schmidt, Campos, Godoi, & Ikeda, 2000;Figure 1b). From there, the NBC meets the northern branch of the SEC, turning northwestward. Just a minor portion of the water transported westward by the Southern SEC branch turns south near 10°S, feeding the Brazil Current (BC).
Between this latitude and the Abrolhos Plateau, local circulation and the southernmost branch of SEC that hits the littoral at approximately 17°S make the coastal water transport complicated, at some points with almost no transport or others with northward transport (Peterson & Stramma, 1991).

| Temporal context
As mentioned above, the AP region (Abrolhos Plateau and Vitória-Trindade ridge) can affect passive transport of marine organisms (e.g., Francisco et al., 2018; for theoretical models for colliding oceanic currents creating distributional patterns, see Gaylord & Gaines, 2000). At the same time, changes in the CW parameters would lead to local adaptation or to more or less permeable barriers. These changes also should take place at the AP region because the SACW influence limit occurs there. The presence of SACW is related to the presence of a fully operating cold Malvinas (Falkland) Current (MC), which in its turn is related to the opening of the Drake Passage. Despite superficial water circulation that may have started 36-28 Ma, the opening between South America and Antarctica became effective for deep waters not earlier than the middle Miocene (Crame, 1999). Martínez and del Río (2002), however, showed that mollusks from latitudes as south as 42°S were inhabitants of warm waters by the late Miocene (Tortonian, up to 7.25 Ma), suggesting that the MC still had a very low activity in the area during that time.
The minimal age for this system must correspond to the develop-

| Halacarid mites as models for marine biogeography
This study aims at testing long-held ideas about biogeographic division of the Brazilian coast using phylogeographic analyses of common marine meiobenthic species as proxies. Those species were chosen because they have several natural history features Note that the posterior probability for the NE and SE clades to be the same species is <0.05 maximizing divergence and maintaining genetic structure (Dawson, 2001): low fecundity, poor dispersal abilities, and intertidal habitat.
Marine mites (family Halacaridae) differ from most common marine benthic organisms by the lack of a dispersing planktonic stageimmature mites look similar to adults and live in the same environment; the mites have low fertility and are unable to swim at any stage of their life cycle (Bartsch, 2004). They rely on passive means of dispersal, like rafting on algae or hitchhiking on large marine organisms covered by algae and debris (Bartsch, 2004). Halacarids have one larval and one to three nymphal stages (proto-, deuto-, tritonymph) before they molt to the adult stage. Although parthenogenesis is known in the family, they usually have separate sexes, with females slightly outnumbering males in a given population (Bartsch, 2006); this is the case in the two linages, Rhombognathus and Agauopsis, included in our study. Most species have a univoltine life cycle, with either short or prolonged periods of reproduction (Bartsch, 2006). In general, the fecundity is low, especially when compared to other marine invertebrates. Free-living species, like those studied here, produce <20 eggs per female lifetime (Bartsch, 2004(Bartsch, , 2006. Because halacarid mites are globally distributed and are an ancient group diverging more the 300 Ma and radiating approximately 270 Ma from a terrestrial ancestor (Pepato, Vidigal, & Klimov, 2018), they can be used for inferring both ancient and recent biogeographic scenarios at local and global scales (Bartsch, 2004).
Halacarids are truly marine animals distinct from, for example, marine Oribatida, whose restriction to littoral zones suggests that F I G U R E 5 Time trees inferred from COI sequences under the multispecies coalescent model in *BEAST for Rhombognathus levigatoides and Agauopsis legionium. Numbers above branches are coalescence time estimates. Node bars represent 95% HPD intervals for time estimates. Vertical gray column indicates approximate period when the Malvinas (Falkland) current (MC) was established 7.25-4.8 Ma.
Inserts are the species trees, indicating the prior (from the multilous analyses) and the posterior on divergence ages. HG, Haplogroups as referred in Figure 2. Color code as in Figures 2 and 3 they have not completely transcended the ecological barrier between marine and terrestrial environments (Procheş & Marshall, 2001). Halacarids inhabiting the upper intertidal fringe apparently evolved tolerance to desiccation and environmental fluctuations characteristic of such habitats (Pugh & King, 1985;Somerfield & Jeal, 1995). The family is exclusively benthic, with species living in interstitial spaces, on marine algae, or sessile invertebrates. Four halacarid genera feed on macroalgae. These were traditionally included in a single subfamily, Rhombognathinae, which included the genus Rhombognathus (Abe, 1998). In recent molecular studies, this subfamily was shown to be diphyletic . In contrast to the algivorous Rhombognathus, species in another lineage included in our study, mites of the genus Agauopsis, were observed preying on copepods, ostracods and other halacarids, capturing and holding the prey between the femoral, genual and tibial crook and the spines of the first leg (Krantz, 1970;MacQuitty, 1984). Spalding et al.'s (2007) ecoregions are in rough agreement with some independently proposed divisions based on oceanographic data, for example, breaks between Northeastern and Eastern ecoregions at 10°S and Northeastern and Amazonian ecoregions at 2°S. However, this view is at odds with the boundary between the Eastern and Southeastern regions. Spalding et al. (2007) placed this boundary at Cabo Frio (23°S), but the oceanographic data suggest that this boundary lies at the Abrolhos region (16°S). Solving this controversy will be the main focus of our study. If the onset of current conditions took place at the Miocene-Pliocene transition, we expect that divergences of at least some sibling species could be relatively old, 4.8-7.5 Ma.

| Taxon sampling and sequencing
Sampling localities are reported in Figure 1 (Bartsch, 2000(Bartsch, , 2003 based on the fact that the differences between the two taxa are slight (e.g., metric characters related to the first pair of dorsal setae and sclerotization).
Irrespective of the status of R. levigatoides from Brazil, its populations lack morphological differences that would justify further splitting of this morphospecies.
To infer absolute time estimates of divergences of our target species (as detected by COI data) versus closely related species, we included sequences of three genes [COI, small (18S) and large (28S) subunit ribosomal genes] from previous studies (Dabert, Proctor, & Dabert, 2016;Klimov et al., 2018;Pepato & Klimov, 2015;) and a total of 23 terminals for which the sequences were generated as part of this study (Appendix S1: Table S1).

| Phylogeographic and population genetics analyses
Our COI alignment did not have indels or stop-codons, which are indicative of pseudogenes. Gene trees were inferred after choosing the best-fitting model and partition scheme in Partition Finder 1.0.1 (Lanfear, Calcott, Ho, & Guindon, 2012).
SAMOVA 2 was employed to test our biogeographic hypotheses (Dupanloup, Schneider, & Excoffier, 2002). Our analyses were based on 1,000 simulated annealing steps and the prior definition of the number of groups, K, ranging from 2 to 10. For each analysis with increasing K, we examined the proportion of genetic variance due to differences between groups, FCT, and searched the range of K for which FCT was the largest and statistically significant. We also calculated BIC (Bayesian information criterion) according to Meirmans (2012).
Additional AMOVA analyses were performed to compare our datadriven hypotheses with those constrained to reflect the Brazilian biogeographic provinces and ecoregions of Spalding et al. (2007).
Our COI phylogeny and haplotype networks (see below) suggested the presence of two distinct cryptic species in Rhombognathus levigatoides. To evaluate this result, we conducted a barcoding gap discovery analysis in ABGD (Puillandre, Lambert, Brouillet, & Achaz, 2012) and Generalized Mixed Yule Coalescent species delimitation analysis in bGMYC (Reid & Carstens, 2012).
For bGMYC analyses, all 99 haplotypes were included in a standard BEAST analysis employing the Yule process for modeling lineage divergence. Runs employing lognormal relaxed and strict clock models led to choosing the former due better fit according to AICM (ΔAICM = 6.18). We used a subsample of 500 BEAST stationary trees; for each of these subsampled trees, we calculated species membership probabilities in the R package bGMYC (available at https ://sites.google.com/site/noahm reid/software).

ABGD was run with the default values (web version), except by
Nb bins = 64 and using Kimura (K80) T s /T v ratio = 4.304 as calculated in MEGA7. The sequences were directly submitted to the web server instead of submitting a precalculated distance matrix.

| Absolute time estimation for mite divergences
A time calibration analysis was performed in BEAST (Bouckaert et al., 2014). We used a representative 169-taxon acariform taxonset, allowing time calibration using paleontological data (12 calibration points, Appendix S1: Table S2). Outgroup sampling, calibration points, methodological details, and results of time-calibrated analyses using multiple loci are reported in Appendix S1.
Time-calibrated COI gene trees were inferred in *BEAST, using the coalescent model and the species divergence priors informed by a separate 169-taxon, time calibration analyses based on fossils (see above, Appendix S1). For these secondary calibration points, we used lognormal priors (Morrison, 2008,
Maxent relies on the occurrence data to infer the niche parameters that shape species distributions (Phillips et al., 2006 Schoener's D metrics were used for niche comparisons (Schoener, 1970); this varies from 0 (no overlap) to 1 (complete overlap). Warren, Glor, and Turelli (2008) case) and RM is set to 1. Data were partitioned using the n-1 jackknife method, in which each occurrence is used for testing, while all others are used for training in that iteration (Muscarella et al., 2014).
The program was run iteratively: First, all layers selected by excluding highly correlated variables were used, but those variables with no contribution to the model along the range of settings examined were pruned. Models with the lowest AICc value and those differing by less than 2 (ΔAICc < 2) were considered as equivalent. If more than one set of parameters lie in this interval of AICc values, the choice was based on the maximum value of AUC test, that is, on our ability to predict the distribution of independent occurrences in the training area (Warren & Seifert, 2011). These two criteria are based on different theoretical approaches; the former is rooted in information theory and model fit to the data, and the latter is based on the performance of independent test data. We consider these criteria complementary to each other.
Moreover, we ran analyses taking into account the effect of sampling bias by feeding into Maxent a bias file in which the cell values reflect sampling effort and give a weight to random background data used for modeling (Elith, Kearney, & Phillips, 2010

| Agauopsis legionium phylogeography
The longest edge in the network (Figure 2)   In the preferred K = 4 scheme, 59.23% of the variation was observed among groups, 12.77% among populations within groups, and 28.00% within populations. This contrasts with a scheme con-  , there was no evidence suggesting that the observed correlation is related to isolation by distance.

| Rhombognathus levigatoides phylogeography and species delimitation
Rhombognathus species delimitation analyses with bGMYC and ABGD agree in inferring the two allopatric clades (see above) as separate species: There was a clear gap at p = .90-.12 attributed to the gap between interspecific and intraspecific distances (Figure 4a), and bGMYC analyses resulted in two distinct species, with the posterior probability for the alternative, single-species delimitation being p < .05 (Figure 4b). Hereafter, these species are referred to as NE and SE species of the Rhombognathus levigatoides species complex.

| Coalescent times among Rhombognathus and Agauopsis species haplotypes
Agauopsis legionium has a relatively recent history, originating 2.01, The Rhombognathus levigatoides species complex has a deeper history than that of A. legionium (Appendix S1, Figure 5). The diver-

| Species distribution modeling
Models accounting for sampling bias consistently had smaller AUC values and predicted distribution ranges broader as compared to models that do not incorporate sampling bias. These two main sets of models also pointed to different variables as main predictors.
The Northeast species of the R. levigatoides complex had its pre-

| D ISCUSS I ON
In the marine environment, two key aspects may influence the presence of genetic structuring or lack thereof: (a) hydrodynamics, such as currents, affecting recruitment and/or availability of passively transported propagules and (b) gradients in ecological factors that can impose a limit on geographic distributions (Andrade et al., 2011). The halacarid mite populations studied here do show a remarkable genetic structuring along the Brazilian coast. In contrast, a substantial number other coastal organisms inhabiting the same area show slight if any genetic structure, for example, mangrove crabs, Ucides cordatus (Oliveira-Neto et al., 2007), land blue crabs (Oliveira-Neto et al., 2008), interstitial ribbon worms, Ototyphlonemertes erneba and O. evelinae (Andrade et al., 2011) and littorinid snails, Nodilittorina lineolata and Littoraria flava (Andrade et al., 2003). Marine mites lack any long-distance dispersal stage and arguably do not extensively engage in passive dispersal by rafting on drifting algae (Bartsch, 2004). Therefore, the presence of considerable geographic genetic structuring in both species of marine mites is not surprising.
Populations of A. legionium and the two allopatric species of the

| Temporal context
The two species of the R. levigatoides complex diverged about 7.2, 1.69-16.28 Ma (Appendix S1, Figure 5) we note that the divergence of these two cryptic species (Appendix S1, Figure 5) matches well with the dates of the presumed origin of

| General considerations
It is remarkable that the characteristic allopatric distributions observed in the NE and SE species of the R. levigatoides complex are mirrored by the distributions of two king weakfish sibling species of the genus Macrodon (Santos et al., 2006). Similarly, many community-level studies (Barroso et al., 2016;Floeter et al., 2001;Lotufo, 2002;Pinheiro et al., 2018) also point to a break northward relative to Cabo Frio. These data indicate the presence of a general biogeographic pattern, most importantly highlighting the Abrolhos Plateau (AP) as a major biogeographic break and the SACW upwelling as major ecological factor affecting niche specialization.
Here, we suggest modifying currently practiced sampling strategies to include extensive sampling near the AP region and using more sophisticated models for circulation and environmental parameters. Special care should be taken when planning sampling along the Bahia State coast because its southern portion might be the boundary between tropical and subtropical faunas. Furthermore, our study highlights the importance of the complex current circulation pattern at the Abrolhos Plateau and Trindade Ridge (Peterson & Stramma, 1991), which may act as a physical barrier restricting gene flow between populations. This complex pattern is usually ignored in most marine biogeographic studies (e.g., Andrade et al., 2011;Santos et al., 2006).

| CON CLUS IONS
Our study supports the hypotheses that the main break in the distribution of coastal rocky-shore communities occurs at Abrolhos Plateau and Vitória-Trindade Ridge, not at Cabo Frio as traditionally accepted in the literature (Briggs, 1974(Briggs, , 1995Petuch, 2013;Spalding et al., 2007). This region is not only a barrier to gene flow due complex current patterns, but it also is a place where sharp changes in major environmental parameters associated with the

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
ARP, THDAV, and PBK conceived the project and conducted the fieldwork; ARP was mainly responsible for obtaining sequences and analyses; and ARP and PBK wrote the paper with the assistance from THDAV.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https ://doi.org/10.5061/ dryad.80gb5 mkn2.

DATA AVA I L A B I L I T Y S TAT E M E N T
• DNA sequences: Genbank accessions, MH999551-MH999800.