Patterns of population structure at microsatellite and mitochondrial DNA markers in the franciscana dolphin (Pontoporia blainvillei)

Abstract The franciscana dolphin, Pontorporia blainvillei, is an endemic cetacean of the Atlantic coast of South America. Its coastal distribution and restricted movement patterns make this species vulnerable to anthropogenic factors, particularly to incidental bycatch. We used mitochondrial DNA control region sequences, 10 microsatellites, and sex data to investigate the population structure of the franciscana dolphin from a previously established management area, which includes the southern edge of its geographic range. F‐statistics and Bayesian cluster analyses revealed the existence of three genetically distinct populations. Based on the microsatellite loci, similar levels of genetic variability were found in the area; 13 private alleles were found in Monte Hermoso, but none in Claromecó. When considering the mitochondrial DNA control region sequences, lower levels of genetic diversity were found in Monte Hermoso, when compared to the other localities. Low levels of gene flow were found between most localities. Additionally, no evidence of isolation by distance nor sex‐biased dispersal was detected in the study area. In view of these results showing that populations from Necochea/Claromecó, Monte Hermoso, and Río Negro were found to be genetically distinct and the available genetic information for the species previously published, Argentina would comprise five distinct populations: Samborombón West/Samborombón South, Cabo San Antonio/Buenos Aires East, Necochea/Claromecó/Buenos Aires Southwest, Monte Hermoso, and Río Negro. In order to ensure the long‐term survival of the franciscana dolphin, management and conservation strategies should be developed considering each of these populations as different management units.


| INTRODUCTION
Along the distribution range of a species, individuals may conform to different populations with different degrees of isolation in response to environmental and/or behavioral factors (Adams & Rosel, 2006;Hoelzel, Dahlheim, & Stern, 1998;Lázaro, Lessa, & Hamilton, 2004;Lessios, 2008;Sellas, Wells, & Rosel, 2005). Commonly, molecular genetics analyses of population structure are used to understand the dynamic of these populations and to facilitate the formulation of effective conservation strategies. Mitochondrial DNA (mtDNA) and microsatellite markers are the most popular markers used for this purpose (e.g., Baker et al., 1998;Costa Urrutia, Abud, Secchi, & Lessa, 2012;Escorza Treviño & Dizon, 2000;Lyrholm, Leimar, Johanneson, & Gyllensten, 1999;Méndez, Rosenbaum, Subramaniam, Yackulic, & Bordino, 2010;Natoli, Birkun, Aguilar, Lopez, & Hoelzel, 2005;Natoli, Peddemors, & Hoelzel, 2008;Pope, Sharp, & Moritz, 1996;Tonione, Johnson, & Routman, 2011). Due to the easy collection, inheritance, lack of recombination, and fast rates of base substitution, mtDNA has been extensively used as a marker in phylogeographic studies (Avise, 1994;Dowling, Moritz, Palmer, & Rieseberg, 1996;Pope et al., 1996;Tonione et al., 2011). However, due to its maternal inheritance, population studies based only on this locus may be biased to female-mediated processes (Avise, 1994;Pope et al., 1996;Zhang & Hweitt, 2003). On the other hand, microsatellites are present widely throughout the euchromatic portion of genomes, are highly polymorphic and apparently neutral, and biparentally inherited. Their introduction in population genetic studies improved our ability to assess genetic diversity, parentage and relatedness, finescale population structure, and recent population history. However, genealogical patterns of relationships cannot be deduced because of the ambiguity of the ancestral information that they contain (Pope et al., 1996;Zhang & Hweitt, 2003). Therefore, incorporating both types of markers may enhance our understanding on the historical and present demographic events that shape the population structure of a species.
The franciscana, Pontoporia blainvillei, is a small dolphin that inhabits coastal waters and estuaries of the South Atlantic coast of America.
Following Gariboldi's et al. (2015) suggestion of increasing the number of loci and samples to acquire a better understanding of the historical and present processes that shape the genetic structure of the franciscana dolphin from the south FMA IV, in this study we combine microsatellite and mtDNA markers to further analyze the franciscana population structure at the southern edge of the species distribution. Argentina. Since three individuals were by-caught simultaneously in this study, and considering that franciscana dolphins would travel in groups comprised in part of close relatives (Costa Urrutia et al., 2012;Méndez, Rosenbaum, Wells, Stamper, & Bordino, 2010;Valsecchi & Zanelatto, 2003), we analyzed only one of those individuals to reduce a potential bias in our analyses. Furthermore, all incidentally entangled or stranded franciscana dolphins used in this study were collected within each locality in sampling events that took place with at least one month of separation between each other. Tissue samples were preserved in 20% dimethyl sulfoxide (DMSO) and/or 96% ethanol.

| Laboratory analyses
Total DNA was extracted from tissue samples using a proteinase K digestion, extraction of proteins with a phenol-chloroform method, and alcohol precipitation of DNA (Sambrook, Fritsch, & Maniatis, 1989).

| Sex determination
In order to determine the sex of individuals, a multiplex polymerase chain reaction (PCR) that targets ZFX and SRY genes was performed (Rosel, 2003). PCR amplifications were performed in a 20 μl reaction volume consisting of 5 μg/ml of template DNA, Buffer 1X (Promega), 0.2 mmol/L of dNTPS, a variable concentration of each primer (Table   S1), 1.5 mmol/L of MgCl2, and 1.25 units of GoTaq polymerase (Promega). PCR cycling profile consisted on an initial denaturation at 94°C for 30 s, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 51°C for 40 s and polymerase extension at 72°C for 40 s, and a final extension at 72°C for 1 min. As positive controls, individuals with known sex confirmed through necropsy were included.

| Microsatellite loci
Ten microsatellite loci developed for other cetacean species (Table S2) and already utilized for franciscana dolphins (Costa Urrutia et al., 2012;Méndez, Rosenbaum, Subramaniam, et al., 2010) were amplified. All  (Table S2) and polymerase extension at 72°C for 40 s, and a final extension at 72°C for 2 min. PCR products were analyzed for length variation using a 3130xl Genetic Analyzer, GeneScan LIZ 600, and Genotyper software (Applied Biosystems, Inc.).
The procedure was performed twice for each sample to confirm the results.

| MtDNA
A fragment of 530 bp from the mtDNA control region was amplified by PCR using primers THR L15926 (Kocher et al., 1989) and TDKD (Kocher, Conroy, McKaye, & Stauffer, 1993). Each PCR, with a final volume of 50 μl, had the same final concentrations of each reagent used for microsatellite PCRs. The PCR cycling profile was as previously used in Gariboldi et al. (2015). PCR products were purified with a commercial kit (AccuPrep PCR Purification Kit, Bioneer) and sequenced using an ABI 337 Automated DNA Prism Sequencer (Applied Biosystems, Inc.). The procedure was performed twice for each sample to confirm the results. for unrelated individuals and for siblings (P ID-SIB ) for the combined set of loci was assessed using GENALEX v6.5 (Peakall & Smouse, 2012).
Since first-order relatives can lead to an overestimation of population structure, to avoid this potential bias in our analyses we used the maximum-likelihood software ML-RELATE (Kalinowski, Wagner, & Taper, 2006) to identify all potential parent-offspring and full sibling pairs. We first determined the most likely relationship (first order, second order, or unrelated) between individuals within each locality and for the whole sampling area. Then, we performed the specific hypothesis test of relatedness with 2 × 10 5 simulations of random genotype pairs to analyze whether the putative relationship fit the data better than the alternative relationship (p = .05): When the putative relationship was parent-offspring, the alternative relationship was full sibling; when the putative relationship was full sibling, the alternative relationship was half sibling. Considering that no parent-offspring relationship was found and none of the putative full sibling pairs first observed (the proportion of putative first-order related individuals found was <1% in all cases) were supported after the hypothesis test, we included the 69 individuals collected to perform the subsequent analyses.
To estimate current migration rate (m) among the four sampled localities (NC, CL, MH, and RN), BAYESASS v1.3 (Wilson & Rannala, 2003) was used. The program utilizes genotypes to estimate rates of recent migration among putative populations. As the method is a nonequilibrium Bayesian method, it does not require HWE within populations (Wilson & Rannala, 2003). Five independent runs of 5 × 10 7 iterations (sampled every 2,000) with 1 × 10 7 iterations discarded as burn-in were performed. Since more reliable results are obtained from runs when the number of proposed changes for allele frequencies (a), inbreeding coefficient (f), and migration rate (m) is between 40 and 60% of the total chain length (Wilson & Rannala, 2003), delta values were adjusted as: m = 0.5, a = 0.6, and f = 0.6.
Five tests for evidence of sex-biased dispersal were conducted using FSTAT 2.9.5 (Goudet, 2001). The first test consisted in comparing population F ST values between sexes; allelic frequencies of the dispersal sex should show higher levels of homogeneity across populations than the philopatric sex. The second test compared relatedness (r) between males and females; within populations, r should be greater in the philopatric sex than in the dispersing one. The third test compared F IS statistics between sexes; a sex-biased dispersal should be reflected in a statistically significant higher F IS for the dispersal sex (Goudet, Perrin, & Waser, 2002). The other tests consisted in calculating the mean and variance of assignment indices (mAIC and vAIC, respectively) to determinate the probability of a genotype originating from the population in which the individual was collected (Favre, Balloux, Goudet, & Perrin, 1997). Immigrants to a population, and therefore the dispersal sex, are expected to have lower mAIC and higher vAIC values than residents (Dubey, Brown, Madsen, & Shine, 2008;Goudet et al., 2002). All tests were performed based on 10,000 permutations.
A Bayesian clustering approach was used to infer putative population differentiation in our data set and to assign individuals to genetic clusters without a prior definition of putative populations, as implemented in STRUCTURE v2.3 (Falush, Stephens, & Pritchard, 2003;Hubisz, Falush, Stephens, & Pritchard, 2009;Pritchard, Stephens, & Donnelly, 2000). Assuming that loci are at HWE and linkage equilibrium within populations, given the number of populations (K), individuals in the data set are probabilistically assigned to one or more populations. Since individuals within a population may have mixed ancestry, we used the admixture model. K ranging from 1 to 7 was evaluated performing 30 independent Markov chain Monte Carlo (MCMC) runs of 5 × 10 6 iterations following a burn-in period of 1 × 10 6 iterations for each value of K. The method of Evanno, Regnaut, and Goudet (2005), which determines the second-order rate of change of the likelihood function with respect to K (∆K), was used to determine most likely value of K over multiple runs, as implemented in STRUCTURE HARVESTER (Earl & vonHoldt, 2012). In order to assign individuals to clusters, a proportion of membership threshold value of q ≥ 0.9 (Costa Urrutia et al., 2012;Trigila, Gómez, Cassini, & Túnez, 2016), which indicates that ≥90% of ancestry can be attributed to a specific cluster, was chosen.
In order to test for isolation by distance (IBD), a Mantel test was performed. IBD v3.23 (Jensen, Bohonak, & Kelly, 2005) was used to examine the correlation between F ST /(1 − F ST ) and the logarithm of the geographic distance between localities. Using a geographic information system (GIS) in ArcGIS software, geographic distances between localities were calculated as the minimum distance by sea between each other. The rejection of the null hypothesis of a flat or negative slope between genetic and geographic distances was used as evidence of IBD.
By using a MCMC approach, migration rates between the sampled locations were estimated with MDIV (Nielsen & Wakeley, 2001). The program estimates the migration rate per gene per generation scaled by the effective population size (M = 2N e m). The finite sites (HKY) model was used. Ten independent runs of 2 × 10 6 iterations each and a burn-in of 5 × 10 5 iterations were performed. For each parameter, likelihood values with the highest posterior probability were accepted as the best estimates.
Similar to microsatellites, a Mantel test based on the mtDNA data set was performed using IBD v2.23 (Jensen et al., 2005).

| Microsatellite analyses
The microsatellite loci studied varied with respect to the number of alleles observed, from 7 to 13. Although putative populations had different sample sizes, the mean number of alleles was similar between them (Table 1). Additionally, AR was similar between localities (NC = 6.44 ± 1.03, CL = 6.34 ± 0.94, MH = 6.70 ± 1.49, and RN = 6.27 ± 1.45). The locus with highest number of alleles was FB17, while FB5 was the one with the lowest number. Alleles found exclusively in one of the putative populations (i.e., private alleles) were not observed in CL. The number of private alleles ranged from 1 in NC to 13 in MH (Table 1). No evidence of null alleles or large allelic dropout was found in our data set.
The analyzed set of microsatellites was powerful enough to discriminate between individuals within putative populations.
Theoretical P ID and P ID-SIB values within populations varied between 5.8 × 10 −12 − 1.9 × 10 −11 and 5.4 × 10 −5 − 7.2 × 10 −5 , respectively. CL, MH, and RN showed similar levels of heterozygosity, and the mean H O ranged from 0.64 in NC to 0.73 in MH and RN (Table 1). We did not find significant deviation from HWE at any of the analyzed microsatellite loci after B-Y correction (all p ˃ .05). Additionally, we did not find significant LD between all pair of microsatellite loci within putative populations after B-Y correction (all p ˃ .05). Therefore, we considered that all loci were genetically independent.
The global test of genetic differentiation was significant (F ST = 0.050; p < 10 −5 ); the greatest source of variation was found within putative populations (95.01%). All pairwise comparisons were statistically significant, except the one between CL and NC (Table 2).
In accordance with the AMOVA results, the Bayesian model-based clustering method showed a clear population subdivision. Considering a range between 1 and 7 populations, the analysis suggested the existence of three populations genetically differentiated (Figure 2). Assuming K = 3, almost all individuals from RN and MH were assigned to their respective cluster (except for two individuals from RN with q = 0.89 and q = 0.82). In the case of CL-NC, twenty-two individuals (55%) were assigned to a common cluster.
Estimated migration rates were low between all locations, except from CL to NC (Table 3). The latter suggests unidirectional migration between both locations. Multiple runs showed consistent results.
The entire tests conducted to detect sex-biased dispersal were not statistically significant (Figure 3). Therefore, our results do not support the hypothesis that females nor males are more dispersive than the other sex.
A positive relation was observed between the geographic and the genetic distance between putative populations (Figure 4a). However, the correlation was not statistically significant (r = 0.325; p = .333).

| mtDNA analyses
From the analysis of the mtDNA control region of the new samples collected for this study, no new haplotypes were found. Instead, three haplotypes previously described in Lázaro et al. (2004), L5 (n = 2), L15 (n = 1), and L22 (n = 1), were found in RN and another haplotype previously described in Gariboldi et al. (2015), G4 (n = 1), was found in MH. When combining our samples with those previously published (Gariboldi et al., 2015;Lázaro et al., 2004), a total of 23 haplotypes were found between NC and RN (Table 4). MH showed the lowest haplotype and nucleotide diversity values, and NC the highest (Table 4).
The AMOVA showed a significant global genetic differentiation between putative populations (F ST = 0.081; p = .003). As for microsatellites, the greatest source of variation was found within putative populations (91.92%). Results of the pairwise comparisons showed no significant differences between NC and CL, NC and RN, or between CL and RN; all other pairwise comparisons were statistically significant (Table 2).
In accordance with the mtDNA population structure results, higher levels of gene flow were observed between NC and CL (m = 32.48) and, to a lesser extent, between CL and RN (m = 11.04). All other population pairs showed values lower than m = 2.28. Similar to microsatellite loci data set, mtDNA data showed a positive but not significant correlation between the geographic and genetic distances between putative populations (r = 0.128; p = .373) (Figure 4b).

| Genetic diversity
Overall, similar levels of genetic variation for both the microsatellite and mtDNA data were observed between localities in this study (Tables 1 and 4), and these were comparable with those previously T A B L E 2 Pairwise genetic differentiation between putative populations for the microsatellite loci and the mtDNA control region found in the FMA IV (Méndez et al., 2008;Méndez, Rosenbaum, Subramaniam, et al., 2010). We found that MH and RN presented several private alleles, in comparison with NC and CL (Table 1), and that these were found in all, but one case, at the edges of the allele sizes (data not shown). This was an expected result since private alleles at a locus tend to be found in the edges of the allele size distribution, that is, commonly they have very long or short repeat lengths with respect to the other alleles at that locus, and this probability is positively correlated with the pairwise F ST estimates between populations (Szpiech & Rosenberg, 2011).
Moreover, it is important to emphasize that among our microsatellite data set, we did not find null alleles that may lead to an overestimation of population genetic differentiation by reducing gene diversity (Chapuis & Estoup, 2007;Putman & Carbone, 2014).
Although this result should be regarded with caution since the ability of the tests performed to detect the bias in dispersal is limited mainly due to the lack of extreme bias in dispersal, the low to moderate dispersal estimates, and the number of loci and samples analyzed (Goudet et al., 2002), previous franciscana dolphin studies also found a lack of sex-biased dispersal (Costa Urrutia et al., 2012;Méndez et al., 2008), supporting our findings.

| Population structure and gene flow
The analyses based on mtDNA and microsatellite loci showed evidence of population genetic structure at the southernmost portion of the franciscanas' range. In general, the estimates of genetic divergence were generally higher for the mtDNA than for the microsatellite data set (Table 2). This was expected based on the lack of sex-biased dispersal ( Figure 3) and the fact that the mtDNA reduces the effective population size to one-fourth that of nuclear genes due to its matrilineal inheritance pattern and haploid nature (Birky, Maruyama, & Fuerst, 1983;DeSalle, Templeton, Mori, Pletscher, & Johnston, 1987) and therefore changes in population allele frequencies may accumulate faster in this marker than in microsatellites.
Nevertheless, only the microsatellite data set revealed a significant differentiation between RN and all other localities (Table 2).
Considering these results, three plausible and nonexclusionary hypotheses can be made to explain the contrasting pattern between molecular markers observed in this study. First, mtDNA and microsatellite loci may be showing different temporal genetic patterns due to differences in their rate and pattern in mutation (Balloux & Lugon Moulin, 2002;Wan, Wu, Fujihara, & Fang, 2004). Second, based on the allele frequency distributions and the ratio R = F ST,organelle /F ST,nuclear , the differences may be reflecting the use of a single organelle marker (mtDNA) versus 10 nuclear markers (microsatellite loci) (Larsson, Charlier, Laikre, & Ryman, 2009). Finally, based on the migration rates estimates (Table 3 and  Among putative populations, migration is a central driving force in evolution that reduces local adaptation (Abdo, Crandall, & Joyce, 2004;Meirmans, 2014). In accordance with a previous mtDNA control region analysis (Gariboldi et al., 2015), high rates of gene flow (Table 3 and mtDNA results) and a lack of genetic differentiation between CL and NC were observed based on the microsatellite loci and mtDNA analyses ( Table 2), suggesting that individuals from both localities form a panmictic population. This may be due to close geographic proximity between localities which may entail similar resource use, as it was previously suggested for this species (Costa Urrutia et al., 2012;Gariboldi et al., 2015;Méndez et al., 2008;Méndez, Rosenbaum, Subramaniam, et al., 2010) and other cetaceans (Fullard et al., 2000;Hoelzel et al., 2007;Natoli et al., 2005). In fact, individuals from NC and CL were found to fed primary on Loligo sanpaulensis .
In the case of MH, in agreement with its previously reported separation as a genetically different population (Gariboldi et al., 2015), both microsatellite and mtDNA data revealed a significant differentiation between MH and all the other analyzed localities (Table 2).
Based on the close relative geographic proximity between MH and CL and the low levels of gene flow between each other (Table 3 and mtDNA results), the observed genetic differentiation between the two localities may be due to different resource specializations.
Prey distribution and abundance are thought to influence the occurrence and distribution of cetaceans (Hastie, Wilson, Wilson, Parsons, & Thompson, 2004) and may lead to intraspecific differentiation by means of resource specialization (Hoelzel, 1998), as it has been previously suggested for other species of dolphins (e.g., Bilgmann, Möller, Harcourt, Gales, & Beheregaray, 2008;Bilgmann, Möller, Harcourt, Gibbs, & Beheregaray, 2007;Möller, Wiszniewski, Allen, & Beheragaray, 2007;Sellas et al., 2005). Since individuals from MH and CL were found to have different diet preferences-individuals from CL feed primary on Loligo sanpaulensis and MH majority preys are Cynoscion guatucupa and Artemesia longinaris -it is feasible to consider that genetic differentiation between these geographically close localities could have been promoted and maintained over time by prey specialization. In fact, Gariboldi et al. (2015) recently proposed that, after its colonization by few maternal lineages, the maintenance in MH of a reduced mtDNA diversity and a relative constant size over time is due to low levels of gene flow between this and other geographically close localities likely promoted by resource specializations. In the case of MH and the other localities, genetic differences may be related to a nonrandom mating between individuals due to the species' relative small home ranges Wells, Bordino, & Douglas, 2013) and the geographic distance between localities, as previously suggested (Gariboldi et al., 2015).
Additionally, although the genetic population structure of many species may be characterized by a pattern of IBD (e.g., Ansmann, Parra, Lanyon, & Seddon, 2012;Hoelzel et al., 2007;Natoli et al., 2005), we did not observe a significant correlation between genetic differentiation and geographic distance among populations (Figure 4).
If equilibrium between the loss of alleles due to genetic drift and their replacement by gene flow between populations exists, the genetic distance between populations will increase with geographic distance due to the changing influence of gene flow and genetic drift as populations become more or less geographically separated (Hutchison & Templeton, 1999). However, this pattern may be affected by time and restrictions in dispersal within an area. Therefore, if time since colonization of a given area is relatively short and gene flow remains strong relative to genetic drift, a pattern reflecting panmixia would persist, as it may be reflecting NC and CL. If lower levels of gene flow exist, a stronger influence of genetic drift through time will occur, and eventually higher genetic differentiation will be observed, as it may be the case of RN. In the case of MH, changes in environmental conditions and resource specializations may conduct to a stronger influence of genetic drift than gene flow, despite the closer geographic distance between this population and others. These cases for deviation from equilibrium gene flow/genetic drift may explain the lack of evidence of IBD in our study area.
Bayesian clustering analyses such as STRUCTURE (Pritchard et al., 2000) can be used to evaluate breaks in allele frequencies, but can overestimate genetic structure in data sets that are characterized by IBD (Frantz, Cellina, Krier, Schley, & Burke, 2009). However, the use of ∆K (Evanno et al., 2005) that may reduce the number of artificial clusters when compared with other Bayesian clustering methods (Frantz et al., 2009), the lack of strong patterns of IBD among our study area ( Figure 2) (Frantz et al., 2009;Meirmans, 2012), and the correspondence with the genetic pattern obtained through the AMOVA and pairwise values of genetic differentiation (  et al. (2004) and included in our study. Therefore, it is possible that NC/CL would not differ from BASW and also that RN would comprise a population genetically distinct from NC/CL/BASW.

| Management and conservation implications
In this study, we have reported a fine-scale genetic structure for the franciscana dolphin over the southernmost portion of the species range, uncovering a new genetic distinct population, RN. Based on these results and previous studies (Gariboldi et al., 2015;Méndez et al., 2008;Méndez, Rosenbaum, Subramaniam, et al., 2010), five populations are found within Argentina: SW/SS, CSA/BAE, MH, NC/CL/BASW, and RN. Furthermore, our study highlights the need to perform multilocus analyses to identify genetically distinct populations since allele frequency distributions, rates of gene flow, mutation rates, and effective population sizes may affect the statistical power of molecular markers (Larsson et al., 2009) and lead to a misinterpretation of the true genetic relationships among populations.
Incidental annual mortality of the franciscana dolphin in Argentina represents up to 2%-5% of its abundance in the area (e.g., Cappozzo et al., 2007;Crespo et al., 2010;Negri et al., 2012) which, according to the International Whaling Commission Scientific Committee (Donovan & Bjørge, 1995), may not be sustainable over time. Additionally, considering our results, MH may have become isolated from geographically close populations due to specializations over resources, whereas RN might have diverged recently from CL-NC due to the geographic distance between them. Consequently, if these populations were severely impacted by certain factors, such as high incidental bycatch, genetic depletion may not be able to be counteracted by gene flow.
Therefore, reformulating Secchi′s et al. (2003) FMAs' division is necessary. The development of conservation and management plans should take into account each genetically distinct population found in Argentina as different management units (sensu Moritz, 1994Moritz, , 2002, considering the mechanisms that may have prompted genetic differentiation between them, as well. However, conservation and management strategies need to be developed upon reliable demographic data. In this regard, the abundance of the species along the Argentina coast was estimated in a single study (Crespo et al., 2010). Additionally, although incidental mortality rates have been assessed previously (e.g., Cappozzo et al., 2007;Crespo et al., 2010;Negri et al., 2012), there are some areas, such as RN, with no information, whereas information needs to be updated in others. Therefore, it is crucial to first carry out demographic studies within each management unit.

ACKNOWLEDGMENTS
We are grateful to all artisanal fishermen, the staff of Fundación