Genetic Consequences of Forest Fragmentation in a Widespread Forest Bat ( Natalus mexicanus , Chiroptera: Natalidae)

: Recent historical and anthropogenic changes in the landscape causing habitat fragmentation can disrupt the connectivity of wild populations and pose a threat to the genetic diversity of multiple species. This study investigated the effect of habitat fragmentation on the structure and genetic diversity of the Mexican greater funnel-eared bat ( Natalus mexicanus ) throughout its distribution range in Mexico, whose natural habitat has decreased dramatically in recent years. Genetic structure and diversity were measured using the HVII hypervariable domain of the mitochondrial control region and ten nuclear microsatellite loci, to analyze historical and contemporary information, respectively. The mitochondrial and nuclear results pointed to a differential genetic structuring, derived mainly from philopatry in females. Our results also showed that genetic diversity was historically high and currently moderate; additionally, the contemporary gene ﬂow between the groups observed was null. These ﬁndings conﬁrm that the effects of habitat fragmentation have started to be expressed in populations and that forest loss is already building barriers to contemporary gene ﬂow. The concern is that gene ﬂow is a process essential to ensure that the genetic diversity of N. mexicanus populations (and probably of many other forest species) distributed in Mexico is preserved or increased in the long term by maintaining forest connectivity between locations.


Introduction
The fragmentation of natural habitats is a key issue for biodiversity and poses a threat to the genetic diversity of multiple species [1][2][3][4]. Fragmentation is a process of the change in the spatial structure from a relatively homogeneous environment to one with a progressively less homogeneous structure that is ultimately transformed into a heterogeneous habitat. This can reduce the total area of a given habitat type, splitting the remaining habitat, and even increasing the isolation of remnants [5][6][7]. Habitat fragmentation disrupts the connectivity among populations of various taxa, reducing population genetic diversity and increasing population structuring [8][9][10][11], due to the genetic drift associated with low gene flow [12].
Bats are among the most abundant and diverse groups of mammals in tropical forests, playing a central role in pollination, regulation of insect populations, and seed dispersal [13,14]. Despite their ability to fly, bats are vulnerable to the loss of genetic variation in response to anthropogenic fragmentation in tropical forests [15][16][17]. molecular markers with different inheritance modes (mtDNA and microsatellites) that will allow the evaluation of the historical and contemporary genetic structure and diversity of natalid bats in Mexico.

Sampling
Between 2004 and 2014, tissue samples were collected from Natalus mexicanus specimens inhabiting 21 locations throughout their geographic range in Mexico (Table S1 and Figure 1). Bats were captured using harp traps and mist nets, and wing membrane biopsies were collected with a 3 mm biopsy punch (Fray Products Corp., Buffalo, NY, USA). Tissue samples were preserved in 70% of ethanol and deposited at −20 • C at the tissue collection of the Laboratorio de Biología y Ecología de Mamíferos de la Universidad Autónoma Metropolitana-Iztapalapa (UAMI), Mexico. The captured bats were released, except for some individuals who were preserved as vouchers and deposited at the Mammal Collection of the UAMI (catalog numbers: RLW300713Nme3-RLW300713Nme10). throughout its distribution range in Mexico, where the natural landscapes have been fragmented as a result of human activities, and to test whether population genetic isolation is occurring due to the lack of dispersion. To this end, we have integrated information from molecular markers with different inheritance modes (mtDNA and microsatellites) that will allow the evaluation of the historical and contemporary genetic structure and diversity of natalid bats in Mexico.

Sampling
Between 2004 and 2014, tissue samples were collected from Natalus mexicanus specimens inhabiting 21 locations throughout their geographic range in Mexico (Table S1 and Figure 1). Bats were captured using harp traps and mist nets, and wing membrane biopsies were collected with a 3 mm biopsy punch (Fray Products Corp., Buffalo, NY, USA). Tissue samples were preserved in 70% of ethanol and deposited at −20 °C at the tissue collection of the Laboratorio de Biología y Ecología de Mamíferos de la Universidad Autónoma Metropolitana-Iztapalapa (UAMI), Mexico. The captured bats were released, except for some individuals who were preserved as vouchers and deposited at the Mammal Collection of the UAMI (catalog numbers: RLW300713Nme3-RLW300713Nme10). Specimen collection protocols and animal handling followed the Institutional ethical guidelines set by the American Society of Mammalogists [52] and the ethical guidelines of the División de Ciencias Biológicas y de la Salud, Universidad Autónoma Metropolitana-Iztapalapa [53] (Project "Biology and Ecology of bats in Mexico," approved by the Consejo Divisional de Ciencias Biológicas y de la Salud. Session 17.18. Date 28 November Specimen collection protocols and animal handling followed the Institutional ethical guidelines set by the American Society of Mammalogists [52]

Mitochondrial DNA
Total DNA was extracted from 245 Natalus mexicanus specimens following the protocol of the WizardSV Genomic DNA Purification System (Promega) kit. For these 245 individuals, a 331-bp fragment of the HVII domain of the mtDNA control region was amplified via polymerase chain reaction (PCR), using the primers L16517 and HSC [54] and following the conditions and the modifications of the primer HSC by [50]. Sequencing was performed with the Big Dye Terminator Kit (Perkin-Elmer, Norwalk, Connecticut) on an ABI 3130xl automatic sequencer (Applied Biosystems, Foster City, California). The sequences were edited and aligned with Geneious v. 5.6.4 [55] using the ClustalW algorithm and were subsequently adjusted visually.

Microsatellite Loci Amplification
Microsatellite loci amplification was made from tissue samples of 171 individuals from 11 localities of Mexico by means of ten dinucleotide microsatellite primers previously developed for the species (Nm1-Nm10- [56]) that were used along with the PCR conditions described above. Fragments were read on an ABI PRISM Genetic Analyzer 3130XLl sequencer, with a LIZ (GeneScan™ 500 ® LIZ Size Standard) as the allele-size standard. Allele size was estimated using GeneMarker v. 2.4.2 (SoftGenetics, LLC, State College, PA, USA).

Genealogical Analysis
Genealogical relationships between haplotypes were determined by a network of haplotypes through the median-joining method, with the software Network v. 4.6.1.3 [57]. Loops were resolved according to the criteria of [58]. Genetic distances between the genealogical groups (haplogroups) obtained were calculated with MEGA v. 5.0.5 [59], using the Tamura-Nei model (TrN). Haplotype diversity (h) and nucleotide diversity (π) were estimated for each locality, as well as for the groups obtained, using DnaSP v. 5 [60].
The genetic structure between and within localities was determined with molecular analysis of variance (AMOVA) with the software Arlequin v. 3.5.1.2 [61], run at two levels: Nongrouped and among the four groups obtained through the haplotype network (see Section 3). The program Barrier v. 2.2 was used to highlight likely geographic areas of genetic discontinuity [62].
The population dynamics was evaluated with an extended Bayesian skyline plot (EBSP) analysis, using mtDNA control region with BEAST v. 1.8.4 [67] on the CIPRES web portal (specialized in phylogeny); the analysis was run twice, each for 30 million generations, using a coalescent Bayesian skyline model, and an uncorrelated lognormal relaxed clock model. The optimal evolutionary model was estimated with jModelTest v. 2.1.6 (Pacific, Gulf of Mexico and Yucatan Peninsula: Tamura Nei 93). A substitution rate ranging from 0.01 to 0.025 substitutions per site per million years (s/s/my) was used following [68]; an Excel graph was produced. All analyses were performed for each genetic haplogroup obtained; the exception was San Sebastián (SS), given the low sample size (see results) and was included in the haplogroup Pacific.

Historical Gene Flow
The relative mutation-scaled migration rates (M) between the four mitochondrial groups obtained in Network and the relative effective population size (θ) were estimated using Markov chain Monte Carlo simulations in Migrate-n v. 3.7.2 [69] under a Bayesian inference model and with a constant mutation rate. A random tree was used as the baseline genealogy. The parameters of the first run were used as baseline values for the subsequent run until a converging result was obtained. The Markov chain length was set as 10,000 steps with 1000-step increments. An adaptive 4-chain heating scheme was set at temperatures of 1.0, 1.5, 3.0, and 1.000. A total of 10,000 trees per chain were discarded.

Microsatellite Data Analysis
The presence and frequency of null alleles were confirmed by locus and locality with MICROCHECKER v. 2.2.3 [70]. To confirm that the presence of null alleles has no effect on the results, we calculated F ST and genetic distance values with and without ENA correction (estimation of null alleles) using the software FREENA [71] and performed a Student's t-test with NCSS v. 11 [72]. The deviation from the Hardy-Weinberg equilibrium (HWE) and the ligation imbalance between pairs of loci were calculated using GENEPOP v. 4.0 [73] and applying the sequential Bonferroni correction to the significance level of p < 0.05 [74].

Population Structure and Genetic Diversity
The genetic structure was evaluated through a Bayesian clustering analysis with STRUCTURE v. 2.2 [75], under conditions of 1,000,000 burn-in and 500,000 Monte Carlo Markov chain, testing clusters from K = 2 to 11, with 20 replicates per K. The most likely number of genetic clusters (K) was determined by estimating Delta K (∆K) and the logarithmic probability of K, In P (K) = L (K) [76] using the Structure Harvester website [77]. On the other hand, the distribution of genetic variation between and within populations was analyzed using molecular analysis of variance (AMOVA), based on F ST and R ST with 30,000 permutations using Arlequin v. 3.5.1.2 [61]. Genetic diversity by locality was obtained using GenAlEx v. 6.3 [78], estimating the number of alleles (Na), exclusive alleles (NP), observed heterozygosis (HO), and expected heterozygosis (HE).

Contemporary Migration Rates
Gene flow among the four groups (see Results) identified by STRUCTURE was estimated in BAYESASS v. 3.0.4 [79]. These programs use different models to estimate gene flow rates. BAYESASS uses an assignment method and does not incorporate genealogy; besides, it reflects the gene flow that occurred only in the past 1-3 generations. The BAYESASS analysis was first run with microsatellite data using the default delta values for allelic frequency, migration rate, and inbreeding. Subsequent analyses incorporated different delta values to ensure that the proposed changes between chains at the end of the run were between 40% and 60% of the total chain length [79]. Once the delta values (∆A = 0.40, ∆m = 0.45, and ∆F = 0.60) were within the accepted proportion (∆A = 0.15, ∆m = 0.15, and ∆F = 0.14) for four genetic groups, analyses were run three additional times (10 million iterations, one million burn-in, and a sampling frequency of 5000) with different random seeds. All parameter estimates converged.

Mitochondrial DNA Data Analysis
The 245 sequences of the Natalus mexicanus mtDNA control region had 331 bp with no tandem replicates. They showed a base composition of T: 23.8%, C: 26.6%, A: 32.5%, and G: 17.1%, with 271 conserved sites and 60 variable sites, 44 of which are parsimoniously informative sites.

Population Structure and Genetic Diversity
When groups were not defined, the AMOVA analysis showed a higher percentage of genetic variation (71%) between localities and a high differentiation (F ST = 0.71; p < 0.05); the differentiation values per group were also high and significant (F CT = 0.446, p < 0.05) ( Table 1). The analysis using pairwise F ST distances in the Barrier software detected three geographic barriers separating localities into three groups (GM, PM, and SS). These barriers are located in the Sierra Madre Oriental and in the central valleys of Oaxaca, separating GM and PM, and PM and SS, respectively.

Demographic Analysis
The mismatch distribution (SSD) and the Harpending's raggedness index (HRI) of the groups PM-SS (SSD = 0.0011, p = 0.705; r = 0.0042, p = 0.809) and PYUC (SSD = 0.0042, p = 0.727; r = 0.0799, p = 0.472) showed a unimodal distribution; for the GM group, the curve was not strictly unimodal (SSD = 0.0095, p = 0.572; r = 0.0186, p = 0.755) (Figure 3), although it was consistent with recent population growth. In all cases, Fu's F tests were negative and significant, while Tajima's D tests were negative but nonsignificant (Table 2).   The extended Bayesian skyline plot analyses indicated that groups PM and GM increased their effective population size from 125,000 years ago, while PYUC remained constant through time (Figure 3).

Microsatellite Data Analysis
There were 198 alleles with ten loci recorded in 171 individuals. Fifty-five unique alleles were found, where the locality Valle de Bravo (VB) showed a higher number of alleles, but with a low frequency (0.031-0.094), while the locality Pe showed a lower number of alleles, but with a higher frequency (0.031-0.188) ( Table 3). Table 3. Genetic diversity statistics with microsatellites for each locality. Sum and average alleles per locus (Na); unique alleles (Np); observed heterozygosis (HO), expected heterozygosis (HE). Abbreviations in each locality are as in Figure 1. Null alleles occurred in most loci of individuals from at least two localities; however, as F ST and genetic distances (with and without ENA correction) did not differ significantly, no loci were excluded. Eight loci (Nme1-Nme4, Nme6, Nme8-NM10) deviated from the Hardy-Weinberg equilibrium in three or more localities. No disequilibrium in the linkage between loci was detected.

Population Structure and Genetic Diversity
The STRUCTURE analysis identified four genetic groups (K = 4), which differed in composition from those obtained with mtDNA: (1) Pe locality (Baja California Peninsula);  Table 4). HO values ranged from 0.404 to 0.681 and HE from 0.585 to 0.787.    Table  5). The highest migration rates were observed from Group 1 (site Pe: Baja California peninsula) to Group 4 (sites Co, TG, Bo, and Cva: Pacific and Gulf of Mexico). The migration rates from the groups were very low (<5%) and not significantly different from zero. All estimates of migration (m) between groups had 95% confidence intervals that approached zero, indicating little to no recent migration between genetic groups (Table 5).

Contemporary Migration Rates
BAYESASS runs yielded low levels of contemporary gene flow among groups ( Table 5). The highest migration rates were observed from Group 1 (site Pe: Baja California peninsula) to Group 4 (sites Co, TG, Bo, and Cva: Pacific and Gulf of Mexico). The migration rates from the groups were very low (<5%) and not significantly different from zero. All estimates of migration (m) between groups had 95% confidence intervals that approached zero, indicating little to no recent migration between genetic groups (Table 5). Table 5. Estimates of contemporary migration rates (95% confidence intervals) based on microsatellite data among groups of Natalus mexicanus.

Discussion
Our results show that contemporary levels of genetic diversity in N. mexicanus are moderate, and gene flow values between groups are either low or nil, in parallel with high values of population genetic differentiation. These data suggest a reduction in effective population size with isolated populations. A study based on microsatellites showed that, similar to N. mexicanus, the papillose woolly bat (Kerivoula papillosa) currently thriving in a fragmented landscape showed parallel reductions in population density and genetic diversity [16]. Small-sized bats, like N. mexicanus and Kerivoula papillosa, with relatively low mobility, may be more severely affected by landscape alterations regardless of a wide geographic distribution [80].
Our findings also show a differential genetic structure for the mitochondrial control region and nuclear microsatellites, suggesting female philopatry (e.g., [81]). Mating behavior and philopatry affect the population structure in bats [82,83]; besides, bat species with limited long-distance flight capacity demonstrate a greater population structuring relative to species with greater mobility [15]. Although poorly documented, N. mexicanus may display sexual segregation, with females remaining in the cave during the gestation and lactation stages, while most males leave the cave at this time [38,84]. They do not seem to show massive migrations, but they migrate locally in search of the most favorable daytime shelters [85].
The patterns detected using mtDNA and microsatellites showed no genetic differentiation between Natalus mexicanus populations living in northern and southern Mexico [48], neither into two reciprocally monophyletic nor deeply divergent groups, as proposed by [86,87]; conversely, this finding is consistent with the observations previously reported by us [49].
The mitochondrial genealogical analyses identified four lineages (GM, PM, PYUC, and SS), consistent with the geographic structure based on the cytochrome b gene [49]. These lineages respond to historical processes and probably evolved due to the effect of barriers restraining dispersal during the Pleistocene, including mountain ranges, depressions, and lowlands in the Isthmus of Tehuantepec, which were partially revealed by the Barrier software (although it failed to detect the PYUC group).
The intraspecific divergence between groups GM and PM was fostered by mountain ranges such as Sierra Madre Oriental and Sierra Madre Occidental, followed by a subsequent expansion, as evidenced by the mismatch analyses; however, a moderate gene flow between the two groups was recorded. The locality TG (Chiapas, group PM) showed haplotypes shared with group GM (H73, H75, and H77), explained probably by an incomplete lineage sorting or retention of ancestral polymorphisms, similar to the pattern shown by other bat species [88][89][90]. The genetic diversity statistics and demographic testing indicated population expansion in groups PM and GM, i.e., these groups experienced an increasing effective population size from 125,000 years ago, similar to reports for other mammal species [91,92], despite the significant climate changes recorded in this period [93,94]. The SS group is located in the central valleys of Oaxaca, a region with lower altitudinal ranges but surrounded by mountains with altitudes above two thousand meters [95], which would explain the isolation of this population from the group PM; however, high levels of PMto-SS historical gene flow were observed. In addition, the highest levels of historical gene flow were recorded from PYUC to SS, which could be a reflection of a historical dispersal route through the Isthmus of Tehuantepec, as has been documented in birds [96,97]. The separation between GM/PM and PYUC probably results from the influence of the Isthmus of Tehuantepec, which has functioned as a geographic barrier for flying organisms such as bats [50,51,98] and birds [99,100]; this hypothesis was supported by the moderate gene flow values obtained here. For the group PYUC, signatures of demographic stability over time were observed, a finding also supported by paleontological information; these observations suggest that the general climate of the region did not change drastically from the end of the Pleistocene to the present [101,102].
Our analyses based on microsatellite data revealed a pattern inconsistent with the distribution of groups based on mtDNA. Discrepancies in population structure derived from markers with different inheritance patterns have been observed in several organisms, including bats [103][104][105], birds [106,107], reptiles [108], and amphibians [109,110], among others.
This work detected a marked contemporary genetic structure with four genetic groups, none of which are consistent with the groups observed based on mtDNA (Figures 2 and 5). Group 1 consists of PE (located on the Baja California Peninsula), which unlike the results based on mtDNA, showed a vicariant process and became separated from the continental genetic signature. Similar results have been observed in different vertebrate species [111,112].
Microsatellite group 2 includes localities of mitochondrial groups GM and PM, while group 4 clusters localities of mitochondrial groups GM, PM, and PYUC. In these localities, geographic barriers (e.g., Sierra Madre Oriental, Sierra Madre Occidental, Isthmus of Tehuantepec) do not seem to hamper connections between the Gulf of Mexico and Pacific slopes, as well as with the Yucatan Peninsula; nonetheless, the gene flow is low (Table 5). Natalus mexicanus thrives in the interior of forests, although also being able to prosper in the remnants of tropical forests by using resources in the coastal corridors that stretch across landscapes [7,40]. Thus, we can assume that this species may be migrating locally through the Balsas depression or the Isthmus of Tehuantepec lowlands, both suggested as biological corridors for other bat species [50,[113][114][115].
Group 3 includes localities VB and SS, located in the trans-Mexican volcanic belt and the Sierra Madre del Sur physiographic provinces, respectively. It is surprising that these two localities are clustered in the same group despite being more than 400 km apart; both are located in conifer and oak forests within two valleys, one in the State of Mexico and the other in Oaxaca. Isolated populations located within valleys have also been recorded for other mammal species [116,117]. In both localities, mitochondrial and nuclear genetic diversity is relatively low, likely due to genetic drift and inbreeding [15] related to isolation.
The genetic structure based on microsatellites appears to match vegetation types of accord INEGI [118] (Table S1), as reported for birds [119,120]. Accordingly, the group 1 locality in the Baja California peninsula is characterized by sarcocaul shrubland; the localities of groups 2 and 4 have a secondary shrub vegetation of low deciduous forest, and group 3 has pine and oak forests. These results suggest that the genetic differentiation of N. mexicanus in Mexico could be related to the great diversity of habitats where it thrives and are in agreement with those of a previous report [48].
In this context, the high deforestation rate in the habitats of this species is cause for concern; the low deciduous forest alone-the main habitat of the species-loses 650 thousand ha annually [33]. Approximately only 27% of the original cover of seasonally dry forest in México remain as intact forest; if the current trends in deforestation continue, the remaining forest will be heavily reduced and degraded in the near future [32].
Although genetic diversity in areas inhabited by the funnel-ear bat is currently moderate, contemporary gene flow is virtually zero among most groups and low between group 1 (Baja California) and group 4 (GM, PM, and PYUC individuals). This may be a consequence of habitat fragmentation, which should be interpreted as a warning signal, given that the loss of genetic variation and flow can reduce the ability of individuals to adapt to a changing environment, resulting in endogamic depression [121], lower reproduction [122,123], and a higher probability of extinction [124,125]. The information obtained for Natalus mexicanus in this study is also alarming because the current status of most of its populations is unknown [40], and during the development of our field work, we have been able to verify that some populations have either declined or completely disappeared due to human disturbances [126].

Conclusions
This work reports the first population genetics analysis of the Mexican greater funneleared bat (Natalus mexicanus) using mitochondrial and nuclear markers, with contrasting results in terms of genetic structure between both molecular markers. This analysis advances our understanding of the underlying evolutionary processes, revealing historical isolation events resulting from geographic barriers, although with some degree of gene flow, as well as an almost null contemporary gene flow and local effects of affinity to the habitat. We predict that this is possibly due to local dispersal by males through biological corridors of great conservation value for the species.
As the populations studied are located in the main habitats in which the species currently thrives and that present low levels of genetic diversity, our results also support the hypothesis that the increasing fragmentation and exploitation of Mexican tropical forests is affecting the levels of current diversity and contemporary genetic flow between populations. Tropical forest remnants are used intensively by many insectivorous bats, so our findings also support the thesis that forest remnants have considerable conservation value probably for many forest species; therefore, their conservation should have a high priority to keep isolation levels low and thus maintain or restore the genetic diversity of many species linked to this particular habitat.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
List of localities, acronyms, haplotypes of mtDNA control region sequences (number of individuals), and GenBank accession numbers of the samples used in this study.