Population genetic structure, introgression, and hybridization in the genus Rhizophora along the Brazilian coast

Abstract Mangrove plants comprise plants with similar ecological features that have enabled them to adapt to life between the sea and the land. Within a geographic region, different mangrove species share not only similar adaptations but also similar genetic structure patterns. Along the eastern coast of South America, there is a subdivision between the populations north and south of the continent's northeastern extremity. Here, we aimed to test for this north‐south genetic structure in Rhizophora mangle, a dominant mangrove plant in the Western Hemisphere. Additionally, we aimed to study the relationships between R. mangle, R. racemosa, and R. × harrisonii and to test for evidence of hybridization and introgression. Our results confirmed the north‐south genetic structure pattern in R. mangle and revealed a less abrupt genetic break in the northern population than those observed in Avicennia species, another dominant and widespread mangrove genus in the Western Hemisphere. These results are consistent with the role of oceanic currents influencing sea‐dispersed plants and differences between Avicennia and Rhizophora propagules in longevity and establishment time. We also observed that introgression and hybridization are relevant biological processes in the northeastern coast of South America and that they are likely asymmetric toward R. mangle, suggesting that adaptation might be a process maintaining this hybrid zone.

tolerance, specialized pneumatophore roots, and seawater-based seed or fruit (propagule) dispersal (Ball, 1988;Duke, Ball, & Ellison, 1998;Tomlinson, 1986). Additionally, at the molecular phenotype level and, more specifically, at the transcriptome level, there are exceptional similarities in the gene expression profiles of mangrove species. These similarities suggest the occurrence of parallel evolution in two distantly related mangrove lineages to cope with a common environment (Dassanayake, Haas, Bohnert, & Cheeseman, 2009).
The resemblances among distantly related mangrove species are not limited to adaptations to an environment that is influenced by both land and sea. Different species also share similar genetic structure patterns at a wide range of geographic scales in many mangrove forest regions. Such similarities are apparent in the Indo-West Pacific biogeographic region (IWP), an area rich in mangrove species in the Eastern Hemisphere (Duke, Lo, & Sun, 2002;Duke et al., 1998). In South-East Asia, there is a clear divergence between populations from the western and eastern coast of the Malay Peninsula observed for Ceriops tagal (Perr.) C.B. Rob. (Rhizophoraceae) (Huang et al., 2012), C. decandra (Griff.) W. Theob. (Tan et al., 2005), Bruguiera gymnorrhiza (L.) Lamk. (Rhizophoraceae; Minobe et al., 2010), Excoecaria agallocha L. (Euphorbiaceae; Zhang et al., 2008), Lumnitzera racemosa Willd. (Combretaceae; Li et al., 2016), Rhizophora apiculata Blume (Rhizophoraceae; Yahya et al., 2014;Yan, Duke, & Sun, 2016), and R. mucronata Lam. (Wee et al., 2014). Although this pattern is not universal-for instance, it was not detected in other Rhizophora (Yan et al., 2016) and Sonneratia  species-it indicates that common extrinsic factors, such as superficial ocean currents (Wee et al., 2014) and sea-level fluctuations (Yan et al., 2016), have shaped the organization and distribution of genetic diversity.
In the Atlantic East Pacific region (AEP), a biogeographic region with lower mangrove species diversity than the IWP, a similar phenomenon was observed at a comparable geographic scale. Along the eastern coast of South America, there is a pattern of genetic subdivision between populations north and south of the northeastern extremity of the South America continent (NEESA). This divergence pattern is shared by the sea-dispersed plants Rhizophora mangle L. (Pil et al., 2011;Takayama, Tamura, Tateishi, Webb, & Kajita, 2013), Avicennia germinans L., A. schaueriana Stapf and Leechman ex Moldenke (Acanthaceae;  and Hibiscus pernambucensis Arruda (Malvaceae; Takayama, Tateishi, Murata, & Kajita, 2008). As in the IWP, the regional genetic structure in the AEP can be explained by superficial marine currents and responses to Pleistocene climate variations (Mori, Zucchi, Sampaio, & Souza, 2015;Pil et al., 2011;Takayama et al., 2008). Despite the genetic structure shared by distantly related species, restricted geographic sampling and variation among different molecular marker sets can potentially influence the description of intraspecific genetic structure.
Here, we tested for the north-south genetic structure pattern in this biogeographic region using R. mangle, a widespread, dominant mangrove species in the AEP biogeographic region, as the biological system. To assess the repeatability of previous results (Pil et al., 2011;Takayama et al., 2013), we developed a new set of microsatellite markers of the same class used previously to describe this genetic structure in this system. Additionally, we sampled plants from four and seven mangrove forest regions north and south of the NEESA, respectively, to improve the representation of the northern "population." Finally, because AEP Rhizophora species exhibit semipermeable species boundaries (Cerón-Souza et al., 2010Takayama et al., 2013), it is possible that ancient and ongoing interspecific hybridization is a relevant evolutionary process in the northern population (Pil et al., 2011). As the existence of such a process could influence the observed patterns of genetic structure, we also tested for the presence of a hybrid zone. In some areas to the littoral north of the NEESA, the three AEP Rhizophora species R. mangle, R. racemosa Meyer, and R. × harrisonii Leechman (pro sp.) are sympatric, although the distributions of the last two species are disjunct (Menezes, Berger, & Mehlig, 2008). Additionally, to the best of our knowledge, R. racemosa and R. × harrisonii have not been recorded south to the NEESA.

| Plant material and sampling strategy
We distinguished the species in the field primarily based on reproductive morphological traits to minimize misidentification issues.
At one extreme, R. mangle inflorescences exhibit 2-5 flowers each F I G U R E 1 Geographic distribution of Rhizophora samples along the northwestern South American coast. Acronyms indicate the locations where plants morphologically identified as R. mangle, R. × harrisonii, and R. racemosa were sampled, according to Table 1. Circles, squares, and triangles represent sampling sites where only R. mangle was collected, where only R. racemosa and R. × harrisonii were sampled, and where the three species were obtained from, respectively. The average speeds of marine currents based on The Global Drifter Program (US National Oceanic and Atmospheric Administration-Atlantic Oceanographic & Meteorological Laboratory) are represented using arrows, whose slope indicates current direction, while the magnitude indicates current velocity. Additionally, we represent current velocity with colors, with higher temperatures indicating higher velocities for each cell and no more than two orders of bifurcation. At the other extreme, R. racemosa inflorescences present up to 128 flowers each, with up to seven orders of bifurcation (Cerón-Souza et al., 2010;Tomlinson, 1986). Furthermore, the latter species displays rounded flower buds, whereas the former shows pointed ones. R. × harrisonii identification was based on its three to five ordered branched inflorescences with up to 32 flowers and its apically pointed flower buds. This species has a morphology intermediate between those of R. racemosa and R. mangle (Tomlinson, 1986), which is likely the outcome of ancient and ongoing gene flow between R. racemosa and R. mangle morphotypes (Cerón-Souza et al., 2010;Takayama et al., 2013).
We collected visually healthy leaf samples from 318 specimens of R. mangle from 11 sites and 33 samples of R. racemosa and 37 samples of R. × harrisonii from two localities along the Brazilian coast, covering more than 4,900 km of coastline ( Figure 1 We deposited voucher specimens from every location in the University of Campinas (UEC) and Embrapa Amazônia Oriental (IAN) herbaria, which are both in Brazil. For genetic analyses, we sampled leaves from flowering trees that were at least 20 m from any other sampled trees and maintained the leaves in resealable zipper plastic bags containing silica gel. The leaf material was lyophilized and stored at −20°C prior to DNA isolation. for 1 min and 72°C for 2 min; and 72°C for 5 min. The amplified samples were genotyped by vertical electrophoresis using 6% denaturing polyacrylamide gels, and the bands were visualized using silver nitrate (Creste, Neto, & Figueira, 2001). The sizes of the resulting fragments were estimated by comparison with a 10 bp DNA ladder (Thermo-Invitrogen, Carlsbad, CA, USA).

| Hybrid detection and characterization
To maximize the number of evaluated loci, we considered different microsatellite sets for each analysis considering their amplification Due to the complex hybridization and introgression processes We used the Bayesian method implemented in STRUCTURE 2.3.4, assuming a model with correlated allele frequencies and admixture (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000). For each number of groups (K) considered, which ranged from 1 to 10, we performed 30 independent Markov Chain Monte Carlo (MCMC) runs with 500,000 iterations, following a burn-in period of 100,000 steps. We determined the K value that best explained our data at the uppermost hierarchical level, considering the log likelihood of each K (lnL) (Pritchard et al., 2000) and the ad hoc statistic ΔΚ (Evanno, Regnaut, & Goudet, 2005). Postprocessing of this model- To identify and classify eventual two-generation hybrids as F 1 s, F 2 s or backcrosses, we used a different model-based method implemented in NewHybrids 1.1 (Anderson & Thompson, 2002). We performed 10 independent MCMC runs with 500,000 steps following a burn-in period of 100,000 iterations for different prior distribution combinations.
We performed DAPC, STRUCTURE, and NEWHYBRIDS analyses considering all samples and those from the north coast of Brazil, because of the observed genetic structure pattern at its uppermost hierarchical level (see Results). Additionally, to evaluate the patterns of introgression between the clusters or hybrid classes inferred by these methods, we employed multinomial regression to estimate individual hybrid indices, corresponding to individual clines in genotype frequency for each marker along a genomic admixture gradient (Gompert & Buerkle, 2009). This regression was implemented using the INTROGRESS R package (Gompert & Alex Buerkle, 2010). As we obtained similar, but partially contrasting results between the STRUCTURE analysis and the DAPC and NEWHYBRIDS analyses (see Results), we estimated the hybrid indices for individuals collected from regions of sympatry that were identified as hybrids using at least one of these methods. As nonadmixed individuals, we considered R. mangle with STRUCTURE-estimated ancestry values (Q-values) higher than 0.80 for the typical R. mangle inferred cluster for the most likely scenario, K = 2 (see Results) and individuals identified as R. racemosa in the field (that were assigned to the alternative cluster for the K = 2 scenario).
Microsatellite alleles were combined into classes of alleles with frequency differentials between the nonadmixed groups (i.e., species) equal to those obtained when each allele was considered separately.
This approach reduces data complexity without distorting the genetic similarity between species or losing information (Gompert & Buerkle, 2009 T A B L E 3 Analysis of molecular variance (AMOVA) for the a priori and a posteriori hypotheses regarding the Rhizophora species complex Finally, due to the complex pattern of hybridization and introgression among these Rhizophora species, we used Arlequin 3.5 (Excoffier & Lischer, 2010) to conduct hierarchical analysis of molecular variance (AMOVA) (Excoffier, Smouse, & Quattro, 1992). We tested competing a priori hypotheses based on the morphological identification of individuals and mangrove genetic structure patterns in the Neotropics (Mori, Zucchi, Sampaio, et al., 2015;Pil et al., 2011;Takayama et al., 2013). Additionally, we tested the a posteriori hypothesis that R. mangle is a separate group and that R. racemosa and R. × harrisonii constitute a different group (see Results) as well as the hypothesis of an additional subdivision within R. mangle (Table 3). As criteria for identifying the scenario that was best supported by our data, we employed the maximum variance among groups (Φ CT ) and a significant departure from a random distribution obtained after 10,000 permutations.

| Intragroup genetic structure
After describing the genetic structure at the interspecific level, we performed intraspecific genetic structure analyses of each inferred group (see Results). As we used a different dataset for each inferred group to maximize the genetic information evaluated, for each group, we investigated the presence of LD and microsatellite technical artifacts (stuttering, large allele dropout and null alleles) using the previously described methods (Excoffier & Lischer, 2010;Van Oosterhout et al., 2004). Again, we did not observe consistent evidence of LD nor technical artifacts across the sampled populations; therefore, we used the entire datasets for the subsequent analyses.
To describe how the observed genetic variation was organized, we applied DAPC, using BIC to identify the K value that best explained our data and the optim.a.score function to avoid overfitting (Jombart et al., 2010). Additionally, we applied the clustering method implemented in STRUCTURE 2.3.3, assuming a model with correlated allele frequencies and admixture (Falush et al., 2003;Pritchard et al., 2000) and considering K values from 1 to 10. Thirty independent MCMC runs with 500,000 iterations and a burn-in period of 100,000 was performed. The K value that most efficiently summarized our data was determined based on lnL (Pritchard et al., 2000) and ΔΚ (Evanno et al., 2005), and the STRUCTURE results were postprocessed using CLUMPAK (Kopelman et al., 2015).
For the R. mangle group, we also applied spatial principal component analysis (sPCA; Jombart, Devillard, Dufour, & Pontier, 2008), implemented in ADEGENET 2.0 (Jombart & Ahmed, 2011), to study the genetic variability in its spatial context. sPCA detects spatial genetic structure patterns that simultaneously present spatial autocorrelation and high variation without assuming Hardy-Weinberg equilibrium or linkage equilibrium among loci (Jombart et al., 2008). It also allows for identifying "positive" or "negative" autocorrelation, which respectively emerge from genetic similarity or spatial gradients that result from barriers to gene flow, isolation by distance or adaptation ("global structure") and greater genetic differences among neighbors than expected ("local structure") (Jombart et al., 2008). To test for global and local structure, we performed 10,000 permutations. We determined the connection network using the minimum distance that retained all sampling localities in this network.

| New microsatellite isolation and characterization
We  Figure S1).

| Hybrid detection and description
Using seven microsatellite markers (Table 1) the simulated results, we observed that our experimental results best fit the F 2 and backcross simulated cases (Figure 3). These results suggest that contemporary hybridization between nonadmixed parental species is uncommon. Additionally, they indicate that F 1 hybrids are viable and that generally, two or more generations have passed since hybridization started.
F I G U R E 2 Hybridization and historic introgression in the Rhizophora species complex considering the entire sample (a) and only individuals from the northern coast of Brazil (b). From left to right: Scatterplot of the discriminant analysis of principal components (DAPC), where individuals morphologically identified as R. mangle, R. × harrisonii, and R. racemosa are represented as blue circles, red squares, and red triangles, respectively (a), or by colors (b). Bar plots showing STRUCTURE 2.3.4 Bayesian group assignments, in which each bar represents a single individual and colors represent proportional membership coefficients, for two (K = 2), three (K = 3) and four (K = 4) groups, considering the whole sample, or for K = 2 and K = 3, considering only the northern samples. Bar plots depict the NEWHYBRIDS 1.1 Bayesian hybrid class assignment, in which bars represent individuals; black represents one "pure species," and gray represents the alternative "pure species (NH)." F I G U R E 3 Triangle plot summarizing interspecific heterozygosity against the hybrid index for putative, inferred and simulated Rhizophora hybrids. The hybrid index was measured as the proportion of Rhizophora racemosa alleles (i.e., 0 indicates individuals with alleles with R. mangle ancestry, and 1 denotes R. racemosa ancestry). As admixed individuals (gray circles), we considered those individuals from the sympatry zone presenting Q-values lower than 0.80 for the typical R. mangle or R. racemosa inferred clusters in the STRUCTURE K = 2 results and those that were morphologically identified as R. × harrisonii in the field F I G U R E 4 Genetic structure patterns of Rhizophora racemosa and R. × harrisonii. Top: (a) Scatterplot of the discriminant analysis of principal components (DAPC), with the first two principal components represented with red circles and light orange squares representing R. racemosa and R. × harrisonii, respectively. Bottom: Bar plots of STRUCTURE 2.3.4 Bayesian population assignments, with each bar representing a single individual, and each color referring to one inferred group, for two groups (K = 2) Given these partially contrasting results, we performed hierarchical AMOVA to rank models of grouping that could explain the division of independently evolving metapopulation lineages. We grouped individuals based on one a priori hypothesis (R. mangle, R. racemosa, and R. × harrisonii comprise three different groups) and three a posteriori hypotheses (Table 3). When we organized individuals according to their morphological identification carried out in the field (model A), we observed that this grouping explained approximately 66.5% of the total genetic variation. Model B tested the grouping based on DAPC, NEWHYBRIDS, and STRUCTURE K = 4 scenario results (i.e., R. racemosa and R. × harrisonii compose a group separate from R. mangle). This model yielded slightly better performance than model A, explaining 67.5% of the total genetic variability.
We created models C and D based on model B to contemplate the expected and observed geographic structure of R. mangle. We con- America was more genetically similar to samples from the southern populations (see below). Models C and D explained 61.8% and 65.6% of total genetic variation, respectively. Thus, according to our criteria, despite the slight differences among models, the hierarchical AMOVA analysis supports the grouping described by two groups: One composed of individuals identified as R. mangle and one composed of individuals with R. racemosa and R. × harrisonii morphological traits (Table 3).

| Intragroup genetic structure
Due to the genetic structure pattern we observed in the three morphologically different Rhizophora species from the Neotropics, we conducted intragroup rather than intraspecific analyses. At the intragroup level, there was no evidence of systematic stuttering, null alleles, large allele dropout or LD across samples. For R. racemosa and

R. × harrisonii, STRUCTURE and DAPC clustered individuals in groups
that largely corresponded to their corresponding morphotypes. The Bayesian clustering method indicated that the K = 2 scenario better summarized the information than the other K scenario ( Figure   S2), with considerable admixture among inferred groups. Although there was some correspondence between species and assigned group, R. × harrisonii individuals were more genetically homogeneous, whereas specimens identified as R. racemosa presented a wider variation of genetic composition, with some presenting the typical R. × harrisonii genotype (Figure 4). This admixture pattern was clearer in the DAPC results. According to the multivariate clustering method, a scenario with K = 4 is a reasonable representation of our data. The inferred population structure showed one cluster that was exclusively composed of R. racemosa individuals and a second cluster of R. × harrisonii plants. The two remaining inferred groups overlapped and both comprise of both morphological species (Figure 4).
The genetic diversities of the R. racemosa and R. × harrisonii populations did not present substantial variation according to the genetic diversity indices (Table 4). In contrast, for R. mangle, we observed large differences between samples from the northern and southern coasts of Brazil, similar to results reported previously for this species (Pil et al., 2011). Samples from the northern populations exhibited higher genetic diversity summary statistics and a higher frequency of private alleles (Table 4). However, the genetic structure presented some incongruities with previous studies Pil et al., 2011;Takayama et al., 2013).
The BIC scores from the DAPC analyses showed that K = 3 was a reasonable group number for our data. The genetic structure pattern yielded by the multivariate clustering method showed that individuals from the northern populations are more variable than are those from localities south of the NEESA (Figure 5), supporting the genetic diversity indices (Table 4). However, DAPC also clustered together individuals from both regions, indicating that there is no abrupt break between them, as supported by STRUCTURE results. Estimates of lnL and the ad hoc statistic ΔΚ presented somewhat discordant but complementary results ( Figure S2). ΔΚ suggested two groupings that supported our data, K = 2 and K = 4, whereas lnL smoothly increased from K = 4 to K = 7 and then decreased from K = 8 to K = 10 ( Figure   S2). As a conservative estimate of population structure, we chose four clusters. Considering both scenarios, K = 2 and K = 4, there is a pattern of genetic structure in which populations south of the NEESA are more homogeneous than those from the northern coast ( Figure 5).
Additionally, the STRUCTURE results showed that within the northern coast, the composition of "southern" gene pools was reduced westward of NEESA.
This westward differentiation was also observed when we explicitly considered the spatial locations of our samples for sPCA. The sPCA results indicated no local structure (p = .522) but some global structure (p < .001); furthermore, when we plotted the sPCA first positive eigenvalues on a map, the pattern supported the westward gradient, we observed in the northern populations based on the STRUCTURE results ( Figure 5). There was an abrupt decrease between the first and second positive eigenvalues, which we interpreted as evidence of strong spatial structure and therefore a reasonable representation of the information we obtained.

| DISCUSSION
Our independent sampling and use of a new set of molecular markers enabled us to build on previous studies on the genetic diversity of Rhizophora mangle along the Atlantic coast of South America (Pil et al., 2011;Takayama et al., 2013). The results indicated a genetic subdivision between the populations sampled north and south of the NEESA. This pattern has also been observed in two species of Avicennia , which is a true mangrove genus (Tomlinson, 1986), and H. pernambucensis (Takayama et al., 2008), a mangrove associate (Tomlinson, 1986). Our results are consistent with these previous findings, and the present study complements and extends them by analyzing previously unsampled populations along the northern coast and individuals identified as R. racemosa and R. harrisonii using a newly developed set of microsatellite markers.

| Rhizophora mangle genetic structure
As expected based on previous findings, we observed that R. mangle belongs to two major populations. The subdivision that we observed did not show complete nonoverlap between the populations north and south of the NESSA, with a few admixed individuals found in the northern populations, as observed for A. germinans or, more subtly, for A. schaueriana . This trend is evident in the map of the sPCA first eigenvalue scores, which gradually increased westward from the southern population to the northwestern population ( Figure 5). Considering the STRUCTURE results, this tendency was also clear for all evolutionary scenarios, that is, considering K = 2, K = 3, and K = 4 ( Figure 5). Recovering this difference between Avicennia  and Rhizophora, genetic structure pattern was possible only because our sampling scheme included two additional localities within the north coast compared with previous studies (Pil et al., 2011;Takayama et al., 2013). Regardless of the method used, genetic structure consisted of an admixture gradient between two gene pools within the northern coast, whereas populations from the southern coast were homogeneous with low genetic diversity ( Figure 5 and Table 4).
The major north-south subdivision observed in the genetic variation of R. mangle is likely maintained by superficial ocean currents, as previously suggested for this species (Pil et al., 2011) as well as A. germinans, A. schaueriana, , and H. pernambucensis (Takayama et al., 2008)  and R. mangle ( Figure 5). Rhizophora propagules are larger, live longer, survive for at least 1 year in fresh or salt water and require a longer time to become established than Avicennia propagules (Rabinowitz, 1978 The findings presented herein, coupled with previously published research, suggest that oceanic current bifurcation is a key environmental driver that maintains and shapes the genetic diversity of R. mangle between the southern and northern populations and within the northern group, corroborating previous studies Pil et al., 2011;Takayama et al., 2008Takayama et al., , 2013. The processes that originated this pattern of genetic structure, however, are likely to be much more intricate. One possible explanation for the genetic structure and the lower genetic diversity of the southern populations is that northern populations may represent descendants of older populations (refugia) that expanded southwards following postglacial climate changes (Pil et al., 2011). Although our R. mangle findings do not reject this process, our analyses considering not only this species but also R. racemosa and R. harrisonii suggest that an additional process may be influencing R. mangle genetic structure. Interspecific hybridization and introgression with R. racemosa and R. harrisonii may also be a key biological process that influences genetic variation and structure along the Atlantic coast of South America.

| Rhizophora species complex in Northeastern South America
The genetic differences among Rhizophora mangle, R. racemosa, and R. × harrisonii differences are complex (Cerón-Souza et al., 2010;Takayama et al., 2013) and much more extensive than those observed in AEP Avicennia species (Mori, Zucchi, Sampaio, et al., 2015;Nettel, Dodd, Afzal-Rafii, & Tovilla-Hernández, 2008). Although there are morphological (Cerón-Souza et al., 2010;Tomlinson, 1986) and physiological (Cerón-Souza et al., 2014) differences between these taxa, species boundaries are permeable (Cerón-Souza et al., 2010Takayama et al., 2013). Our results showed that individuals identified in the field as R. mangle comprise a relatively well-defined genetic cluster that is different from the cluster composed of individuals with R. racemosa and R. × harrisonii morphological traits.   Figure 2). Additionally, the hierarchical AMOVA results supported this grouping (model B) as the best representation of our data among the a priori and a posteriori hypotheses that we tested. However, this model was only slightly better than model A, which separated R. mangle, R. racemosa, and R. × harrisonii, and model D, which accounted for the previously discussed R. mangle geographic subdivision (Table 3).
At a finer scale, regarding the genetic structure within the group composed of R. racemosa and R. × harrisonii, there are two main clusters that correspond to the morphologically identified individuals and two intermediate clusters (Figure 4). This clustering suggests that the traits that are traditionally used to differentiate these morphotypes might be genetically and evolutionarily meaningful. However, due to the porous species boundaries observed in this species complex based on our genetic data and on previous genetic studies (Cerón-Souza et al., 2010;Takayama et al., 2013), the systematic relationships between these species should be considered with caution.
Under the "unified species concept," which allows for decoupling of species conceptualization and species delimitation (De Queiroz, 2007), it is difficult to distinguish Rhizophora species from the AEP biogeographic region because of their nonreciprocal distinctiveness based on morphological and genetic information. As interspecific hybridization and introgression may interfere with species integrity, and because these processes are pervasive in the studied species complex, we advocate for caution when using these morphological traits for species identification due to the extensive admixture between the R. racemosa and R. × harrisonii groups and considering the entire species complex. Therefore, a synthetic reevaluation of current AEPregion Rhizophora taxonomy would be beneficial, preferably simultaneously considering morphology and genetic data and additional lines of evidence to corroborate or reject the current classification.
The incongruence between species boundaries and morphological traits and its associated taxonomic uncertainty lead to conservation and management issues. In groups without clear delimitations, as it is the case for the AEP Rhizophora species complex, it is intrinsically difficult to identify and record individuals using an objective strategy.
Thus, we recommend that different scenarios of classification, namely individual morphological "species," hybrids, groups of species and the entire species complex, should be used as units for applied purposes (Fitzpatrick, Ryan, Johnson, Corush, & Carter, 2015 Acknowledging that introgression and hybridization are relevant biological processes in this genus would be a useful feature for consideration by conservationists and managers for efficient long-term conservation plans. Using each grouping may present desirable and unsatisfactory consequences depending on the geographic and temporal scales considered, and thus a case-by-case approach is advisable based on conservationists and managers' objectives.
Our results highlight that hybridization and introgression are major biological processes shaping the genetic diversity of the Rhizophora species complex in the Atlantic basin of South America. When we considered the K = 2 scenario of the STRUCTURE results, admixture was evident between R. mangle and the other two species, as was more clearly shown by the INTROGRESS results ( Figure 3). This finding suggests that there might be asymmetric introgression toward R. mangle, supporting the adaptive evolution hypothesis of this species complex based on physiological traits (Cerón-Souza et al., 2014). In the Pacific basin, R. mangle and F 1 hybrids presented higher salinity tolerance than did R. racemosa, suggesting that introgression and hybridization could be maintained by adaptive evolution (Cerón-Souza et al., 2014). In the Atlantic basin, R. racemosa and R. × harrisonii are more geographically restricted than R. mangle, such that these species are not recorded south to the NEESA, and they are not homogenously distributed in the northern coast of Brazil, where they are found only in estuaries with substantial freshwater input (Menezes et al., 2008). Therefore, our results support the hypothesis that salinity tolerance might be relevant in maintaining this species complex (Cerón-Souza et al., 2010. In addition to blurring species boundaries, a consequence of ancient and ongoing hybridization and introgression could be an increased genetic diversity in populations where the parental species coexist. For example, in East Pacific Central America, ancient interspecific introgression associated with intraspecific population divergence following secondary contact likely increased the genetic diversity of A. germinans (Nettel et al., 2008). Likewise, in the Rhizophora species complex studied herein, interspecific gene flow, or more specifically, asymmetrical hybridization and introgression toward R. mangle, could increase the genetic diversity of the R. mangle northern populations.
Indeed, we observed higher levels of genetic diversity in these populations (Table 4), supporting independent previous findings (Pil et al., 2011). Thus, although we cannot reject a postglacial southward expansion of mangrove forest distributions in the Atlantic basin of South America as an important evolutionary process shaping the genetic diversity of R. mangle (Pil et al., 2011), our results suggest that the evolutionary scenario is more complicated. Our findings highlight that current and historical hybridization and introgression may have increased the level of genetic diversity in northern R. mangle populations and altered the distribution of genetic variation between and within populations, possibly creating a bias in R. mangle demographic analyses that do not account for these evolutionary processes.
In this study, we recovered and refined Rhizophora mangle northsouth genetic structure along the northeastern coast of South America.
This subdivision is shared among R. mangle (Pil et al., 2011;Takayama et al., 2013), A. germinans, A. schaueriana , and H. pernambucensis (Takayama et al., 2008), suggesting that similar extrinsic evolutionary and ecological processes influence how these plants are distributed. Our work highlights the importance of oceanic currents in the dispersal of sea-dispersed organisms such as R. mangle. Finally, we found that past hybridization and introgression processes play important roles in the Rhizophora species complex along the northeastern coast of South America. This interspecies gene flow is likely asymmetric, suggesting that adaptation may play a role in maintaining this hybrid zone. We also thank Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for the research fellowship granted to APS and the research award granted to GMM (448286/2014-9).

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTION
Conceived and designed the experiments: PMF, GMM,and APS.
Performed the experiments: PMF and FMA. Analyzed the data: GMM.
Contributed reagents/materials/analysis tools: APS. Wrote the article: GMM. All authors revised, edited and approved the final version of the article.