Discordance in maternal and paternal genetic markers in lesser long-nosed bat Leptonycteris yerbabuenae, a migratory bat: recent expansion to the North and male phylopatry

Sampling

Samples were taken from 23 localities in roosting caves and mist nest on field feeding-sites along L. yerbabuenae distribution from 2014 to 2016 (Fig. 2), with sampling permit Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT) SGPA/DGVS/07161/15 (Supplemental File S1) following the Animal Care and Use protocols of the American Society of Mammalogists (Sikes, 2016). Tissue samples were taken with a three mm² biopsy wing punch in an area of the wing with no blood capillaries or nerve terminals. Wing biopsies were fixed in 90% ethanol at environmental temperature and then stored at −20 °C until DNA extraction. Additionally, we used tissue samples from wing punches collected from 2001 to 2003 (Morales-Garza et al., 2007) and tissue samples from wing punches from the tissue collection of the Laboratorio de Ecología y Conservación de Vertebrados Terrestres (LECVT), Instituto de Ecología, UNAM (Table 1).

DNA extraction and amplification

Total genomic DNA was extracted following Paboö’s modified protocol (Gasca Pineda, 2015). Tissue was digested for 12 h at 40 °C in Paboö lysis solution (100 mM NaCl, 100mM Tris HCl and 2mM EDTA, pH8.0) with 20 mg/ml proteinase K, 2% SDS and 0.04M DTT, followed by a phenol: chloroform protocol for DNA isolation (Gasca Pineda, 2015). The quality and amount of extracted DNA was visualized in a 1% agarose gel. We performed gel electrophoresis at 90 V for 30 min; gel was stained with Midori green advance solution and visualized in UV light. We sequenced two mitochondrial DNA regions, cytochrome b (Cyt-b) and control region (D-loop), which are maternally inherited and variable in most mammal species. We also sequenced the region of Dead box associated to the Y chromosome gene (DBY) which is located in the Y chromosome and is paternally inherited, this marker has been used before for phylogeographic studies in phyllostomid bats (Clare et al., 2011).

We obtained 1121 bp of Cyt-b with primers L14125 5′TGAAAAAYCATCGTTGT 3′and H15915 5′TCTTCATTTYWGGTTTACAAGAC 3′(Steppan et al., 1999). For the D-loop region , we amplified 828 bp with primers L15933 5′-CTCTGGTCTTGTAAACCAAAAATG-3′and H637 5′-AGGACCAAACCTTTGTGTTTATG-3′(Oshida et al., 2001). PCR for mitochondrial markers were performed in a final total reaction volume of 15 µl, and contained 2 µl of DNA, 2 U of Taq polymerase (GoTaq Flexi DNA Polymerasa, Promega, USA), 0.4 µM of each primer (10 µM), 1x Taq buffer, 2.5 µM of MgCl₂ (25 µM), 0.2 µM of dNTPs (10 µM) and 7.325µl of H₂O. The PCR profile for Cyt-b was: 5 min of initial denaturation at 95 °C, followed by 35 cycles of 30 s at 96 °C, 1 min at 53 °C, 2 min at 72 °C, and a final extension of 7 min at 72 °C; and for D-loop: 3 min at 95 °C, followed by 35 cycles of 30 s at 94 °C, 45 s at 54 °C, 2 min of 72 °C, and a final extension of 10 min at 72 °C in an ABI Veriti 96-Well Thermal Cycler (Model: 9902; Thermo Fisher Scientific Inc.).

We amplified a 450 bp fragment from the Dead box associated to the Y chromosome gene (DBY) with primers 5′-CCGTTACTTCCATTTTCAAAA-3′and 5′-GCTAAAACCAACGAGATTGGT-3′(Lim, 2007; Lim et al., 2008); the reaction mixture of 15µl total volume contained 2µl of genomic DNA, 2 µM of each primer (10 µM), 200 µM dNTPs (10 µM), 1.5 mM MgCl₂ (25 µM), 2.5 U of Taq DNA polymerase (Promega) and 7.7µl of H₂O. Amplification was carried out as follows: 10 min of an initial denaturation at 94 °C, 36 cycles of 45 s at 94 °C, 30 s at 54 °C, 2:30 min at 72 °C, and a final extension of 5 min at 72 °C in an ABI Veriti 96-Well Thermal Cycler (Model: 9902; Thermo Fisher Scientific Inc.). We sequenced each genetic region with forward and reverse primers at Macrogen USA’s Maryland headquarters (http://www.macrogenusa.com).

Differences in the final sample number among markers resulted from PCR artifacts, Cyt-b was successfully amplified, while for D-loop sequences were low quality for several individuals and were discarded from analyses. We amplified a total of 213 individuals for Cyt-b, 152 for D-loop and 130 for DBY (Table 1; Supplemental File S2). All sequences are available in NCBI GenBank (accession number: Dloop: MT790834–MT790986; Cytb: MT859334–MT859403; DBY: MT913638–MT913767).

Data analysis

Historical analyses

Historical demography analysis

To estimate the demographic dynamics of L. yerbabuenae through time we constructed Bayesian skyline plots with BEAST 1.10.4 (Suchard et al., 2018). We calculated coalescence times for each locus separately considering all individuals (Cyt-b n = 213, D-loop n = 156 and DBY n = 130) using GTR+G+I in a Piecewise-linear nucleotide substitution model (Suchard et al., 2018).

Genealogies and model parameters for each lineage were sampled every 50,000 iterations for 5 × 10⁸ generations under a relaxed lognormal molecular clock with uniformly distributed priors and a pre-burn in of 1000. Demographic plots for each analysis were visualized with Tracer 1.7.1 (Rambaut et al., 2018). To scale the time of the Bayesian coalescence of the Skyline (evolutionary time between real time), we used the last divergence time of the branches in the calibrated tree. We also calculated Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) with a coalescent algorithm in DNAsp 5.10.01 (Librado & Rozas, 2009), as an independent test for demographic inferences and to provide further support to skyline plot results.

Genetic diversity

We analyzed a total of 250 individuals from 23 sites across Mexico (females = 120, males = 130) (Table 1; Fig. 1). We obtained a total of 1,128 bp for Cyt-b (n = 213; 76 haplotypes); 828 bp for D-loop (n = 152; 43 haplotypes) and 450 bp for DBY (n = 130; 70 haplotypes). Genetic diversity (H_d) values were lower for maternally inherited mitochondrial regions (Cyt-b = 0.757; D-loop = 0.8082) than for paternally inherited DBY (0.91; Table 2), a similar pattern was observed for nucleotide diversity (Table 2).

Genetic structure

The haplotype networks (Fig. 3) exhibited differences among them. Mitochondrial Cyt-b and D-loop networks showed a main group with most haplotypes from all localities connected, and smaller groups consisting mainly of members from the Peninsula of Baja California-West of Balsas localities (Figs. 3A–3B). The DBY haplotype network consisted of three-star groups with geographic congruence (West of Balsas, Central Mexico and Isthmus of Tehuantepec) connected by several mutational steps (Fig. 3C).

Pairwise F_ST showed low genetic differentiation between localities for the mitochondrial markers (Tables S1–S2) but higher genetic differentiation for DBY (Table S3). For Cyt-b and D-loop, the results suggest that locations from the West of Balsas region are moderately differentiated from East of Balsas, with values below F_ST = 0.4. In contrast, the paternal lineage (DBY marker) exhibited a well-defined genetic structure between the Isthmus of Tehuantepec-East of Balsas and West of Balsas (Pacific Coast) with values ranging from 0.6 to 0.98.

Bayesian population analysis (BAPS) for Cyt-b (log likelihood = −10512.121) and D-loop (log likelihood = −10990.658) detected two genetic clusters (K = 2; Fig. S1), consistent with East of Balsas vs. West of Balsas (or Pacific Coast) groups. In the case of the DBY gen, BAPS detected three different genetic groups, congruent with the geographic distribution of sampled localities (K = 3; log likelihood = −6397.346; S2. Fig. S1), corresponding with West of Balsas, Central Mexico, and Isthmus of Tehuantepec zones (Fig. 1).

AMOVA analyses revealed that variance was better explained without considering a hierarchical grouping for both mtDNA regions, while for DBY clustering provided a slightly better explanation of genetic variance (Table 3). In all cases, most of the genetic variation was found within-localities (55.2% for Cyt-b, 75.7% for D-loop and 66.6% for DBY) (Table 3).

Historical analyses

Divergence times

The Bayesian genealogy for Cyt-b (Fig. 4) and D-loop (Fig. 5) show that L. yerbabuenae is a monophyletic clade consisting of several lineages. The chronogram of the most probable tree constructed with mitochondrial Cyt-b (Fig. 4) indicated that L. yerbabuenae originated 4.03 mya (95% HDP, 2.27–8.63 mya). In contrast, for mitochondrial D-loop (Fig. 5), we obtained a deeper date of 10.26 mya (95% HDP, 8.73–11.82 mya). The chronogram based on chromosome Y DBY gene (Fig. 6) suggests a similar (but not overlapping, according to the 95% HDPs) date for the origin of L. yerbabuenae, 12.23 mya (95% HDP, 9.99–13.78 mya).

Regarding the geographical distribution of maternal lineages, for Cyt-b, individuals from SJR, DF and SSF are paraphyletic and show long branches denoting a possible bottleneck, and the rest of the individuals form a monophyletic clade with short branches denoting a possible population expansion (Fig. 4). For D-loop, there is no clear geographic structure, but we can also see that some groups have longer branches suggesting that different lineages might have undergone different demographic trajectories (Fig. 5). For the paternal lineages, the genealogy for the DBY gene shows several groups with a few geographical correlations among haplotypes (Fig. 6). The main group contains haplotypes from States of Pacific coast of Mexico and some members from Oaxaca and Morelos (3.06 mya; West Balsas, Pacific Coast and Northern States, Fig. 6). Regarding individuals from Central Mexico, we can see that most individuals from DF form a monophyletic clade with long branches (3.25 mya; Fig. 6), and most individuals from SJR in Oaxaca, form another clade (0.96 mya; Fig. 6). There is a group exclusively composed by members from Jalisco (1.05 mya; Fig. 6).

Historical demography analysis

Bayesian skyline plots (SLP) indicate a late Pleistocene demographic expansion, starting at ∼130,000 ya for Cyt-b (Fig. 7A) and ∼500,000 ya for D-loop (Fig. 7B). In contrast, for the DBY gene (Fig. 7C), SLP basically exhibits demographic stability, with an older and constant but slight demographic expansion, about 6-7 mya, coincident with its own evolutionary history and divergence times dated in the gene genealogy above.

Tajima’s D was negative and significant for both mitochondrial genes (Cyt-b Tajima’s D = −1.73, p = 0.007; D-loop Tajima’s D =  − 1.6, p = 0.023; Table S4), suggesting a possible historical population expansion, while Tajima’s D for the DBY was positive and non-significant (Tajima’s D = 4.49, p = 0.99; Table S4), suggesting demographic stability at deeper times. Fu’s neutrality tests exhibited negative but non-significant values for the mitochondrial regions (Cyt-b = −5.86, p = 0.19; D-loop = −2.05, p = 0.37799 and was positive but also not significant for DBY Fu’s = 0.44, p = 0.87).

Past distribution models

We retained seven bioclimatic layers (Isothermality, Min Temperature of Coldest Month, Temperature Annual Range, Mean Temperature of Coldest Quarter, Precipitation of Wettest Quarter, Precipitation of Warmest Quarter, Precipitation of Coldest Quarter) to build distribution models. The potential distribution models for L. yerbabuenae (Figs. 2; 8) showed stability and good support, for each distribution model the area under the ROC curve (AUC) was >0.8. The projection to the Holocene (Fig. 8B) suggests that the distribution area of species has been stable during this period. Nevertheless, LIG (Fig. 8D) and LGM (Fig. 8C) models do not recover suitable environmental conditions for the presence of the species at higher latitudes, such as Arizona, USA. However, environmental suitability is predicted in West and East Balsas and the Isthmus of Tehuantepec for all periods (Fig. 2).

For the current model, the area with the greatest environmental suitability for the presence of the species is located close to West Balsas (along Sierra Madre Occidental), south of the peninsula of Baja California, the states of Oaxaca and Chiapas on the Isthmus of Tehuantepec zones (coast of the Pacific Ocean states and portions of Central Mexico; Fig. 8A). For the Holocene, L. yerbabuenae had a marked contraction specially in Central America, in contrast with the expansion observed in Central Mexico, West Balsas (and Pacific Coast of Mexico) and south of Baja California Peninsula (Fig. 8B). During the LGM, suitable environmental conditions occurred along the Pacific coast including West-Balsas to East-Balsas comprising areas of Central Mexico towards the southern area of the Chihuahuan Desert (Fig. 8C). The distribution of L. yerbabuenae apparently contracted during the LIG period but exhibited suitable conditions for the presence of species in West Balsas and the Isthmus of Tehuantepec. Finally, during the transition between the Holocene and Current periods there is an evident expansion towards the north, particularly in Sonora and Arizona (Figs. 8A–8B).

Genetic diversity

Levels of genetic diversity as measured by haplotype diversity for the three molecular markers (Table 2) were consistent with previous studies of this species (Wilkinson & Fleming, 1996; Morales-Garza et al., 2007; Arteaga et al., 2018; Menchaca et al., 2020), and of other American bat species (i.e., Tadarida brasiliense (Russell, Medellín & McCracken, 2005), Sturnira parvidens (Hernández-Canchola & León-Paniagua, 2017), and lower than values reported for Artibeus jamaicensis (Ruiz, Vargas-Miranda & Zúñiga, 2013).

The highest haplotype diversity for both mitochondrial markers was found in the Central Mexico near the Isthmus of Tehuantepec. DBY showed the highest genetic diversity in populations from the Pacific Coast (including Baja California Sur = 0.90; Juxtlahuaca, from Guerrero = 0.94 and Arandas, Jalisco = 0.9789) (Table S5). The conflicting patterns of genetic diversity between DBY and mitochondrial markers could be related to different possible causes, including their mode of inheritance, patterns of molecular evolution and distinct evolutionary and demographic history (Tosi, Morales & Melnick, 2000), and to L. yerbabuenae’s female migratory behavior and male philopatry (Herrera, 1997; Rojas-Martínez et al., 1999); there is a sex bias during spring-summer seasons, females are more abundant at northern latitudes when they live in roosting caves and males leave mating caves at the southern area of their distribution (Herrera, 1997; Rojas-Martínez et al., 1999; Medellín et al., 2018).

Genetic structure

Analysis of genetic differentiation, haplotype network, BAPS and AMOVA for mtDNA showed lower genetic structure for female inherited genetic markers than for the male inherited marker, coherent with reports in other bat species (i.e., Tadarida brasiliensis (Russell, Medellín & McCracken, 2005), Miniopterus schreibersii (Bilgin et al., 2008b), Myotis capaccinii (Bilgin et al., 2008a) Myotis nattereri (Rivers, Butlin & Altringham, 2005)). In particular, for the DBY gene, the analyses indicate strong genetic structure with geographical consistency, comprising three groups: West of Balsas, Central Mexico-Balsas, and Isthmus of Tehuantepec (Fig. 1).

Our results support the previous inferences that L. yerbabuenae females move more than males, and that males show philopatric behavior, conducting only local and altitudinal movements (Cockrum, 1991; Herrera, 1997; Rojas-Martínez et al., 1999; Tellez et al., 2000; Medellín et al., 2018). Migratory females apparently promote gene flow among caves minimizing genetic differentiation, while male philopatry and perennial presence promote genetic differentiation in the DYB gene.

Historical analyses

The relationships among haplotypes within L. yerbabuenae clade were different in the three gene genealogies; these differences were marked between the paternally (DBY) and one of the maternally (Cyt-b) inherited molecular markers, although there are some similarities. Genealogical relationship among L. yerbabuenae and sister species is not clear, with D-loop and DBY markers of L. curasoae and L. nivalis forming a monophyletic group. While in the case of the genealogy based on Cyt-b, L. nivalis and L. curasoae form independent branches, and L. nivalis is established as the closest species to L. yerbabuenae clade. In addition, we observe some disagreement in haplotype relationships depicted in the maternally inherited Cyt-b and D-Loop genealogies. In this sense, high homoplasy has been reported for the D-loop in several species (Martes pennanti in Finnilä, Lehtonen & Majamaa, 2001; Homo sapiens in Knaus et al., 2011).

The major difference occurred with calculated divergence times, for DBY marker, divergence times shows 12.01 mya [9.62–13.43 95% HDP] for the origin of L. yerbabuenae. These results reflect the evolutionary history for each gene, in particular differences in mutation rates (Allio et al., 2017), as well differences in the evolutionary and ecological history of each sex. Some authors have stablished that Y-chromosome genes are less polymorphic than mitochondrial genes and this could be a reason for discordant divergences times (Boissinot & Boursot, 1997).

Divergence times of Leptonycteris yerbabuenae for D-loop (10.26 mya [95% HDP, 8.73–11.82 mya]) and the chromosome Y associated DBY marker (12.23 mya [95% HDP, 9.99–13.78 mya]) were older than those obtained for Cyt-b (4.03 mya [95% HDP, 2.27−8.63 mya]). Dates of divergence for DBY and D-loop are consistent with the origin of ecologically related taxa, such as Agave sensu lato (4.6–12.3 mya; Flores-Abreu et al., 2019; and more recently Jiménez-Barrón et al., 2020 reported 9 mya for this group). It is interesting to note that date of divergence for Cyt-b is consistent with the divergence of A. lechuguilla (2.47–6.71 mya; (Scheinvar et al., 2017). This is relevant because Agave and Leptonycteris, are closely associated and they share close evolutionary trajectories (Rocha, Valera & Eguiarte, 2005). Times of divergence are consistent with two pulses of acceleration in the diversification rate of Agave sensu lato, first 8-6 mya, and second 3−2.5 mya (Good-Avila et al., 2006) and with the most recent report with a pulse at 6.18 mya and the second one at 4.91 mya (Jiménez-Barrón et al., 2020). Moreover, these dates of divergence also coincide with a temperature decrease for tropical wet climates in Mexico during glacial periods between 5.3 and 1.8 Mya (Van Devender, 2000). It is worth mentioning that in all genealogies L. yerbabuenae consists of several lineages that may have originated by the isolation of populations during the Pliocene, in particular in the area of Central Mexico. Nevertheless, an analysis of dates of divergence based on genome wide data will allow a better understanding of the processes leading to lineage divergence within Leptonycteris yerbabuenae.

The demographic history of populations also influences the shape of ultrametric trees, where a genetic bottleneck can increase the rate of coalescence of lineages, while a demographic expansion could produce isolated long branches (Ho & Shapiro, 2011; Gattepaille, Jakobsson & Blum, 2013). Accordingly, during the Pleistocene, the species may have undergone geographic expansion that erased geographic structure. In particular, one shows (Cyt-b) signals of demographic expansion, probably associated with Pleistocene climate changes as supported by demographic analyses.

High haplotype diversity and medium-low nucleotide diversity, in addition to a star-like haplotype network, further suggests population expansion (Slatkin & Hudson, 1991; Avise, 2000). SkyLine Plots for mtDNA suggest that a population expansion began approximately 130,000 and 500,000 years ago, during the Pleistocene; Tajima’s D results also support demographic expansion for mtDNA. These results are further supported by past distribution models showing an expansion from the LIG through the Holocene.

Pleistocene climate changes influenced the distribution of a variety of living groups. Overall, highland biota was fragmented during warmer interglacial periods (McDonald, 1993; Metcalfe et al., 2000; León-Paniagua et al., 2007; Ruiz et al., 2010; Bryson et al., 2011; Gugger et al., 2011), followed by expansions related to an increase in temperature (Hundertmark et al., 2002; Hofreiter & Stewart, 2009; De Bruyn, Hoelzel AR & Hofreiter, 2011). The transition between the LIG and Holocene distribution models are consistent with warmer conditions in Mexico 15,000–12,000 years ago, and with the colonization of ecologically related groups such as Cactacea and Agave towards the north of Mexico (Metcalfe et al., 2000). Furthermore, current distribution for L. yerbabuenae agrees with the distribution of several species of the genus Agave (Scheinvar et al., 2017; Scheinvar, 2018).

Female migration in L. yerbabuenae might have occurred following the changing climatic conditions of the Pleistocene. Our models for LIG, Holocene and current periods suggest an expansion to northern areas along the Pacific Coast towards Sonora. This area is home to numerous maternity colonies of L. yerbabuenae, including the important Pinacate roost cave (Medellín et al., 2018). The volcanic record indicates that the most recent volcanic activity in the area occurred 12,000 (+/- 4,000) years ago, while the origin of the current ecosystem has been dated 9,000 years ago (Marshall & Blake, 2009). These dates for geological and biological events coincide with our estimated expansion model for the Holocene. These models are also supported by Arteaga et al. (2018), who reported a local demographic expansion for L. yerbabuenae in northwestern Mexico associated to Pleistocene climate changes.

Demographic analyses among genetic markers are contrasting, suggesting a demographic expansion for the female inherited markers, but not in the male lineages. This recent expansion is further supported by past projection models. Moreover, molecular data and past distribution models suggest that L. yerbabuenae female’s migration to northern maternal roost has been an ecological and dynamic process originated in the Pacific Coast zone of Mexico and Central-South area (Isthmus of Tehuantepec; Figs. 8B–8D). This hypothesis is reinforced because the Pinacate area became ideal for the arrival of females only ∼8000–16,000 years ago, due in part to more favorable climatic and geologic conditions for the bat species during the climatic transition from LIG-Holocene. Moreover, genetic diversity, Bayesian genealogies for the three markers and the model for more stable thermal conditions from the Pleistocene to Current time, provide information to assume that the site of origin of the species could have been located in the Balsas zone of Mexico and Central-South (Isthmus of Tehuantepec; Fig. 1). According to the basal lineage in the gene genealogy, an ancestor of the living populations of L. yerbabuenae existed south of its current distribution. This origin is also supported by the greatest diversity of lineages currently present in the zone of the Isthmus of Tehuantepec. Nevertheless, these biogeographic hypotheses should be tested with other type of molecular markers, such as SNPs, with higher phylogenetic resolution and including the three species of the genus Leptonycteris.