Genetic demography at the leading edge of the distribution of a rabies virus vector

Abstract The common vampire bat, Desmodus rotundus, ranges from South America into northern Mexico in North America. This sanguivorous species of bat feeds primarily on medium to large‐sized mammals and is known to rely on livestock as primary prey. Each year, there are hotspot areas of D. rotundus‐specific rabies virus outbreaks that lead to the deaths of livestock and economic losses. Based on incidental captures in our study area, which is an area of high cattle mortality from D. rotundus transmitted rabies, it appears that D. rotundus are being caught regularly in areas and elevations where they previously were thought to be uncommon. Our goal was to investigate demographic processes and genetic diversity at the north eastern edge of the range of D. rotundus in Mexico. We generated control region sequences (441 bp) and 12‐locus microsatellite genotypes for 602 individuals of D. rotundus. These data were analyzed using network analyses, Bayesian clustering approaches, and standard population genetic statistical analyses. Our results demonstrate panmixia across our sampling area with low genetic diversity, low population differentiation, loss of intermediate frequency alleles at microsatellite loci, and very low mtDNA haplotype diversity with all haplotypes being very closely related. Our study also revealed strong signals of population expansion. These results follow predictions from the leading‐edge model of expanding populations and supports conclusions from another study that climate change may allow this species to find suitable habitat within the U.S. border.

where it has since evolved into unique viral lineages within many bat species, including D. rotundus (Streicker et al., 2010). Significant economic damage is caused by D. rotundus rabies virus variants through loss of livestock throughout Latin America and the need for prophylaxis in human and domestic animals (Acha, 1967;Acha & Malaga-Alba, 1988;Anderson et al., 2014; see review in Streicker et al., 2010). Thus, any range expansions of rabies virus vectors, particularly D. rotundus, could lead to increased transmission of the disease and subsequent economic impacts (Anderson et al., 2014).

Desmodus rotundus ranges throughout most of South and Central
America and into northern Mexico in North America (Greenhall, Joermann, & Schmidt, 1983;Wilson & Reeder, 2005). Within the past 5 years, D. rotundus has been recorded within 50 km of the United States (U.S.) border in the Mexican state of Tamaulipas (Méxican Secretariat of Agriculture, Livestock, Rural Development, Fisheries, and Food (SAGARPA), unpublished data), which is north of their documented range limit. The fossil record shows that between 30,000 and 5,000 years before present (yBP), D. rotundus ranged within the current borders of the United States (Ray, Linares, & Morgan, 1988).
Fossil records also indicate the presence of multiple Desmodus species (D. stocki, D. archaeodaptes, D. draculae, and D. rotundus) from deposits in present-day New Mexico, Arizona, California, Texas, Florida, and West Virginia (Czaplewski & Peachey, 2003;Ray et al., 1988). It has been suggested that the northern distribution of D. rotundus is limited by the 10°C January minimal isotherm and that temperature and the bats' rate of metabolism are correlated (Lee, Papes, & Van Den Bussche, 2012;McNab, 1973). The 10°C minimal isotherm corresponds to areas in the northern hemisphere that, on average, experience minimum daily temperatures of 10°C during the month of January. It is likely that climate change will push the 10°C minimal isotherm northward and thus shift the northern limits of the D. rotundus distribution (Walther et al., 2002).
As described under the leading-edge model (Hewitt, 1993(Hewitt, , 1996(Hewitt, , 1999(Hewitt, , 2000, populations that have expanded rapidly from a core population are expected to experience a loss of genetic diversity due to genetic drift operating through repeated bottlenecks or founder events. Under this model of range expansion, population expansion occurs through long-distance dispersal and subsequent exponential population growth (Hewitt, 1999(Hewitt, , 2000. Such populations may have relatively little variation within the region of population expansion (Hewitt, 2000;Ibrahim, Nichols, & Hewitt, 1996). This model characterizes genetic diversity of many temperate species populations that expanded into newly suitable habitat after the retreat of glaciers (e.g., Arbogast, 1999;Lessa, Cook, & Patton, 2003;Stone & Cook, 2000 (Ditchfield, 2000;Martins et al., 2007Martins et al., , 2009. Some of these studies conclude that based on the distribution of genetic diversity across the sampled range of D. rotundus, that this species may have existed as isolated, refugial populations in South America during the Pleistocene (Martins et al., 2007(Martins et al., , 2009. However, none of these studies (Ditchfield, 2000;Martins et al., 2007Martins et al., , 2009  was thought to be uncommon (personal communications of local producers in Tamaulipas State) suggest that this process may be in progress. The potential establishment of populations in the US is of concern given that a rabies epidemic in livestock (vampire-bat variant) has expanded from the southeast to northwest across the states of San Luis Potosí and Tamaulípas, Mexico since 2001 (SAGARPA).
Population genetic analyses can be used to test predictions from the leading-edge range expansion model and provide an additional test of the potential for D. rotundus to expand its geographic distribution across the Mexico/U.S. border.
In this study, we use population genetic approaches to examine populations of D. rotundus at the north eastern edge of the species' distribution. We evaluate signals of demographic and spatial population expansion and levels of genetic diversity. We assessed demographic events with more rapidly evolving markers than used in previous studies of D. rotundus (Ditchfield, 2000;Martins et al., 2007Martins et al., , 2009), by investigating sequence variation in the mitochondrial DNA (mtDNA) hyper-variable I segment of the control region (HV1) and at 12 nuclear microsatellite loci (Piaggio, Johnston, & Perkins, 2008). If the leading-edge model accurately describes D. rotundus populations at the north eastern edge of its range, then we will find low genetic diversity and signals of recent population expansion.
Purified products (1 μL) were used in 10 μL cycle sequencing reactions for each primer (P* and E) with 1 μmol/L primer, 0.25 μL Sequences and genotypes were run on either an ABI 3130 or an ABI 3130xl genetic analyzer and visualized using GenemAPPer v.4.0 software (Applied Biosystems Carlsbad, California USA). Sequences were reviewed, edited, and aligned with sequencher v.5.2 (Gene Codes Ann Arbor, Michigan USA). Microsatellite genotypes were converted from GenemAPPer into GenePOP files for downstream conversion into appropriate input files for various analytical software through GmcOnvert (Faircloth, 2006).

| Statistical analyses-DNA sequences
Network analyses were performed to assess relationships among D. rotundus mtDNA haplotypes using netwOrk v.4.613 (available at http://www.fluxus-engineering.com/sharenet.htm). The number of haplotypes was assessed, and a haplotype network was calculated using a median-joining approach (Bandelt, Forster, & Röhl, 1999) with default settings including equal weights for transitions/transversions and Ɛ = 0. dIveIn (Deng et al., 2010) is an online, maximum-likelihood F I G U R E 1 Map of sampling area. Colors denote two separate sampling efforts, the size of the circle indicates sampling effort at that site. In blue, samples were collected from 2003 to 2005, samples were from the states of Tamaulípas, and Nuevo Leon were obtained as wing biopsies from bats captured using mist nets, harp traps, or hand-held nets at roosts and corrals (n = 288). In the second sampling effort (yellow), tissue samples were obtained between 2006 and 2010 during state and federal control efforts for D. rotundus in the states of San Luis Potosí and Tamaulípas (n = 314). The black star indicates where samples from another study (Romero-Nava et al., 2014) overlapped our sampling area; they collected 12 bats from Ciudad Valles, San Luis Potosí, Mexico tool that was used to calculate average sequence divergences among unique haplotypes from both the consensus and the most common haplotype. This program also was used to identify phylogenetically informative sites among sequences. In an expanding leading-edge population, we expect to observe low genetic diversity (Parsons et al., 1997) and a starburst-shaped haplotype network (Russell, Pinzari, Vonhof, Olival, & Bonaccorso, 2015) at a rapidly evolving locus such as HV1.
We used various approaches to test the hypothesis that leadingedge populations harbor signals of population expansion in DNA sequences (Hewitt, 1999(Hewitt, , 2000. The neutrality tests Fu's (1997) Tajima's (1989) D were used to test for population growth, with significance assessed via 10,000 simulated samples using ArlequIn v.3.5.1.2 (Excoffier & Lischer, 2010). A mismatch distribution was also generated in ArlequIn with 10,000 bootstrap replicates to test the hypothesis that polymorphisms among sequences fit a model of rapid demographic or spatial expansion (Rogers & Harpending, 1992;Slatkin & Hudson, 1991).

| Statistical analyses-microsatellite genotypes
Using the entire dataset of 12-locus genotypes for 602 D. rotundus individuals, null allele frequencies were calculated with a 95% confidence interval (Brookfield 1;Brookfield, 1996) Nei, 1987). Tests of HWE were evaluated using Bonferroni's corrections for multiple tests (Rice, 1989). GenAlex v.6.5 (Peakall & Smouse, 2006) was used to examine diversity across loci, including the Information Index (I), which estimates diversity of alleles while accounting for evenness in the representation of those alleles across the samples and correcting for uneven sample sizes (Smouse, Whitehead, & Peakall, 2015). A low value of I suggests that a few alleles dominate while most are rare, which is expected during a population expansion/founder event or the emergence of cladogenesis (Smouse et al., 2015). fstAt v.2.9.3.2 (Goudet, 1995) was used to calculate F IS and assess genotypic disequilibrium among all individuals and across all loci. Two tests of isolation by distance (IBD) were performed: one among roosts and another among individuals (Rousset, 2000). Arle-quIn was also used to assess IBD using a Mantel test for correlation between Slatkin's linear F ST and the log of the straight-line distances among the 18 roosts where we caught ≥ 6 individuals. For this analysis, we considered only individuals captured at the entrance of or within roosts (N = 394) rather than bats caught on the wing. A test of IBD by nonroost capture localities may pool individuals belonging to different populations and thus be susceptible to bias from the Wahlund (1928) effect, which would underestimate IBD. However, we were concerned that the IBD analysis between roosts may not adequately describe the relationships between genetic diversity and geographic distance across the sampling region due to our use of a subset of the total data, so we chose to also assess IBD across the entire sampling region and with all genotyped samples. IBD was calculated among all individuals (N = 602) in pairwise comparisons. This test uses log-transformed geographical distances, which are tested for correlation to probability of identity of genes between individuals (estimator â; Rousset, 2000).
We used bAPs v.6.0 (Corander & Marttinen, 2006;Corander, Marttinen, Sirén, & Tang, 2008) to estimate the number of genetic clusters (K) distributed among samples, and among females separately, informed by spatial information about capture locality (Corander, Sirén, & Arjas, 2008). We assessed population structure among females alone to test for sex-biased genetic clustering. The program was run as an assessment of spatial clustering of individuals with five replications of each value of K from 1 to 20. structure v.2.3.4 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) also performs a Bayesian analysis to estimate genetic clustering but uses HWE as an optimization criterion, whereas bAPs seeks to construct genetically homogeneous clusters (François & Durand, 2010). structure was run with all D. rotundus genotypes and set with a burn-in of 50,000 and 300,000 MCMC repetitions after burn-in. The ancestry model was set to admixture, and allele frequencies were assumed to be correlated among individuals. Values of K from 1 to 20 were explored, with five replications of each K to assess the stability of results between runs (Waples & Gaggiotti, 2006). Population demographics were evaluated from microsatellite data by calculating the imbalance index β (Kimmel et al., 1998). This estimator (β = 1 stable populations; β > 1 recent expansion or recovery from previously reduced population; β < 1 recent expansion from stable population) is based on variance in repeat numbers and allele frequencies and was calculated with 95% confidence intervals using SAS package and a program written and shared by T. Lehmann (Donnelly, Licht, & Lehmann, 2001).

| Genetic Diversity
We generated HV1 sequences of 441 bp (GenBank Accession #-) and 12-locus genotypes for 602 D. rotundus individuals. The 602 sequences comprised only 34 unique haplotypes, which all were closely related to each other (Figure 2). Sequence divergences calculated from the most common haplotype and from the consensus sequence of all haplotypes were very low (both = 0.5%) with only 14 segregating sites among the unique haplotypes, which included 12 transitions and two single-bp indels. Pairwise differences among these haplotypes as a proportion of sequence length ranged from 0.002 to 0.016, with an average of 0.007. The two sampling efforts combined for this study were largely distinct geographically ( Figure 1) and temporally. However, all the haplotypes, from both studies, were closely related with limited geographical association evident (Figure 2).
Evidence of possible null alleles were detected at two microsatellite loci: Dero_G10F_B03R (frequency = 0.07) and Dero_H02F_C03R (frequency = 0.19). Because the frequency of potential null alleles was low (Dero_G10F_B03R) to moderate (Dero_H02F_C03R), we kept these loci in the analyses as they were not expected to bias our estimates because they occurred at < 0.20 (Dakin & Avise, 2004 Index was low to moderate across loci (I = 0.56-2.57) with an average I = 1.37 ± 0.16, suggesting that most loci are characterized by a few very common alleles and more rare alleles. This pattern is evident when graphs of allele frequencies per locus are observed (Appendix S1). F IS across all loci was 0.08 (p = 0.05), which we interpret as being not biologically meaningful. No significant genotypic disequilibrium was detected in pairwise comparisons between all loci.

| Population structure
Microsatellite-based tests of IBD considering genetic and geographic distance between roosts were significant (r 2 = 0.21; p < 0.05). An analysis of IBD between individuals for the full dataset was highly significant (p < 0.00001). Bayesian clustering analyses in both bAPs and structure provided the strongest support for a single genetic cluster (K = 1) among all individuals (Figure 3). Although the average lnL scores for K = 2 were not significantly different than for K = 1 in structure results, bAPs showed that K = 1 was the best explanation of genetic clustering across the samples. Furthermore, with two or more clusters, assignment of individuals often was apportioned into more than one cluster, and there was not geographical thus not biologically meaningful coherence. Separate analyses of females in bAPs also resulted in the most likely clustering as K = 1. To test that the null alleles and subsequent HWE violations at Dero_G10F_B03R and Dero_H02F_C03R did not bias our results, we reran two analyses that assume HWE without these loci. First we estimated F IS (0.03; p = 0.05) and then we ran structure (K = 1, average lnL = −16974.7; K = 2, average lnL = −17177.4), which showed stronger support for a single cluster after removing loci with null alleles.
As in Bayesian clustering analyses, we found no indication of separate clustering of individuals based on shared alleles in the FCA analysis. FCA results (Figure 4) show a single tight cluster; thus, there no significant pattern of allelic distributions among individuals.

| Molecular demography
Desmodus rotundus from north eastern Mexico are characterized by a shallow evolutionary history and low genetic diversity with evidence of rapid population expansion. The starburst pattern of the mtDNA network is characteristic of rapid expansion evidenced by only a few founder haplotypes, followed by the recent evolution of closely related haplotypes (Slatkin & Hudson, 1991). This demographic model is further supported by a significantly negative Fu's F S , mismatch distribution analyses, and estimates of Kimmel's β > 1, showing significant recent population growth.
The control region of the mitochondrial DNA genome is known to accumulate mutations more rapidly than cytochrome b (Meyer, 1993;Parsons et al., 1997). With only 34 haplotypes ( Figure 2)  F I G U R E 3 Log-likelihood values (mean ± SD) for genetic clustering analyses of microsatellite data using structure. Loglikelihood values (mean ± SD) for K = 1 through K = 10 from structure analyses are in black; those from the top 10 values across runs of K from bAPs, which was only K = 1 and K = 2 are in red F I G U R E 4 Factorial Correspondence Analysis plot of categorical data exploring correspondence between individuals and alleles. Tight clustering of individuals and alleles indicates low diversity and a lack of population structure. The initial sampling effort is indicated in yellow, while the second sampling effort is indicated in blue T A B L E 1 The number of alleles per locus, observed heterozygosity (H O ), expected heterozygosity (H E ), and estimated null allele frequencies (Brookfield 1;Brookfield, 1996) (Martins et al., 2007(Martins et al., , 2009 found low sequence divergence among localities separated by over 400 km, but also found larger genetic differences over short distances. Instead of simple isolation by distance, these authors found low sequence divergence within an ecoregion and large divergences between ecoregions. Therefore, samples collected across short distances but in different ecoregions showed high divergences. Our samples were collected in a larger area than some of the ecoregions described by Martins et al. (2007Martins et al. ( , 2009 and spanned multiple ecoregions in Mexico (Olson et al., 2001), yet we found extremely low mtDNA sequence divergence across our entire study area. This pattern is explained by the leading-edge model where broadly distributed species harbor more genetic diversity in the core of their range (e.g., Ditchfield, 2000;Martins et al., 2007Martins et al., , 2009) than in the leading edge of their range where diversity is expected to be influenced by range contractions and expansions, founder events, and resulting population bottlenecks (Hewitt, 1993(Hewitt, , 1996(Hewitt, , 1999Hewitt 2004 (Grant, 2015). However, examining the mtDNA sequences and microsatellite data reveals that there is a lack of diversity at both markers and loss of intermediate-sized alleles (Appendix S1), which are more likely due to rapid demographic expansion. to thank Matt Hopken and Amy Gilbert for their comments and suggestions that greatly improved the manuscript.
led the writing.