A preliminary assessment of population genetic structure of the common vampire bat (Desmodus rotundus) in Colombia

Introduction

Rabies virus (RABV) is a prevalent and highly impactful pathogen (Anderson et al., 2014; Velasco-Villa et al., 2017) predominantly spread from bats to other mammalian species in Latin America by the common vampire bat (Desmodus rotundus) (Barquez et al., 2015). RABV outbreaks in humans and livestock regularly occur in Latin America (Stoner-Duncan, Streicker & Tedeschi, 2014; Soler-Tovar & Escobar, 2025), where human infections regularly follow outbreaks in livestock (Meske et al., 2021). Spread of RABV among D. rotundus colonies is thought to be a function of individual dispersal and migration patterns between populations (Blackwood et al., 2013; Benavides, Valderrama & Streicker, 2016; Megid et al., 2021) patterns of which can be inferred from population genetic structure. An understanding of basic population genetics and geographic structure of D. rotundus is beneficial for local stakeholders to mitigate RABV’s impact on humans and animals. Assessments of populations structure along landscapes is further justified in areas with pronounced landscape heterogeneity, such as that found in Latin America. Landscape heterogeneity impacts the spatial distribution and overall health of wildlife species and infectious disease dynamics (Biek & Real, 2010; García-Morales, Badano & Moreno, 2013; Becker et al., 2019). Furthermore, assessing population genetic structure can bypass the need for lengthy telemetry, mark-recapture, or longitudinal serology-based assessments of bat movement that would otherwise be needed to elucidate possible transmission routs between populations (Blackwood et al., 2013; Mara, Wikelski & Dechmann, 2014; Schorr et al., 2014; Streicker & Allgeier, 2016; Megid et al., 2021). To assess how landscape heterogeneity may influence population structure of D. rotundus, we conducted a cross-elevational assessment of population genetic structure for this species in Colombia using nuclear microsatellite DNA markers.

Colombia is a country with large environmental gradients within a relatively small area, and has endemic circulation of RABV (Fernando et al., 2017; Cediel, 2020; Rojas-Sereno et al., 2022), which is regularly monitored by the Epidemiological Information and Surveillance System of the Instituto Colombiano Agropecuario (ICA) (Instituto Colombiano Agropecuario (ICA), 2019; Instituto Colombiano Agropecuario (ICA), 2020; Bonilla-Aldana et al., 2022). Bovine rabies in Colombia shows marked regional variation and evolving outbreak dynamics. Bonilla-Aldana et al. (2022) reported 4,888 confirmed cases between 2005 and 2019 in Colombia, with a significant decline after 2014 and a notable concentration of cases in the Department of Sucre (Caribbean Region, North Colombia). Other research has citied low vaccination coverage and implementation difficulties as potential drivers for RABV outbreaks (Álvarez, Buitrago & Sáenz, 2014). Across the last five decades, 2,858 outbreaks of RABV have been reported in Colombia in cattle alone (Soler-Tovar & Escobar, 2025). Furthermore, Colombia has experienced an expansion of D. rotundus-associated RABV into novel areas, and an increase in landscape change since the cessation of armed conflict in 2016 (Clerici et al., 2020; Rojas-Sereno et al., 2022; Bonilla-Aldana et al., 2022). As a result, there is a need for a population genetics assessment for vampire bats in this region, especially within the context of landscape heterogeneity. Elevation can be used as a valuable proxy for environmental variation, since changes in altitude often correlate with changes in numerous environmental factors, such as temperature, vegetation composition, and humidity, making elevation a useful metric for studying landscape heterogeneity impacts (González, Willig & Waide, 2013). We therefore evaluated genetic diversity, population structure and its spatial distribution among sampling locations across an elevational gradient in the Colombian Andes. Our objective was to assess the extent to which population genetic diversity and structure of D. rotundus populations varied across space and were associated with environmental factors, which could be used to infer RABV transmission dynamics. This information could be useful to target RABV surveillance and pre-exposure human and livestock vaccination efforts.

Materials and Methods

Experimental design

To assess spatial patterns of genetic diversity and structure in D. rotundus populations, and to test their association with environmental factors relevant to rabies virus (RABV) transmission, we implemented a multi-step population genetic approach. This included: (1) broad-scale field sampling across an elevational gradient in Colombia, (2) genotyping at highly polymorphic microsatellite loci, (3) quantifying within- and among population genetic diversity and structure, and (4) testing the relationships between genetic metrics and environmental variables including geography, elevation, climate, and landscape features at the locations of capture.

We utilized Colombia’s dynamic topography and ecological diversity to select sampling locations across an elevational gradient in the Andes Mountains, under the assumption that elevation can be used as a reliable proxy for a variety of environmental factors (González, Willig & Waide, 2013; Lyons & Kozak, 2019). Thus, sample sites were selected to encompass elevations from low (<500 m), to moderate (500–1,000 m), to high (>1,000 m) to ensure that D. rotundus individuals were captured from a diversity of environments. Sites were also selected to minimize potential spatial autocorrelation and to reflect home-range distances of D. rotundus (McNab, 1973; Trajano, 1996), and hence were spaced at least 20 km apart. Elevational categories (i.e., low, moderate, and high) were utilized during site selection, but not in subsequent analyses. That is, all assessments were made independent of elevational category so as to not skew any potential relationships identified in subsequent analyses.

We obtained tissue samples from 81 D. rotundus individuals from Colombia in RABV-endemic areas during June and July of 2022 and 2023 (Table S1). Complementarily, five samples from a previous study collected between 2019 and 2022 were provided by the Instituto de Investigación de Recursos Biológicos Alexander von Humboldt (Bogota, Colombia) for the purposes of this study (Table S1), totaling 86 D. rotundus individuals from Colombia. As a genetic comparator, we included six D. rotundus individuals from western Mexico (Piaggio, Johnson & Perkins, 2008) (Table S1, Fig. 1) to contextualize population divergence. Environmental variables used for downstream associations with genetic structure included elevation (m), isothermality (average diurnal temperature range in ^∘C divided by annual temperature range in ^∘C), minimum temperature of the coldest month (^∘C), annual average precipitation (mm), precipitation seasonality (coefficient of variance) and landscape features such as associated land use and vegetation cover (i.e., cattle density or vegetation indices present at the location of capture) for each sample site. Environmental variables were derived from the NASA Socioeconomic Data and Applications Center (SEDAC) (NASA Earth Data, 2022), Gridded Livestock of the World Database (Food and Agriculture Organization of the United Nations, 2022), and the WorldGrids Archived database (Hengl, 2018) at 1-km spatial resolution. Variables were selected based on their ecological importance to D. rotundus and RABV spillover risk as established in previous research (Lyman & Wimsatt, 1966; Zarza et al., 2017; Mantovan et al., 2022; Magalhães et al., 2023; Van de Vuurst et al., 2023; Jones et al., 2024). Variable temporal range of the climate and landscape data were from 2016–2020.

Bat capture and tissue collection

Desmodus rotundus individuals were collected using mist nets (12 × 2.6 m, mesh size 36 mm) for at least five consecutive hours after sunset targeting periods of peak activity during the waning moon phase (Zeppelini, Azeredo & Lopez, 2019). Animal handling followed protocols adhering to recognized bat handling guidelines of the American Society of Mammalogists (Sikes, 2016), as approved by the Institutional Animal Care and Use Committee at Virginia Tech (IACUC approval #21-138). Approval was also obtained from the Institutional Care and Use Committee of the Universidad de La Salle (IACUC approval #087). Our sampling protocol abided by legal regulations of Colombia, including the National Statute for the Protection of Animals, Resolution 8430 of the Ministry of Health: scientific, technical, and administrative standards for health research, the Code of Ethics for the professional practice of veterinary medicine, and Law 1774, amending the Civil Code, Law 84 of 1989, the Penal Code, and the Code of Criminal Procedure. Capture of bats and data collection in Colombia were carried out under environmental permit Resolución 1473 de 2014 of the Autoridad Nacional de Licencias Ambientales (ANLA) (i.e., UniSalle Wild Species Collection Framework Permit). Both capture-and-release sampling and lethal sampling were conducted to collect the tissues analyzed in this study. Tissue samples were collected via pectoralis muscle dissection from individuals that were lethally sampled as part of an independent bat virome study aimed at reconstructing the virosphere of bats in the tropics (IACUC approval #21-138). For released individuals, tissue samples were collected via a 1-mm wing punch. All individuals that were captured during netting were sampled without discrimination based on age class, sex, or location.

To minimize stress and potential harm for captured individuals, bats were removed from nets and handled gently using clean cloth bags and restrained carefully by trained personnel. For euthanasia, bats were carefully placed inside a sealed Ziploc-type plastic bag containing a cotton pad soaked with two mL of USP-grade isoflurane, positioned to prevent direct contact with the animal. Bats were continuously observed for signs of anesthesia (i.e., decreased responsiveness to stimuli, loss of righting reflex, and slowed respiration). Once deep anesthesia was confirmed via absence of reflexes, cardiac puncture was conducted to induce exsanguination, followed by cervical dislocation to ensure death (Sikes, 2016). Death was confirmed by the absence of heartbeat, respiration, and the presence of fixed, dilated pupils. All samples were preserved at −20 °C and stored long-term at −80 °C at the Universidad de La Salle in Bogotá, Colombia.

DNA extraction and microsatellite genotyping

DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen) and quantified with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific) (Qiagen, 2023) at the Universidad de La Salle in Bogota, Colombia. All extractions had a yield of 50 ng/µL or greater. For each individual, we then amplified 12 nuclear microsatellite loci previously developed by Piaggio, Johnson & Perkins (2008) using Promega GoTaq Gold PCR protocols (Promega, 2019). Each 25-µl PCR reaction consisted of 1.0 µl 25 mM MgCl₂, 5.0 µl of 5X colorless GoTaq Flexi buffer (Promega, Madison, WI, USA), 0.5 µl of 10 mM PCR nucleotide mix, 1.0 µl of 10 µM forward and reverse primers, 0.125 µl of GoTaq DNA polymerase, 2.0 µl of template DNA (<0.25 µg/25 µl), and 14.375 µl of nuclease-free water. PCR amplification was conducted using a T-100 thermocycler (BioRad, Hercules, CA, USA) with a profile consisting of initial denaturation at 95 °C for 10 min; followed by 35 cycles of 94 °C for 30s, 52 °C for 45s, and 72 °C for 45s; and a final extension at 72 °C for 10 min. Samples then were held at 12 °C for at least 15 min prior to being placed in storage at 2 °C (Romero-Nava, León-Paniagua & Ortega, 2014). PCR products were sent to the Cornell University Institute of Biotechnology Resource Center (Ithaca, NY) for fragment analysis. We used the Fragman package (version 1.0.9) (Covarrubias-Pazaran et al., 2016) in R statistical software version 4.1.0 (R Core Team, 2021) to automate peak calls for determining amplification fragment sizes for alleles at each of the 12 loci for each D. rotundus individual. Allele lengths were confirmed via a visual inspection of peaks for each locus in GeneMarker version 2.6.3 software (SoftGenetics, 2018). We used Brookfield (Brookfield, 1996; Campagne et al., 2012) and Van Oosterhout methods (Van Oosterhout et al., 2004; Van Oosterhout, Weetman & Hutchinson, 2006; The Earth and Life Systems Alliance, 2023) through the PopGenReport package version 3.1 (Adamack & Gruber, 2014) in R to identify and exclude loci with high null allele frequencies (>0.5), ensuring data reliability.

Population genetics analysis

To assess how genetic variation was structured across space and potentially linked to environmental factors, we conducted the following sequential analyses:

Population structure and clustering.

To infer genetic clusters without prior assumption, we applied two complementary approaches: Discriminant Analysis of Principal Components (DAPC) (Jombart, Devillard & Balloux, 2010) using sequential K-means of the number of populations or conformance of genotype frequencies to Hardy-Weinberg equilibrium (HWE) expectations (Jombart, 2015; Desvars-Larrive et al., 2019), optimized using Bayesian information criterion (BIC) for each increasing number of K. We also applied Bayesian Clustering via STRUCTURE software (Pritchard, Stephens & Donnelly, 2000) applying the log-likelihood approach of Pritchard, Stephens & Donnelly (2000) and the delta K (δK) method of Evanno, Regnaut & Goudet (2005) to identify the best-supported value of K. We applied the poppr package version 2.9.6 (Kamvar, Tabima & Grunwald, 2014) to conduct analysis of molecular variance (AMOVA) within individuals, among individuals within clusters, and between clusters. These methods allowed us to identify putative populations and assess admixture, forming the basis for all population-specific analyses.

Genetic diversity and differentiation.

For each cluster identified by DAPC, we used hierfstat version 0.5-11 (Jerome et al., 2022) and dartRverse version 2.0 (Mijangos et al., 2022) to calculate:

1. Allelic richness, observed and expected heterozygosity, and inbreeding coefficient (F_IS, (Weir & Cockerham, 1984).

2. Pairwise F_ST and G′_ST (Yang, 1998; Hedrick, 2005) as measures of genetic differentiation using dartRverse. Statistical significance was evaluated using 1,000 permutations, with p-values < 0.05 indicating significant genetic differentiation.

3. Deviations from Hardy-Weinberg Equilibrium (HWE) using the poppr package (Kamvar, Tabima & Grunwald, 2014)

Effective population size (N_e) was estimated using NeEstimator (Do et al., 2014) with the Waples correction (Waples & Do, 2010; Waples, Larson & Waples, 2016) to infer local demographic stability.

Isolation-by-distance and sex-based dispersal.

To evaluate whether genetic differentiation was driven by distance or environmental factors, we conducted a Mantel (1967) test using Euclidean distance, F_ST, and G′_ST in the PopGenReport and dartRverse packages (Adamack & Gruber, 2014; Mijangos et al., 2022). A Mantel test also was used to assess isolation-by-distance among our Colombia samples, and the relationship between the elevation of sample sites and genetic distance. Statistical significance was assessed using 10,000 permutations, with significance thresholds set at p < 0.05.

Additionally, sex-specific spatial autocorrelation analyses was conducted to test for fine-scale structure and sex-biased dispersal using F_ST of females vs males, and by comparing genetic distance and geographic distance of each sex (Banks & Peakall, 2012). Autocorrelation values were calculated using Euclidean distance matrices of genetic and geographic distances and visualized using the dartR package (Mijangos et al., 2022). We then estimated 95% bootstrap confidence intervals (100 bootstrap replications) around the r values for each sex at different distance-classes based on our sample sites. Assessment of dispersal by sex for each independent population was not possible due to limited power for Population Three, where only five females were captured. As a result, sex-biased dispersal was conducted using all pooled samples from Colombia.

Inference of migration.

To distinguish between historical migration and contemporary movement, we estimated equilibrium migration rates (m) using the F_ST formula (Wright, 1951; Zhivotovsky, 2015). Values were interpreted based on conventional thresholds (e.g., Nm >1 indicates historic migration between populations). Also, we identified putative recent first-generation migrants using the Bayesian criterion in GeneClass2.0 (Piry et al., 2004) with 10,000 simulated individuals and a significance threshold of p < 0.01. identifying individuals likely to have moved among populations (Piry et al., 2004). This dual approach allowed us to infer both long-term and recent connectivity among D. rotundus populations in Colombia.

Linking genetics and environment.

Finally, we tested for associations between genetic structure and environmental characteristics, both independently and simultaneously, using multiple t-tests and a multivariate analysis of variance (MANOVA) test. For the purpose of this comparison, the comparator group from Mexico was not considered. This allowed us to assess whether environmental variables differed between populations within Colombia without confounding the comparison with data from Mexico. Populations were compared for differences in elevation (m), isothermality (average diurnal temperature range in ^∘C divided by annual temperature range in ^∘C), minimum temperature of the coldest month (^∘C), annual average precipitation (mm), precipitation seasonality (coefficient of variance), cattle density (individuals per km²), human population density, and annual average enhanced vegetation index (EVI) (factors with potential ecological relevance to RABV exposure and transmission risk). We confirmed that the selected environmental characteristics varied across the elevational gradient used for sample site selection using a multiple linear regression model. A p-value < 0.05 was considered statistically significant. We utilized the rstatix package (version 0.7.2) and base R to complete these analyses in R statistical software version 4.1.1 (R Core Team, 2021; Kassambra, 2023).

Results

Populations cluster analysis

Analysis of Bayesian Information Criterion (BIC) for increasing numbers K of discriminate analysis of principal components (DAPC) clusters indicated greatest support for three clusters within the dataset (Fig. S2). This inference was supported by the results of STRUCTURE analysis, with delta K most strongly supporting presence of three clusters among the sampling sites (Fig. S3). The clusters were partitioned geographically (Fig. 2). Population One encompassed the Mexico samples provided by A. Piaggio (United States Department of Agriculture) used as a comparator. Two populations were distributed across Colombia, with most of Population Two occurring in high-elevation portions of Colombia, east of the central Andes Mountain passage (Figs. 2 & 3) (807.5 m average elevation of capture). Population Three, in contrast, was distributed across central and western portions of the study area, slightly lower in elevation (Fig. 2) (575.2 m average elevation of capture). One sampling site (Agua de Dios) with the largest sampling effort (40 individuals) was split between populations two and three. We used the three clusters identified from the DAPC analysis as hierarchical strata (i.e., populations), while the respective sampling sites were considered subpopulations to calculate descriptive genotypic metrics (e.g., allelic variation, heterozygosity, and genetic distances) across both sampling sites and inferred populations. Observed heterozygosity (H_o) was less than expected heterozygosity (H_e) for all sampling sites and populations (Table 1).

Genetic variation and migration between and among populations

Results of AMOVA analysis showed significant genetic variance (27.7%, p = 0.01) within the respective sampling sites (Table 2). We also identified significant variance between sampling sites within populations (5%, p = 0.01), and between populations in Colombia (14.9%, p = 0.01) (Table 2). Populations Two and Three from Colombia had similar mean numbers of alleles per locus (7.4 and 7.1, respectively) and allelic richness (4.9 and 4.5) (Table 1). The inbreeding coefficient F_IS for all populations was between 0.31 and 0.33, indicating an excess of homozygotes above Hardy-Weinberg equilibrium expectations (Campagne et al., 2012), suggesting inbreeding or family structure. Estimated effective population sizes for the populations in Colombia were 16.9 for Population Two and 41.6 for Population Three, but were unbounded (Table 1). High genetic distances (F_ST∼0.30, G’_ST∼0.80) were observed between Population One in Mexico and Populations Two and Three in Colombia (Table 3), as was expected since Population One was at a greater distance from the others, had a small sample size, and was utilized as a comparator. Genetic distances between Populations Two and Three in Colombia were comparatively low (F_ST = 0.12, G’_ST = 0.34). The estimated equilibrium rate of migration across an ecological time-scale (i.e., genetically effective migration) between Populations Two and Three was 1.83 individuals per generation. Estimates of genetically effective migration between sites indicated that there was historic migration (>1 individuals per generation) between sites that were closer together (<120 km) and located on either side of the Río Magdalena valley, Llanos Orientales, and Orinoquia of the central Andean corridor (i.e., Ibague and Chaparral, Tolima, Yopal, Casanare, Agua de Dios, Cundinamarca, and San Martín, Meta) (Table S2). GeneClass2 direct migrant assessment revealed 15 likely first-generation migrants. Ten of these likely migrants were captured in Population Three (Table S3) and the remaining five in Population Two. Migrants were captured at sites Agua de Dios and Puente Quetame in Cundinamarca, Coello and Piedras in Tolima, Yopal in Casanare, and at San Martín and Pipiral in Meta. Identifying first-generation migrants between local collection sites was not possible due to limited sample size from some such sites. F_STacross all females sampled for this study was higher (0.24) than for all males sampled (0.14). Only females at the lowest distance class (60.4) calculated for our samples demonstrated significant positive spatial autocorrelation (r = 0.31, p = 0.001) (Fig. 4). Males had positive spatial autocorrelation at larger distance classes (438.8, r = 0.06), but this was not significant per our 95% confidence interval (p = 0.08). Females had stronger fine-scale genetic structure than males at short geographic distances. Assessment of differential dispersal by sex for each independent population was not possible due to limited power for Population Three, where only five females were captured.

Table 2:

Analysis of molecular variance (AMOVA) between and among collection of D. rotundus populations and sampling sites analyzed.

Variance and significance values within and between both Colombia populations delineated by DAPC and subpopulation (i.e., sample site locations). Significant variance was present across all strata tested (i.e., p values < 0.05). Most covariance originated from within the samples.

Source	Degrees of freedom	Sum of squares	Mean square	Sigma	% Covariance	ϕ	pvalue
Between populations	2	164.93	82.47	1.28	14.89	0.15	0.01
Between sample sites within populations	20	239.34	11.97	0.43	5.01	0.06	0.01
Between individuals within each sample site	69	640.05	9.28	2.38	27.70	0.35	0.01
Within individuals	92	414.81	4.51	4.51	52.40	0.48	0.01
Total	183	1,459.14	7.97	8.61	100	–	–

DOI: 10.7717/peerj.20306/table-2

Table 3:

Metrics of genetic distance between populations.

F_ST (below diagonal) and G’_ST (above diagonal) metrics of genetic distance between populations identified by DAPC. Pairwise F_ST (Yang, 1998; Hedrick, 2005) was assessed using dartRverse. Greater genetic distance is signified by large numbers. The largest genetic distance was present between populations one and three. Less genetic distance was present between Colombian populations two and three, indicating greater gene flow between the two populations.

	Population 1	Population 2	Population 3
Population 1	–	0.78	0.83
Population 2	0.28	–	0.34
Population 3	0.31	0.12	–

DOI: 10.7717/peerj.20306/table-3

Discussion

This study’s identification of two genetically distinct but interconnected D. rotundus populations in Colombia provides a significant insight into the potential mechanisms of RABV spread. The relatively low genetic differentiation between these populations based on both the F_ST (0.12) and G′_ST (0.34) metrics (Table 3), despite considerable geographic distances (21.6–696.7 km) supports the inference of both long-term migration and recent dispersal, potentially along low-elevation corridors such as the Rio Magdalena valley. The presence of first-generation migrants, especially individuals originating in the east but captured in the west, implies that current movement patterns may mirror historical migration, sustaining connectivity across the region. The presence of two populations at one sampling location (Agua de Dios), which rests at the center of the Central Andean corridor, supports this interpretation.

Given that D. rotundus dispersal is known to be limited to relatively short range (10–20 km) (Trajano, 1996; Streicker et al., 2012b; Rocha & Dias, 2020), and short-range dispersal has been related to slow wavefronts of outbreaks of D. rotundus- associated RABV (Benavides, Valderrama & Streicker, 2016), the observed genetic connectivity suggests that even moderate dispersal rates may enable RABV to cross regional ecological and geographical barriers. This dispersal is particularly concerning in areas where livestock densities are high, as infected bats could transmit RABV across populations despite physical separation by elevation or landscape. Overall, the consistency of this dispersal with previous findings and the importance of low-resistance corridors strengthen the ecological inference that migration and, potentially, viral flow in D. rotundus is governed more by geographic permeability than strict distance or habitat type (Pinto, 2009; Das et al., 2019; Seetahal et al., 2024)

Allelic richness, a metric of genetic diversity (Foulley & Ollivier, 2006), also showed similar genetic diversity within the two Colombian populations (4.91 and 4.49) (Table 1). The small number of populations identified, coupled with the limited genetic distance between these two populations indicates that D. rotundus colonies in Colombia are interconnected by migration (Martins et al., 2009; Streicker et al., 2012a; Streicker et al., 2016). Other research on genetic structure of D. rotundus in French Guiana using microsatellites and mitochondrial DNA found a greater degree of historic migration at a similar geographic scale to that in our study (i.e., regional) when compared to larger-scale comparisons between North and South America (Huguin et al., 2018). The spatial scale and geographic locations of previous studies may account for any differences between our results and those for studies done in different locations.

Our estimate of the equilibrium rate of migration (m = 1.83) showed long-term connection between Colombian populations, and our direct migrant identification assessment identified many first-generation migrants within the sampled generation (n = 15 of 86 individuals assessed). Most first-generation migrants (n = 10) originated in the eastern population in Colombia (Population Two), but were captured in the western population (Population Three) (Fig. 4). All migrants were captured at lower elevations (<1,000 m). The demonstration of interconnected bat populations across wide regions despite presumed physical barriers (e.g., the Andes Mountains) implies that regional rabies control cannot be confined to localized hotspots. Surveillance systems should be designed to track RABV not just where outbreaks occur, but in regions of potential bat movement, especially along valleys and lowland corridors.

No significant differences in environment or climate of occupancy were apparent between the two inferred populations, nor were there differences in anthropogenic factors such as agricultural activity (i.e., cattle density) or human population density (p > 0.05). Greater genetically effective migration rates were inferred between sites surrounding the Río Magdalena valley of the central Andean corridor, with sites to the east (Population 2) having the most apparent migration (Fig. 4, Table S2). The two inferred populations were geographically separated by the Andes Mountains to the east of the valley (Fig. 4). Other studies of D. rotundus-associated RABV have found that this bat may disperse through low-resistance areas such as valleys (Benavides, Valderrama & Streicker, 2016; Streicker & Allgeier, 2016), a finding convergent with our results. As such, our results support the potential role of elevation or topography in shaping D. rotundus dispersal, and hence we suggest that livestock vaccination campaigns prioritize areas adjacent to high-migration zones, such as the Río Magdalena valley.

Previous research using mitochondrial DNA also has found low genetic divergence among D. rotundus populations over large distances in some regions (i.e., Central America and the Pantanal region of Bolivia), while other regions in the same study had large genetic divergence over short distances (i.e., Atlantic forests in southern Brazil) (Martins et al., 2009). The authors hypothesized that greater genetic divergence may be driven by differences in ecoregion, rather than by distance per se for this species due to its volant nature. These results are consistent with those of our study, as distance among sites did not relate significantly to genetic differentiation. Nevertheless, other studies using the same microsatellite loci as this study have found relatively low genetic differentiation across different ecoregions for D. rotundus at the leading edge of the species’ broader geographic range in Mexico (Piaggio et al., 2017). In contrast, another study in Mexico using methods similar to those of our study found a higher degree of genetic diversity across a larger sample size and geographic area (Romero-Nava, León-Paniagua & Ortega, 2014). The lack of consensus among previous studies, and the lack of significant geographic or landscape patterns between populations in this study, reflects variation of geographic population genetic patterning in D. rotundus across its range. Results from this study may support previous hypotheses of male-biased dispersal for D. rotundus (Martins et al., 2009; Streicker & Allgeier, 2016). Females showed positive spatial autocorrelation as lower distance-classes than males, and positive spatial autocorrelation for males was not significant even at greater distances (Fig. 4). Our finding of related males potentially being collected at locations much farther apart than females echoes results from previous studies (Streicker & Allgeier, 2016). Nevertheless, previous studies have assessed inbreeding coefficients (F_IS) by sex and by population to further assess male-biased dispersal (Streicker & Allgeier, 2016). Assessment of F_IS by sex and by population was not possible here, however, due to the limited number of females—just five—collected in Population Three. Dispersal of infected D. rotundus individuals could facilitate RABV transmission between subpopulations (Streicker et al., 2012b; Bakker et al., 2019; Becker et al., 2021; Horta et al., 2022).

A limitation of this study is the lack of RABV prevalence data for the individuals collected. Hence, we are unable to directly assess the effect of genetic structure of D. rotundus in Colombia upon RABV spillover. RABV spillover occurrence data are available from Colombia, and are collected via the Epidemiological Information and Surveillance System of ICA (Instituto Colombiano Agropecuario (ICA), 2019). This organization performs surveillance for RABV outbreaks in livestock and reports suspected cases from local-level stakeholders (Rojas-Sereno et al., 2022; Instituto Colombiano Agropecuario (ICA), 2023). Recent reports have identified a persistent concentration of cases in the Northern Caribbean region of the country, mostly in cattle, and an expansion of RABV cases into previously naive areas (Rojas-Sereno et al., 2022). Between 1982 and 2023, Instituto Colombiano Agropecuario (ICA) (2024) recorded 3878 rabies outbreaks in animals (89% in cattle), distributed across 31 of Colombia’s 32 departments. Cases peaked in 2014 (542 cases), but have declined to as few 43 cases reported in 2019 (Bonilla-Aldana et al., 2022). The Caribbean region remains at high risk, however, with constant impact in the departments Cesar, Magdalena, Bolívar, and Sucre (Bonilla-Aldana et al., 2022). In the Orinoquia region (east of the Andes Mountains), 13 spatial clusters were identified in Vichada, in addition to cases in Casanare and Meta (Bonilla-Aldana et al., 2022). During the 2005–2019 period, the departments with the most reported cases of RABV were Cesar (18.6%), Magdalena (11.1%), Arauca (9.3%), Casanare (8.5%), and Sucre (6.6%) (Bonilla-Aldana et al., 2022). In 2020, most bovine rabies outbreaks were concentrated in Sucre, indicating focalized viral circulation that requires continuous epidemiological monitoring (Instituto Nacional de Salud, INS). Our findings of broader interconnection between D. rotundus populations in Colombia across many of these regions could help explain how RABV is spread across the country, particularly from eastern locations with higher cattle density to other, less agriculturally dense areas (Rojas-Sereno et al., 2022). To build upon the assessment conducted in this study, future assessments with a broader sampling base, potentially in the Caribbean region where cases are concentrated, could connect D. rotundus populations genetics information with RABV spillover data and coupled antibody tests for field samples.

Our sample size was smaller than those of previous assessments (Romero-Nava, León-Paniagua & Ortega, 2014; Piaggio et al., 2017), though our samples were collected over a narrower time-range. Assessment of a greater sample size across a broader geographic area or more diverse environments may further elucidate the impacts of environmental variation upon population genetic patterning of D. rotundus in Colombia. Utilizing only amplified fragment-size analysis for scoring microsatellite markers may have limited our observation of genetic variation as well. Other methods, such as microsatellite genotyping by sequencing (SSR-GBS) or whole-genome sequencing for single nucleotide polymorphisms (SNPs), could have shown more detailed patterns of genetic structure (Tibihika et al., 2019). Future research in this system should work to increase the geographic scope of sampling to capture more ecoregions or fine-scale landscape factors that may shape genetic differentiation for D. rotundus populations. Future research should also couple genetic assessments with RABV seroprevalence assessment to more directly test hypotheses of RABV transmission driven by dispersal of infected D. rotundus individuals.

Conclusions

This study aimed to establish a baseline assessment of population genetic structure of D. rotundus individuals across an elevational gradient in Colombia, with the broader goal of evaluating how genetic connectivity may influence RABV transmission dynamics. Screening 11 microsatellite loci developed by Piaggio, Johnson & Perkins (2008) and applying standard population genetic analyses, we quantified genetic variation and assessed the geographic distribution of genetic clusters. Comparative analysis with a reference population from Mexico allowed us to contextualize our findings within a broader regional framework. Our results revealed substantial genetic connectivity among Colombian D. rotundus populations, even across geographically complex terrains. This finding supports the hypothesis that interpopulation bat movement contributes to the dissemination of RABV. While environmental and anthropogenic barriers showed limited impact, natural geographical features such as valleys appeared to provide important movement corridors. These results indicate that it is difficult to predict vampire bats’ gene flow and dispersal, especially within the context of predicting RABV transmission between populations, as has been noted in previous research (Streicker et al., 2012b; Blackwood et al., 2013).

Our results contribute to current knowledge of the regional distribution of genetic variation for D. rotundus. Our results demonstrate that D. rotundus populations in Colombia exhibit significant genetic connectivity, even across complex landscapes, supporting the hypothesis that interpopulation bat movement can facilitate RABV spread. These findings advance our understanding of the spatial dynamics of vampire bat populations and underscore the relevance of incorporating ecological and genetic data into public health planning. Effective rabies control strategies should consider the role of landscape connectivity in shaping bat movement and virus spread. A more detailed assessment of the role of elevation is warranted to fill remaining knowledge gaps; specifically, future studies should incorporate denser sampling along elevational gradients, finer-scale genotyping (e.g., SNP-based approaches), and integration of environmental niche modeling to better understand how elevation and landscape features influence population structure and viral transmission risk. While environmental and anthropogenic differences between sites occupied by the respective populations were minimal, geographical corridors such as valleys play a key role in movement. Such integrative approaches can support more ecologically informed public health policies and practices, including targeted vaccination, strategic land-use management, and interdisciplinary surveillance.

Supplemental Information