Abstract

Coatis (Procyonidae; Nasua) are considered the only truly social meso-carnivore mammals in Neotropical forests. In Mexico, white-nosed coatis (Nasua narica) are suspected to have undergone population reduction due to habitat loss and fragmentation and led to a lack of genetic adaptability and genetic isolation throughout its range. We examined patterns of genetic diversity and connectivity of five populations of Nasua narica distributed throughout Mexico (n = 60) by sequencing an ≈ 800 bp fragment of the mitochondrial cytochrome-b gene and also by screening 12 microsatellite loci. We found moderate to high levels of genetic variability for both genetic markers. We recorded twenty-two different cytochrome-b haplotypes throughout the 5 sampled areas and found that each of the sampled population of white-nosed coatis in Mexico harbors unique haplotypes and only three haplotypes were shared among two different populations that were closer geographically. All populations had high haplotype diversity (h) (0.968 ± 0.008 (SD)) but lower levels of nucleotide diversity (π) of 0.007 ± 0.001 (SD). All microsatellite loci were polymorphic in all of the populations and the mean number of alleles per locus was 5.033 ± 1.545 (SD) with expected (HE) and observed (HO) heterozygosity values of 0.774 and 0.664, respectively. However, low Wright F statistic values suggest the existence of a reduced heterozygosity (FST = 0.203, FIS = 0.134 and FIT = 0.310). Significant differences between the five populations confirmed isolation by distance, which suggests genetic structure among five subpopulations.

Keywords

Alleles, Cytochrome b, Heterozygosity, Microsatellites, Nasua narica, Procyonids

Share and Cite:

1. Introduction

Processes that generate genetic structure of species are important for understanding their ecology and evolution. In mammalian carnivores this information is crucial for understanding the mechanisms responsible for adaptive divergence and speciation [1] [2] . Genetic structure is defined as the non-random distribution of genetic variation among individuals within populations as the result of behavioral traits (forming herds, flocks or colonies; [3] [4] ), geographical features (spatial distance or topographical barriers; [5] [6] ), or historical factors (colonization, range expansion or isolation; [7] ). Among these factors, geographical features appear to be one of the major forces shaping the genetic structure of several species, topographical barriers as mountain ranges [8] , water bodies [9] , steppes and deserts [10] , and human infrastructure [11] [12] can become determining factors in structuring populations due to their role in restricting gene flow [6] . Those geographical barriers impede the movement of genes among populations, which together with the dispersal capabilities and ability to occupy a variety of habitats define the genetic variation within subpopulations [2] [13] . Moreover, it has been stated that successful conservation of species needs to maintain several populations across its range, in representative ecological and evolutionary conditions [14] , which it will only be achieved through understanding the genetic structure of populations [13] [15] [16] .

Paradoxically, even though mammals are among the most studied organisms on the planet [17] , there are information gaps for some mammalian species, especially in terms of their degree of variation and genetic diversity [18] . Such is the case of many mesocarnivores species, that despite their ecological importance, they have been understudied in many genetic aspects [17] [19] [20] . Nonetheless, there are studies in mesocarnivores that have found genetic structure in a variety of ecological conditions, such as American mink (Neovison vison; [13] ), coyotes (Canis latrans; [11] ), fishers (Martes pennanti; [21] [22] ), Eurasian otters (Lutra lutra; [23] [24] ), Neotropical otters (Lontra longicaudis; [25] ), Tasmanian devils (Sarcophilus laniarius; [15] ), and wolverines (Gulo gulo; [26] ).

Members of the family Procyonidae, to which the white-nosed coati (Nasua narica) belongs, have been thought to be the least studied of the caniform (dog-like) carnivore families [17] . Genetic markers have been used in populations of procyonids for studies of analysis of relatedness [27] - [32] , phylogeography [17] [33] [34] [35] , and even taxonomic boundaries [36] . However, genetic markers, even widely used microsatellites, have not been employed to describe the variability and genetic structure of different N. narica populations. This mesocarnivore is highly gregarious, forming social groups or “bands” of up to 20 - 30 females and juveniles [27] [37] [38] . Males are solitary individuals and disperse once they reach sexual maturity, joining bands briefly only during mating season [27] [28] [37] [38] . Females are highly philopatric and their home ranges are close to and may even include their birth area [27] [39] [40] . They do not migrate and their home ranges can change between seasons and years, for both bands and males [28] [39] [41] [42] .

It is considered that higher levels of philopatry as well as limited dispersal capabilities, lead to more pronounced fine-scale genetic structure [43] [44] . If animals do not disperse far away from their natal site, closely related individuals could be form part of the same bands [44] . Moreover, it is considered that large scale habitat loss and fragmentation are threatening N. narica [45] . Therefore it is suspected from that in Mexico N. narica populations have been severely reduced and they may have even been extirpated in certain areas [45] [46] . All the aforementioned combined with the peculiar orographic characteristics of Mexico, with a wide variation of reliefs and mountain ranges [47] , certainly could lead to genetic isolation of these populations throughout the country.

In this study, we investigated the magnitude and spatial distribution of genetic diversity of N. narica populations in Mexico. Our objective was to explore the patterns of genetic diversity and connectivity of N. narica populations using mitochondrial and nuclear DNA markers. We hypothesized that genetic structure would be found, and that could be related to the presence of the geographical barriers across the country.

2. Materials and Methods

2.1. Study Area

We used samples of N. narica from five study sites in Mexico (Figure 1), where the species is known to be relatively abundant, that were chosen based primarily on their geographic location because they are separated by the main mountain ranges of Mexico (Figure S1) and by insurmountable distances by the migratory patterns of the species: 1) Punta Raza (21˚00'15''N, 105˚17'20''W), Nayarit (n = 11); 2) Reserva de la Biosfera de Chamela-Cuixmala (19˚22'03''N, 104˚56'13''W), Jalisco (n = 10); 3) Parque Nacional El Tepozteco (18˚53'20''N, 99˚02'00''W), Morelos (n = 13); 4) Parque Museo de La Venta (17˚59'21''N, 92˚55'41''W), Tabasco (n = 12); and 5) Puerto Morelos (20˚51'13''N, 86˚53'55''W), Quintana Roo (n = 14). We considered each sampling site as a separate population of the studied species.

2.2. Physical Restraint and Chemical Immobilization

The animals were captured using Tomahawk traps (Tomahawk Live Trap LLC.) measuring 32 × 10 × 12 inches (81.3 × 25.4 × 30.5 cm) and baited with sardines

and/or fruit (banana, pineapple, mango, etc.). We used an average of 25 traps in each study site, which were laid out in variable length transects; the choice of transects was based on the record of signs (footprints and feces) on the dirt roads and trails of the study areas, traps were set in one particular sampling transect at a study site for 4 to 8 days, or until there were no captures for 3 consecutive days. Traps were patrolled twice a day, after de dawn and before sunset. Captured animals were immobilized with the use of Zoletil 100® (Virbac México S.A. de C.V.) administered at 10 mg/kg body weight [48] [49] . Captures were done under capture license FAUT-0251 extended by the Mexican government to the third author. All samples were collected following guidelines approved by the American Society of Mammalogists [50] .

2.3. Sample Collection

Captured animals were sexed, weighed, measured, sampled, marked, and released within 35 min. Tissue samples were obtained through a small cut from the edge of the ear (approximately 10 mm²) and were fixed in ethanol (96%). Most of the sampling took place between spring 2011 and autumn 2012.

2.4. Laboratory Procedures

2.4.1. DNA Extraction, Amplification and Sequencing

Genomic DNA from tissue samples was extracted according to the universal and rapid salt-extraction protocol [51] . Quality of DNA was assessed by electrophoresis on 1 percent agarose gels in combination with molecular weight standards.

A fragment of ≈800 base pairs (bp) of the mitochondrial cytochrome-b (cyt-b) gene was amplified with two universal primers L14724 and H15915 [52] . PCRs were performed in a volume of 25 µL with 0.2 µL Platinum Taq DNA polymerase (Invitrogen^TM), 2.5 µL 10x RXN buffer, 2.5 µL of dNTP’s, 2 µL of each primer 10 mM, 1.5 µL of 50 mM MgCl₂, 5 µL of DNA, and 9.3 µL of bidistilled water. PCR protocol consisted of 3-min denaturation at 94˚C, followed by 30 cycles of 94˚C for 30 s, 55˚C for 30 s, and 72˚C for 30 s; with a final 6-min elongation step at 72˚C. PCR products were confirmed by electrophoresis in 2 percent agarose gels stained with ethidium bromide using a 100 bp ladder and a negative control. All reactions were performed on a Techne thermal cycler TC-3000. Mitochondrial DNA was purified using the Illustra GFX PCR DNA and gel purification kit (GE Healthcare Co.), cycle sequenced using the BigDye Terminator ver.3.1 (Applied Biosystems Inc.), and done using Applied Biosystems 3730xl automated sequencers (Life Technologies, Thermo Fisher Scientific Inc.).

2.4.2. Microsatellite Amplification and Genotyping

The set of 12 microsatellite primers used in this study was designed specifically for N. narica (Nasua 71, Nasua 193, Nasua 504, Nasua 833, Nasua 895, Nasua 931, Nasua 1050, Nasua 1251, Nasua 1757, Nasua 1785, Nasua 1990, and Nasua 2841― [53] ). For these primers, the polymerase chain reactions (PCRs) were performed in a volume of 25 µL with 0.2 µL Platinum Taq DNA polymerase (Invitrogen^TM), 2.5 µL 10x RXN buffer, 2.5 µL of dNTP’s, 2 µL of each primer 10 mM, 2.5 µL of BSA 100x, 1.5 µL of 50 mM MgCl₂, 4 µL of DNA, and 9.8 µL of bidistilled water. PCR protocol beginning with a 5-min denaturation at 94˚C, followed by 30 cycles of: 94˚C for 30 s, 55˚C for 30 s, and 72˚C for 30 s; and a final 6-min elongation step at 72˚C. All reactions were performed on a Techne thermal cycler TC-3000. PCR products were confirmed by electrophoresis in 1 percent agarose gels stained with ethidium bromide using a 100 bp ladder and a negative control. Products were electrophoresed through an ABI 3730 genetic analyzer and fragment size analysis was performed using GENEMAPPER (Applied Biosystems Inc.), each genotype was confirmed by visual inspection of the electropherograms.

2.5. Data Analysis

2.5.1. Phylogeographic Analysis Based on mtDNA

For the 800 bp fragment of the cyt-b gene, sequence alignment was done manually with BIOEDIT 7.1.3 [54] . To determine the appropriate model of sequence evolution, we used the Akaike Information Criterion (AIC) implemented in JMODELTEST 2.1.4 [55] . A tree of maximum likelihood (ML) was reconstructed with PHYML 3.0 [56] , statistical confidence was estimated by bootstrap resampling with 1000 replications, using a heuristic search with NNI (nearest neighbor interchange) branch-swapping. A Bayesian inference (BI) was performed in MRBAYES3.2.1 [57] , 1,000,000 steps of the MCMC algorithm (with trees sampled every 100 generations), and the posterior probabilities were calculated discarding 25 percent of the initial iterations as burn-in, after the stabilization of log-likelihood values. The phylogenetic trees were rooted using the cyt-b haplotype of Nasua nasua as outgroup, since this is the sister species of N. narica that inhabits South America [33] .

2.5.2. Population Genetic Diversity and Structure Analyses Based on mtDNA

Basic statistics of DNA diversity, nucleotide (π) and haplotype (h) diversities were estimated using DNASP 5.10.01 [58] . We also constructed a haplotype network using the median-joining method implemented in the program NETWORK 4.6 [59] . The subdivision and population structure were examined through paired F_ST tests [60] and analysis of molecular variance (AMOVA; [61] ) in ARLEQUIN 3.5.1.2 [62] . Finally, genetic haplotype of mtDNA equilibrium was tested by Tajima’s D [63] with ARLEQUIN.

2.5.3. Population Genetic Diversity and Structure Analyses Based on Microsatellites

For the microsatellites, genetic diversity was measured by the proportion of polymorphic loci (P) and the mean number of alleles per locus (A) by using Genetic Data Analysis (GDA) software 1.1 [64] . We calculated expected heterozygosity (H_E), and observed heterozygosity (H_O) according to Nei (1987; [65] ), as implemented in the program GENEPOP 4.0 [66] . Also, we performed an exact test for Hardy-Weinberg equilibrium (HWE) per population, by using Markov Chain Monte Carlo with 10,000 dememorization steps, and estimated deviations from HWE for each population through F_IS [67] , as well as linkage disequilibrium using the program GENEPOP with the Bonferroni correction [68] . We used MICROCHECKER 2.2.3 [69] , to detect loci containing errors due to stuttering, short allele dominance, and null alleles through the Brookfield (1996; [70] ) method. To assess genetic differentiation of populations, the Wright’s (1951; [71] ) fixation indices were calculated with FSTAT 2.9.3.2. [72] .

AMOVA was performed using ARLEQUIN, with 1000 repetitions and confidence intervals based on 20,000 repetitions, to assess the proportion of variation found between individuals within populations, followed by intra-population variation. The degree of genetic differentiation among the predefined geographical groups was measured by pairwise F_ST and R_ST and tested for statistical significance using the ARLEQUIN permutation test (10,000 permutations per comparison). We computed a statistical certainty of assignment or exclusion for each individual by using the exclusion simulation in GENECLASS2 [73] , according to the partial Bayesian criterion of Rannala and Mountain (1997; [74] ).

Isolation by distance was evaluated using the subroutine ISOLDE in GENEPOP, which computed a regression of R_ST on geographic distance [75] , and applied the Mantel test with 10,000 permutations per comparison [76] . The effective number of migrants (Nm) was calculated according to the private allele method in GENEPOP [77] . To determine the degree of population genetic structure, we used the program STRUCTURE2.2 [78] . This program applies Bayesian analysis to cluster individuals into subpopulations, with no prior information and no known populations (K). An admixture model was employed and correlated allele frequencies were assumed. The mean posterior probability was calculated for each set of ten runs of K and employed to determine an optimal K, using 1,000,000 MCMC iterations and a burn-in period of 500,000 steps. We analyzed the mean posterior probabilities [Ln(P|D)] for each K, and accepted the value that provided the best fit to the data [78] .

Structure analyses were also performed with another Bayesian method that inferred population genetic structure (BAPS 5.4; [79] ). This method utilizes stochastic optimization to infer the posterior mode of the number of populations. The program was employed for cluster individuals using spatial mixture. We ran the software with a predefined maximum of K = 2 - 10, and repeated the runs five times in order to check the consistency of the results. For each run, the program reports the probabilities for different numbers of subpopulations, K ≤ maximum, and we averaged probabilities over the five runs.

3. Results

3.1. Levels of mtDNA Genetic Diversity and Distribution

We obtained an alignment of 800 bp of the cyt-b gene from 50 sequences gathered from specimens of N. narica; 52 variable sites were recorded, of which 42 were parsimonious informative sites. The nucleotide and haplotype diversity were h = 0.968 ± 0.008 and π = 0.007 ± 0.001, respectively (Table 1). We found a total of 22 haplotypes (Figure 2), haplotype V was the most common and found in 4 individuals. Haplotypes H and N were recovered from two localities (Puerto Morelos and La Venta). Haplotype L was found in three individuals from two localities (Punta Raza and Chamela-Cuixmala). Haplotypes A, M, Q and U were only found one time on individuals from Chamela-Cuixmala, La Venta, El Tepozteco and Puerto Morelos, respectively. The localities with the highest number of haplotypes were Puerto Morelos with seven and La Venta with six, these were the areas with the largest number of collected samples for this genetic marker.

^aOutside of Hardy-Weinberg Equlibrium; ^bShare one haplotype between populations; ^cShare two haplotypes between populations.

3.2. Phylogenetic Relationships and Phylogeographic Analyses

The Hasegawa-Kishino-Yano model (HKY; [80] ) provided the best fit for the data set (−lnL = 1489.026). In general, topology estimates from ML and BI were very similar (Figure 3) both produced a phylogenetic reconstruction with high resolution, showing seven clades formed for haplotype clusters according to their membership in the sampled areas. On one end of ML the phylogenetic reconstruction are the haplotypes from populations of the Pacific Coast (Punta Raza and Chamela-Cuixmala), and on the other end of the ML reconstruction are haplotypes from central Mexico (El Tepozteco) and populations of Southeast Mexico (La Venta and Puerto Morelos). The BI showed a very slight difference, haplotypes from the Chamela-Cuixmala were closely related to those from the Southeast, while haplotypes from central Mexico remains below the southern populations. Our haplotype network, like the phylogenetic reconstructions, showed clusters according to their membership in the study areas (Figure 4): Punta Raza (blue), Chamela-Cuixmala (red), El Tepozteco (green), La Venta (yellow), Puerto Morelos (pink), and shared haplotypes between La Venta and Puerto Morelos (aqua) and Chamela-Cuixmala and Punta Raza (orange).

The AMOVA results for cyt-b indicated that most of the genetic variability (88.4%) was found between individuals within populations, while between-pop- ulation variation had a relatively low value (11.6%). Population pairwise F_ST for cyt-b (Table 2) registered the lowest value (0.061; P < 0.05) between Puerto Morelos and La Venta, both populations of the Yucatán Peninsula (ca. 707.1 km apart). The highest value was recorded between Chamela-Cuixmala on the Pacific Coast and El Tepozteco in central Mexico (0.125; P < 0.05) (ca. 620.8 km apart). All comparisons except one showed significant differences (P < 0.05). The value of Tajima’s D was 1833 ± 0912 (SD, P < 0.05) for all individuals.

Abbreviations: PR = Punta Raza, RBCC = Reserva de la Biosfera Chamela Cuixmala, PNET = Parque Nacional El Tepozteco, PMLV = Parque Museo de La Venta, PM = Puerto Morelos.

3.3. Population Genetic Diversity and Structure Analyses Based on Microsatellites

We successfully extracted the DNA of 60 samples from individuals from the five different Mexican locations. We amplified the loci from each of the samples, and identified that all loci were polymorphic in all of the populations (P = 1.000), and the mean number of alleles per locus (A) was 5.033 ± 1.545 (SD) (Table 1). The mean H_E was 0.774 and the mean H_O was 0.664. Only one population, El Tepozteco, was in HWE disequilibrium (X² = 41.151, d.f. = 24, P = 0.016). One pair of microsatellites (Nasua 504, Nasua 1251 correct superindex) were found in linkage disequilibrium (P = 0.029).

Results of AMOVA for microsatellites showed that the highest proportion of variation was found within populations (67.04%), followed by variation between-population (21.52%) (Table 3). Fixation indices supported the presence of genetic structure among populations (F_ST = 0.203), with a relatively high inbreeding coefficient among them (F_IS = 0.134 and F_IT = 0.310) (Table 3). Population pairwise R_ST and F_ST exhibited the lowest value between Punta Raza and Chamela-Cuixmala (0.668 and 0.138 respectively; P < 0.05) (Table 2 and Table 4), both are located in the Pacific Coast and are the closest sampled locations (ca. 182.5 km). The highest value for pairwise R_ST (0.995; P < 0.05) was recorded between Chamela-Cuixmala on the Pacific Coast and Puerto Morelos in Quintana Roo (ca. 1897.6 kma part). The highest value for pairwise F_ST was recorded between Punta Raza and El Tepozteco (0.258; P < 0.05), populations separated by 686.7 km. All the comparisons showed significant differences (P < 0.05).

The assignment-exclusion test performed in GENECLASS confirmed the membership of each individual in its respective population, with a single exception: the sample RBCC-004 was assigned to two populations: Punta Raza and Chamela-Cuixmala (10,000 replicates, CI 95%). The results of the Mantel test suggested the presence of isolation by distance between populations of N. narica in Mexico (r² = 0.232, y = 0.023 + 0.038x, P < 0.05). Also, the Nm had a value of 0.248 dispersing individuals per generation (P = 0.199), which denotes insufficient values of gene flow between populations.

Abbreviations: PR = Punta Raza, RBCC = Reserva de la Biosfera Chamela Cuixmala, PNET = Parque Nacional El Tepozteco, PMLV = Parque Museo de La Venta, PM = Puerto Morelos.

The STRUCTURE program suggested that sampled N. narica populations most likely represented defined groups. Five population clusters [mean LnP(D) = −6405.1] were identified without any prior population information, these clusters corresponded to the geographical groups (Figure 5(a)). BAPS program also yielded results suggesting 5 distinct clusters corresponding with the geographical groups (Figure 5(b)). The highest posterior probability was observed for K = 5 [mean log (marginal likelihood) = −3226.70], but in this case does not denote the sharing of alleles between populations of Chamela-Cuixmala and Punta Raza.

4. Discussion and Conclusions

This is the first phylogeographic and population genetics study of N. narica in Mexico. To date, little is known about the patterns of genetic diversity of this species. Previous genetic research efforts have been focused on identifying kinship among individuals [27] [28] , or included in phylogenetic studies aimed at understanding relationships of different procyonid species [17] [33] .

The haplotype and nucleotide diversity (h = 0.968, π = 0.007) had similar values to those recorded in ring-tailed coatis (N. nasua; h = 0.783, π = 0.017; [35] ), crab-eating raccoons (Procyon cancrivorus; h = 0.832, π = 0.002; [35] ), and individuals of N. narica from the Yucatán Peninsula (h = 0.888, π = 0.009; [81] ). These values of genetic diversity are very similar to ours, and our findings suggest that surveyed populations of N. narica share a similar pattern with other populations of procyonids (i.e. N. narica, N. nasua and P. cancrivorus).

Furthermore, high haplotype diversity coupled with low nucleotide diversity, has been linked with demographic expansions followed by a decline in population size [82] . This inference is also supported by the results obtained by demographic test, Tajima’s D suggests that populations of N. narica coatis in Mexico previously went through a period of population expansion [83] , but now are experience population decline as some authors have proposed [45] .

The microsatellite results showed that expected heterozygosity was generally higher than observed heterozygosity, and the overall genetic variation of N. narica in Mexico ranged from moderate to high (H_E = 0.751 − 0.825, H_O = 0.646 − 0.693; P > 0.05). Studies based on microsatellite data in other mesocarnivores found similar levels of heterozygosity: 0.673 in American mink [13] , 0.675 in crab-eating raccoons [35] , 0.600 in Eurasian otters [23] , 0.623 in fishers [21] , and 0.623 in ring-tailed coatis [35] .

Although we observed moderate to high genetic variation in populations of N. narica in Mexico, our results indicate that a portion of this diversity has been locally lost and is now spatially subdivided. The AMOVA results showed that within populations a considerable degree of variation remains. However, fixation indices indicate a highly significant differentiation between the collected sites, because apparently bands are functioning as inbreeding entities. Additionally, the low number of migrants per generation found in our results (Nm = 0.248), supports the idea of the low vagility/high philopatric behavior of the species. Previous studies have concluded that bands of individuals of the genus Nasua are “extended” families with some unrelated individuals [27] [39] [84] . Although there have been cases where male individuals have moved more than 20 km between years [85] our results suggest that at least in Mexico, gene flow among the sampled populations is rare and therefore the exchange of alleles among these populations is limited (distances between populations comprises the range of 182.5 km - 1916.2 km).

Several mechanisms that prevent inbreeding in the species have been reported: 1) males avoid mating with related females, 2) absence of temporal continuity of a dominant male in the bands, 3) persistence of promiscuity despite “harem” behavior, and 4) variability in mate choice by females [27] [28] [39] . In addition, female migration occasionally occurs, despite the aggressions from the other members of the band due to the absence of kinship [27] [31] [39] . Still, our data suggests that the mating system and overlapping of band and home ranges of males [28] [39] are probably the cause for the observed values of inbreeding coefficients (F_IS = 0.143 and F_IT = 0.318). Although the value of inbreeding at individuals within populations level was relatively low at individuals relative to the total populations was relatively high [86] [87] . This may be an effect of males being the major contributors to gene flow [37] [39] [40] [41] .

Our results clearly indicated that white-nosed coatis in Mexico are currently not a panmictic population. The Mantel test result showed the existence of genetic isolation by distance; this was reinforced by the population pairwise comparisons (F_ST and R_ST) for both genetic markers. Positive results of isolation by distance have been recorded in other mesocarnivores such as American minks [13] , Eurasian otters [23] , fishers [21] , wolverines [26] , and ring-tailed coatis [35] . The ring-tailed coati provides an interesting comparison. It has been documented that they have marked genetic structure [35] , even though the social groups in South America typically include one adult male throughout the year [88] . This has the potential to limit the genetic variation between litters, and finally between populations. Considering this, it is not surprising that Mexican populations have shown marked structure because of their unique social system.

BAPS and STRUCTURE analyses indicated that the five geographical groups are differentiated and form distinct genetic clusters. Results of the assignment- exclusion test from GENECLASS were consistent with the results of STRUCTURE and BAPS, supporting our hypothesis that the genetic population structure reflects the geographical barriers between populations across the country. Moreover, the number and distribution of haplotypes in each population (i.e., the presence of specific haplotypes and the sharing of some of these only between geographically close populations) is consistent with what was found for both types of molecular markers notwithstanding they have different types of inheritance and temporal resolution [89] [90] . This points to the findings that the distribution of genetic variation found in the Mexican populations of N. narica mirrors geographic barriers that have been present for hundreds of thousands of years, such as the mountain ranges of the Sierra Madre Occidental (SMOc), Sierra Madre Oriental (SMOr), and Trans-Mexican Volcanic Belt (TMVB) which originated in the Cretaceous-Cenozoic, late-Eocene and mid-Miocene respectively [6] [91] [92] (Figure S1). The SMOc runs parallel to the Pacific Coast, from the border with the United States until the TMVB and is bounded on the west by the Gulf of California and on the east by the Central Highlands [92] , virtually isolating Punta Raza and Chamela-Cuixmala. On the other hand, the SMOr runs from the TMVB and it is projected continuously northwest, reaching into the center of Coahuila and Chihuahua, up to the northern border of Mexico [91] , forming a barrier between the populations of the southeast and the rest of the country. Finally, the TMVB extends east-west from the Pacific coast in San Blas, Nayarit and Banderas Bay, Jalisco to the Gulf of Mexico coast at Palma Sola, Veracruz [6] , isolating El Tepozteco from the other ones.

The genetic differentiation between the N. narica populations recorded in this study is probably a result of limited gene flow, due to the combination of the demographic decline [45] , dispersal capabilities of the species [40] , and the isolating effects of the geographical barriers [6] [47] . That pattern has been observed for other carnivore species [21] [26] .

The Maximum Likelihood and Bayesian phylogenies revealed similar patterns, with minor changes between the position of clades and support nodes. According to the distribution of haplotypes in the phylograms, the cluster from Punta Raza was the most divergent, followed by Chamela-Cuixmala, and finally the populations of El Tepozteco and Yucatán Peninsula, La Venta and Puerto Morelos. Separations and divergences observed in the phylogenetic reconstructions are underpinned by the records of the origin and dispersion of this family of procyonids [33] [93] [94] [95] . Fossil evidence suggests that these divergences occurred in tropical and subtropical forests of South and Central America while their dispersal occurred from north to south in the Americas [33] [94] and consequently in Mexico, which is supported by the order in which the groups diverge in the phylogenetic trees.

Like other mesocarnivores, N. narica is often viewed as a resource to be harvested―as vermin or agricultural pests whose impacts on human activities are to be mitigated [20] . Nonetheless, the mesocarnivores play an important role in the ecosystem. Given their smaller size and ability to thrive in diverse habitats, they are usually more abundant than large carnivores, yet their impact within communities is generally assumed to be relatively minor. Mesocarnivores can play important roles as drivers of ecosystem function, structure, or dynamics, by dispersing seeds, controlling prey population and inhibiting competitors [20] . Because of their diet, composed mainly of fruits and arthropods [38] [96] , N. narica is an important seed disperser [97] and potentially a controller of arthropods populations. Furthermore, the N. narica is an important component of the diet of apex predators, such pumas and jaguars [98] [99] . For the proper conservation management of carnivores, no matter their size, and their associated potential ecosystem benefits to be realized, it is necessary to accept that ecological function depends on the genetic integrity and social structure of populations [19] . Finally, a better understanding of carnivore interactions and their functional roles within a whole ecosystem context is crucial before wildlife management actions are applied, to avoid unforeseen deleterious effects [19] .

Based on our results, we can conclude that populations of white-nosed coatis retain a significant level of genetic variability and there is strong evidence of genetic structure between populations distributed along various habitats and separated by well-established biogeographic barriers. Our research provides a first description of the variability and genetic structure of N. narica in Mexico, information that can be an initial basis for setting conservation actions for this species aiming to maintain its genetic diversity and can provide a general framework for this purpose in regard to mesocarnivore species conservation in Mexico. Additional samples from the northern part of this species’ range may help to further elucidate the degree and history of separation between the different populations. Future research should focus on obtaining more samples per population and sampling more and closer populations throughout the entire range of the species, to compare heterozygosity patterns, genetic diversity, and observed population structure. If this could be achieved in addition to the implementation of a monitoring protocol for population genetic variables, it could allow for early detection of population declines. Certainly, population genetics studies can have relevant implications for the implementation of conservation strategies in Mexico for this and other mesocarnivore species.

Acknowledgements

The first author was supported by Consejo Nacional de Ciencia y Tecnología (CONACyT; fellowship No. 257904). This work was supported by Cleveland Metropark Zoo Scott Neotropical Fund (2012) and Proyecto CONACyT Ciencia Básica (156725). We thank the Escuela Nacional de Ciencias Biológicas from Instituto Politécnico Nacional (IPN) and Centro de Investigación en Biodiversidad y Conservación-Universidad Autónoma del Estado de Morelos (CIByC-UAEM) for all facilities. We thank Dr. Gerardo Zúñiga and Dr. Norberto Martínez for their valuable insights into the experimental design and data analysis which greatly improved the quality of this manuscript. We also thank Fernando Montiel and Eduardo Sánchez for their field assistance. The scientific collecting permit was supplied by Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT; FAUT-0251).

Conflicts of Interest

The authors declare no conflicts of interest.

References