Subspecific affinities, determination of population boundaries, and levels of population connectedness are of critical importance for the development of management and conservation planning. We used variation at a mitochondrial locus and 5 biparentally inherited nuclear loci to determine partitioning of genetic variation of western big-eared bats (Corynorhinus townsendii) within and among caves occurring in a fragmented landscape of gypsum deposits in western Oklahoma. To accomplish this objective, we first performed a phylogenetic analysis based on the mitochondrial locus of western big-eared bats from a large portion of their range. This analysis indicated that western big-eared bats at the periphery of the distribution in western Oklahoma share phylogenetic affinities with the most geographically restricted subspecies, C. t. pallescens. Because C. townsendii is rare in Oklahoma and is listed as a species of special concern, this finding provides additional support for the continued protection of this species in Oklahoma. Within western Oklahoma, we failed to detect significant differentiation among any caves for the biparentally inherited microsatellite data. However, the mitochondrial locus exhibited significant levels of genetic differentiation among caves, with the highest level of differentiation occurring between caves within the disjunct distributions of gypsum (ΦST = 38.76%). Although a significant amount of genetic differentiation was detected between populations on the 2 disjunct distributions of gypsum deposits, Analysis with the program Migrate suggested high levels of asymmetric gene flow among some populations. Our results provide a greater understanding of the population dynamics of western big-eared bats on the periphery of their range and highlight the importance of continued monitoring and study of this taxon.
Understanding ecology, behavior, and population structure of a species is essential to its proper management and conservation (Moritz 1994). Unfortunately, these types of data are difficult to obtain for most species. It also is difficult to define populations a priori for most organisms, and if population delineations are incorrectly identified, false information about the biology of the species can result (Burland and Worthington Wilmer 2001; Pearse and Crandall 2004). Fortunately, recent advances in molecular genetic techniques and statistical analyses provide new approaches to obtain data critical for understanding the biology of species and thus improving management and conservation programs.
Because of their nocturnal and volant nature, it often is difficult to delineate population boundaries of bats (Mammalia: Chiroptera) a priori. Recent studies have shown that barriers to dispersal exist even though many species of bats are capable of flying long distances. The greater mouse-eared bat (Myotis myotis) is widely distributed throughout much of Europe and the Middle East (Arlettaz et al. 1997; Benda and Horáček 1995) and is capable of covering several hundred kilometers annually between summer and winter roosts (Horáček 1985; Paz et al. 1986). Furthermore, female greater mouse-eared bats make daily movements of up to 25 km between nursery and feeding territories during lactation (Arlettaz 1996, 1999). Despite those behavioral characteristics, the Gibralter Strait, about 14 km of water separating the Iberian Peninsula from the Maghreb Africa, is a complete barrier to gene flow for both male and female greater mouse-eared bats (Castella et al. 2000). Although small-scale microgeographic population structure has been documented for several species of bats (Burland et al. 1999; Kerth et al. 2000; Weyandt et al. 2005; Wilkinson and Flemming 1996), studies on many others are lacking (Burland and Worthington Wilmer 2001).
Corynorhinus townsendii (Chiroptera: Vespertilionidae) occurs in the majority of the western United States, Canada, and Mexico (Fig. 1). C. townsendii has an affinity for caves (particularly maternity colonies) but also will roost in mines, rock crevices, and some human structures (e.g., under bridges—Barbour and Davis 1969; Handley 1959; Humphrey and Kunz 1976; Kunz and Martin 1982). In the southern Great Plains, including Oklahoma, C. townsendii roosts predominantly in gypsum caves, in what has been described as small, isolated, relictual populations (Barbour and Davis 1969; W. Caire, in litt.; Handley 1959; Humphrey and Kunz 1976). In western Oklahoma, maternity colonies of C. townsendii generally consist of 17–40 individuals (Humphrey and Kunz 1976) and have never been found with >80 individuals as in Colorado and New Mexico. Male C. townsendii roost singly during the summer, but their roost habits vary during other times of the year (e.g., single or ≥60 males and females roosting together—Humphrey and Kunz 1976; Twente 1955). In some of the larger gypsum caves, other species of bats such as Brazilian free-tailed bats (Tadarida brasiliensis) and cave myotis (Myotis velifer) roost in colonies 2–4 times larger than those of C. townsendii.
Federal agencies list C. townsendii as a species of special concern or a sensitive species (United States Fish and Wildlife Service 1994). In western Oklahoma, C. townsendii is rare, and listed as a Tier II species of special concern (Oklahoma Department of Wildlife Conservation 2005). C. townsendii always may have occurred in low numbers in the southern Great Plains, where it settled during the western contraction of its range during the late Pleistocene (Humphrey and Kunz 1976). Humphrey and Kunz (1976) estimated that 14,000 C. townsendii occurred in the southern Great Plains but also noted that this was likely a 2-factor overestimation. A 2-year survey of gypsum caves in western Oklahoma found only 167 C. townsendii in 26 of 74 caves surveyed (many caves were revisited); only 3 caves had ≥30 individuals and 23 caves had only 1–8 individuals (W. Caire, in litt.).
Corynorhinus townsendii is relatively sedentary, with migration from summer roosts to hibernacula or transitory stops at gypsum caves occurring within the same area (Barbour and Davis 1969; Humphrey and Kunz 1976). Based on recapture of banded individuals in Oklahoma and Kansas, movements of C. townsendii were generally <1.6 km, with a maximum of 8 km; movements from nurseries to hibernacula in autumn averaged 11.6 km, with 1 individual moving 39.7 km (Humphrey and Kunz 1976). Other telemetry studies of C. townsendii documented movements ≤5 km, with rare movements >5 km (Adam et al. 1994; Clark et al. 1993; Fellers and Pierson 2002; Wethington et al. 1996).
Given breaks in gypsum formations in western Oklahoma (Fig. 1), coupled with the relatively sedentary nature of this species, populations of C. townsendii in western Oklahoma may be isolated from each other. Genetic variation is important for population viability, and loss of this variation can cause decreased fitness in natural populations, or even extinction (Frankel and Soulé 1981; Lande 1988; Saccheri et al. 1998; Westemier et al. 1998). Because C. townsendii is a species of special concern in Oklahoma, and occurs there at the periphery of its geographic distribution in a fragmented landscape (Fig. 1), it may be especially prone to genetic drift, particularly if the less-suitable habitat separating the disjunct gypsum deposits represents barriers to gene flow. Therefore, our objective was to ascertain levels of population structure and hence gene flow of C. townsendii in western Oklahoma. To address this objective, we first used mitochondrial DNA (mtDNA) sequence data to perform a phylogenetic analysis of populations of C. townsendii from throughout the range of the western subspecies (australis, pallescens, and townsendii) to evaluate subspecific affinities of C. townsendii from western Oklahoma. We then used mtDNA sequence data and nuclear microsatellite DNA data to elucidate intra- and intercolony structure and to assess population dynamics and boundaries separating populations of C. townsendii in western Oklahoma.
Materials and Methods
Sampling sites and collection
Mist nests were placed at entrances of 12 gypsum caves in western Oklahoma that had been searched for roosting individuals or colonies of C. townsendii. Mist nets were supplemented with fish netting draped over gaps around mist nets and other openings of the same cave to decrease escape by bats. Netting occurred from late May to mid-August in 2003 and 2004 and resulted in capture of C. townsendii at 8 of the 12 caves. The 8 caves occurred in 4 counties in western Oklahoma: Alabaster Caverns in Woodward County; Saloon, Nescatunga, and Cult caves in Major County; Ake's Cave in Blaine County; and Fink 1, Ratzlaff, and 3 Domes caves in Washita County (Fig. 1). Maximum distance between any 2 of the 3 caves in Major County was 4.6 km, and 3.8 km between any 2 of the 3 caves in Washita County. The 3 caves in Major County are separated from the 3 caves in Washita County by about 108 km (Fig. 1). Distance between Alabaster Caverns and the caves in Washita County was about 144 km. All other distances among the 8 caves were >38 km.
Tissue samples consisted of 3-mm-diameter punch biopsies (1 from each wing) taken from the distal one-third of the plagiopatagium after sterilization with ethyl alcohol (Worthington Wilmer and Barratt 1996). Wing biopsies were stored in lysis buffer (Longmire et al. 1997) until whole genomic DNA was extracted. Procedures for capturing, handling, and obtaining wing biopsies followed guidelines of the American Society of Mammalogists (Gannon et al. 2007), the Oklahoma State University Animal Care and Use Committee, and a previous study of Ozark big-eared bats (C. t. ingens—Weyandt et al. 2005). After handling, all bats immediately took flight and appeared normal. Isolation of DNA was performed using a phenol extraction and salt precipitation procedure (Longmire et al. 1997).
Mitochondrial DNA variation and phylogenetic analyses
Primers were used to amplify via polymerase chain reaction (PCR) about 485 base pairs (bp) of the hypervariable region I of the mtDNA D-loop (Wilkinson and Chapman 1991). Amplifications were performed in 50-μl volumes containing 50 ng of DNA, 50 pmol of each primer, 5.0 μl of Taq buffer, 10 mM of deoxynucleoside triphosphates, 1.5 mM of MgCl2, 1 unit of Taq DNA polymerase, and double-distilled H2O to the final volume. The thermal profile consisted of an initial incubation at 95°C for 2 min followed by 35 cycles of 95°C for 1 min, 55°C for 1 min, and 72°C for 1 min, followed by a final incubation of 72°C for 30 min. Double-stranded amplicons were cleaned before DNA sequencing using the Wizard PCR Prep DNA Purification System (Promega, Madison, Wisconsin). Amplicons were sequenced in both directions using Big-Dye chain terminators and a 377 Automated DNA Sequencer (Applied Biosystems, Inc., Foster City, California). Overlapping sequence fragments were pieced together using AssemblyLIGN 1.0.9 (Oxford Molecular Group, PLC 1998).
For phylogenetic analyses, we downloaded from GenBank 161 orthologous sequences of C. townsendii from the western portion of its range (subspecies australis, pallescens, and townsendii—Piaggio and Perkins 2005). To polarize character-state changes, we also downloaded from GenBank orthologous sequences from C. t. ingens (Weyandt et al. 2005). CLUSTAL X (Thompson et al. 1997) was used to generate a multiple sequence alignment of 169 sequences: 161 sequences from Piaggio and Perkins (2005), 4 from this study, and 4 sequences of C. t. ingens from Weyandt et al. (2005). That alignment subsequently was imported into MacClade 4.0 (Maddison and Maddison 2000) for visual inspection and determination of unique haplotypes using the REDUNDANT TAXA option. Although Piaggio and Perkins (2005:765) stated that “a total of 1091 bp of the control region were sequenced for all 445 samples,” after viewing the alignment of sequences in their samples of C. townsendii, it was clear that not all individuals had their entire control region sequenced. To maximize the number of D-loop sequences for our phylogenetic analysis, we only used those sequences that spanned the fragment of the D-loop that we amplified. After removing an additional 13 bp from the 5′-end of the aligned sequences (those were invariant among our samples as well as those of Weyandt et al. ), we ended up with a multiple alignment of 101 sequences: 93 from Piaggio and Perkins (2005), 4 from this study representing C. townsendii from western Oklahoma, and 4 sequences from Weyandt et al. (2005) representing C. t. ingens from eastern Oklahoma. The 93 sequences deposited into GenBank by Piaggio and Perkins (2005) that we were able to use in our analyses were from Arizona (n = 4), California (n = 1), Canada (n = 31), Colorado (n = 25), Mexico (n = 3), Nevada (n = 3), New Mexico (n = 4), Oregon (n = 11), Utah (n = 7), Washington (n = 1), and Wyoming (n = 3). GenBank accession numbers, general locality information, and subspecific designations are given in Appendix I. The final DNA sequence alignment can be obtained by contacting the corresponding author.
MODELTEST (Posada and Crandall 1998) determined that the most appropriate model of DNA sequence evolution and model parameters were TrN + I + G: base frequencies = 0.3382, 0.2966, 0.1322; NST = 6; Rmat = 1.0000, 53.5752, 1.0000, 1.0000, 38.4058; shape parameter (α) of the gamma distribution = 0.4720; and Pinvar = 0.3962. Those model parameters generated within a maximum-likelihood framework were incorporated into PAUP* (Swofford 2000) to create a neighbor-joining tree. Sequences from C. t. ingens were used as outgroup taxa and support for internal nodes was estimated using a bootstrap approach with 1,000 iterations. As an additional measure of intraspecific divergence, parsimony analysis was performed using the computer program TCS (Clement et al. 2000). The phylogeny of haplotypes was constructed in the most-parsimonious fashion; the parsimony criterion defined each step between haplotypes as a single mutational step. Multiple clades were generated when the number of steps required to unite clades exceeded the number of steps allowed by a 95% confidence interval (95% CI—Crandall and Templeton 1993).
A Bayesian coalescent framework was employed using the program BEAST version 1.4.5 (Drummond and Rambaut 2003a) to estimate time to most recent common ancestor for all 3 subspecies, under the assumption that those taxa last shared a common ancestor approximately 1,000,000 years ago (Piaggio and Perkins 2005). A Bayesian skyline plot (Drummond et al. 2005) was used to calculate a prior distribution on the parameter θ. A demographic assumption of constant population size was employed along with the HKY nucleotide substitution model and a strict molecular clock. Three independent Markov chain Monte Carlo chains of 30,000,000 generations were run with the first 10% discarded as burn-in to optimize scale factors, adjusting scale factors until effective sample size (Drummond et al. 2002) was >100 for all estimated parameters. The final run consisted of a single chain of 60,000,000 generations with the optimized search settings. The program TRACER version 1.3 (Drummond and Rambaut 2003b) was used to examine quality of the final run based on degree of mixing, effective sample size values, and 95% confidence intervals for all divergence dates estimated.
Intra- and interpopulational variation
Because populations of C. townsendii were difficult to delineate a priori, coupled with the small samples sizes at some caves, individuals were combined into groups based on cave proximity. From concurrent telemetry conducted by S. J. Smith and previous studies (Adam et al. 1994; Clark et al. 1993; Fellers and Pierson 2002; Humphrey and Kunz 1976; Wethington et al. 1996), known maximum nightly movement of C. townsendii was determined to be about 5 km. Individuals from caves separated by <5 km were combined into 1 study group. Thus, we considered the 3 caves in Major County and the 3 caves in Washita County as distinct populations and refer to the 4 study areas as Alabaster Caverns, Major County, Ake's Cave, and Washita County.
Arlequin version 3.00 (Excoffier et al. 2005) was used to determine intrapopulational haplotype diversity (h) and nucleotide diversity (π) and global tests of nondifferentiation among caves (or groups of caves). The extent to which mtDNA variation was partitioned among and between caves within the disjunct gypsum deposits (Fig. 1) was analyzed in a hierarchical manner by the Φ-statistics using the nested analysis of molecular variance (AMOVA) option in Arlequin. Significance of variance estimates for Φ-statistics was obtained using a randomization procedure with 10,000 permutations (Excoffier et al. 1992). Arlequin also was used to calculate pairwise ΦST among bats from the 4 study areas. Sequential Bonferroni correction was used to correct for multiple tests (Rice 1989).
We used Migrate 2.3 (Beerli and Felsenstein 1999, 2001) to obtain Bayesian nonequilibrium estimates of effective population sizes (θ = 2Nefμ) and rates of migration (Nefm), where Nef was the effective population size of females, μ was the mutation rate per site per generation, and Nefm was the number of immigrant females per generation. Migrate 2.3 relaxes assumptions of populations having the same size and symmetric migration rates, which are necessary for ΦST-derived estimates of Nm, but Migrate 2.3 assumes no large increase in population size over time. Because we had sequence data for only 7 individuals from Alabaster Caverns, we combined bats from Alabaster Caverns with those from Major County. Thus, this analysis was conducted with 3 groups, Alabaster Caverns + Major County, Ake's Cave, and Washita County.
The following genealogy settings were used: 10 short Markov chains each consisting of 500,000 genealogies that were sampled every 100 trees, followed by 1 long Markov chain with 10,000,000 genealogies sampled every 100 trees. At the beginning of each chain, 10,000 trees were discarded so that the next chain was not biased toward the previously estimated parameters. An adaptive heating scheme with 4 parallel chains and initial relative temperatures of 1.0, 1.2, 1.5, and 3.0 were used. We used a transition to transversion substitution ratio of 10, and base frequencies were estimated from the data. Migrate 2.3 was run 5 times, with the ΦST-based starting parameters and different random seed numbers. Estimates of θ and Nfm from the 5 runs were averaged and used as initial parameters for a final run. Profile likelihoods for all parameters were evaluated at 0.025 and 0.975 percentiles.
Microsatellite variation and analysis
As an additional measure of genetic variation and differentiation for the samples of C. townsendii collected from western Oklahoma, we performed a microsatellite analysis. Using primers isolated from Eptesicus fuscus (Vonhof et al. 2002), 5 dinucleotide microsatellite loci were amplified via polymerase chain reaction. Reactions were performed in a final volume of 15 μl containing 50 ng of genomic DNA, 10 pmol of each primer, 9 μl of ABI Prism True Allele PCR Premix (Applied Biosystems, Inc.), and 3.8 μl of double-distilled H2O. The thermal profile consisted of 12 min at 95°C followed by 35 cycles of 94°C for 30 s, and annealing temperatures of 40°C (EF14), 45°C (EF1, EF6), 53°C (EF21), or 55°C (EF15) for 45 s, extension at 72°C for 45 s, and a final 10-min incubation at 72°C. Upon completion of polymerase chain reaction, 1.5 μl of polymerase chain reaction product was added to 3.5 μl of loading mix, containing an internal size standard (ROX); 1.5 μl of that combination was loaded into a lane of a 6% acrylamide gel. To determine genotypes, polymerase chain reaction products were visualized with an ABI-377 Automated DNA Sequencer and GenScan 2.0 and Genotyper software (all from Applied Biosystems, Inc.). To minimize bias in scoring, 2 people scored all microsatellite data independently.
Average alleles across loci, allelic frequency, deviations from Hardy–Weinberg and linkage equilibrium, and observed heterozygosity (HO) and unbiased expected heterozygosity (HE) were performed using Arlequin. Arlequin also was used to perform both a global test of nondifferentiation and to test for population differentiation at the same hierachical levels examined for the mtDNA data using a nested AMOVA. Population differentiation and spatial structuring were assessed with the Φ-statistics, including pairwise ΦST (with 10,000 steps in the Markov chain). Sequential Bonferroni correction was used to correct for multiple tests of the same hypothesis (Rice 1989).
Bayesian analysis of population structure (BAPS version 3.1—Corander and Marttinen 2005) was used to test for genetically distinct populations or clusters of individuals. Mixture clustering was performed at the group and individual level (Corander et al. 2003). When performing individual-level analysis, we tested for the number of genetically divergent populations (i.e., K), but for group-level analysis, we designated suspected clustering in data using the same groupings as used in the AMOVA. At the individual level, BAPS was performed several times with integers defined for K, from 3 to 25, using different series of numbers, with the same numbers repeated several times within an analysis (the possibility of K = 1 was automatically considered). To prevent K from fixing at a local mode (Corander and Marttinen 2005), analyses were performed with different maximum K integers (e.g., K = 3–15 for 1 analysis and K = 3–9 for another analysis). The mixture analysis compared data from multiple tests by examining the natural logarithm of the marginal likelihood of the data (logml), which was used to determine goodness-of-fit for assigned clustering (Corander and Marttinen 2005). When mixture results were subjected to admixture analysis, we designated 5 as the threshold for the minimum size of a population.
To determine if levels of genetic variability were consistent with historical demographics, BOTTLENECK 1.2.02 (Cornuet and Luikart 1996) was used with the 2-phase model. The 2-phase model incorporated the infinite alleles model and the stepwise mutation model at 70%, 80%, and 90% stepwise mutation model and 30%, 20%, and 10% infinite alleles model, respectively. That analysis was performed for both the 4 study areas and considering caves in the disjunct (northern and southern) gypsum distributions (Fig. 1) as 2 separate populations; all analyses were performed with 5,000 iterations. The Wilcoxon sign-rank test was used to evaluate significance between HO and heterozygosity expected if the population was at mutation–drift equilibrium. Heterozygosity excess (compared with equilibrium) is expected of bottlenecked populations because number of alleles is expected to decrease faster than heterozygosity with decreasing effective population size (Cornuet and Luikart 1996).
Sample sites and collection
Ninety-six C. townsendii were sampled at 8 gypsum cave entrances in western Oklahoma. Individuals captured by cave were 7 from Alabaster Caverns (2 females and 5 males), 18 from Saloon Cave (14 females and 4 males), 5 from Nescatunga Cave (3 females and 2 males), 14 from Cult Cave (9 females and 5 males), 15 from Ake's Cave (10 females and 5 males), 13 from Fink 1 Cave (9 females and 4 males), 6 from Ratzlaff Cave (2 females and 4 males), and 18 from 3 Domes Cave (11 females and 7 males). Colonies of 30–50 individuals were found at Alabaster Caverns, Cult Cave, Ake's Cave, Fink 1 Cave, and 3 Domes Cave. Despite a relatively large colony of C. townsendii at Alabaster Caverns, few individuals were captured because the cave was large and had many openings. No colonies were seen in Saloon or Nescatunga caves, but these also were large caves and colonies of C. townsendii could have been in inaccessible areas. Ratzlaff Cave was thoroughly investigated and although no colonies were found, 2 separately roosting C. townsendii were observed. Of the 8 caves, Saloon, Ake's, and 3 Domes caves were the only confirmed maternity sites (all males captured were juveniles), but demography of C. townsendii in the other caves could not be definitively determined.
Mitochondrial DNA variation and phylogenetic analysis
Mitochondrial sequence data were obtained from 90 of the 96 individuals sampled from the 8 caves. Multiple sequence alignment of 101 individuals resulted in 470 aligned positions of which 132 were variable. The 90 samples of C. townsendii from this study revealed 4 unique haplotypes that were not shared with C. townsendii from any other locality. Weyandt et al. (2005) sequenced the same fragment of the mitochondrial genome for the Ozark big-eared bat and therefore naming of our haplotypes (E–H) followed their system. Representatives of haplotypes E–H have been deposited in GenBank (accession numbers EU700341–EU700344).
The minimum-evolution tree based on maximum-likelilhood distances resulted in a tree of length 0.77563 and bootstrap analysis revealed 3 well-supported clades (Fig. 2). One clade contained individuals from Durango, Mexico, and corresponded to C. t. australis of Piaggio and Perkins (2005). A 2nd clade contained individuals from central and southeastern Colorado, New Mexico, and our samples from western Oklahoma. Although Piaggio and Perkins (2005) did not have individuals from western Oklahoma, this clade corresponded to C. t. pallescens. The 3rd clade contained the remainder of individuals and corresponded to C. t. townsendii (Piaggio and Perkins 2005). Within the clade of C. t. townsendii, we recovered a strongly supported clade containing individuals from western Colorado and Wyoming and a smaller, but also highly supported clade, consisting of individuals from Arizona and Nevada (Fig. 2).
As an independent measure of phylogenetic relationships, we constructed a parsimony haplotype network using the computer program TCS, which resulted in 4 independent, unconnected networks (Fig. 3). Additionally, the following haplotypes could not be connected to any other haplotypes based on the parsimony criterion and 95% confidence limit of 9 steps: Nevada U2 (AY713704) and V2 (AY713714); Oregon X2 (AY713931); Mexico K3 (AY713589), T2 (AY713520), and L3 (AY713594); and Canada Z (AY713584, AY713587, and AY713588). One clade consisted of haplotypes from central and southeastern Colorado, New Mexico, and western Oklahoma and was in agreement with results described above, corresponding to C. t. pallescens (Figs. 2 and 3). The 4 haplotypes (E–H) from western Oklahoma are most closely related to each other and then next most closely related to haplotypes found in the southeasternmost county of Colorado (Baca County). This clade was well supported (bootstrap support = 75%) in the phylogenetic analysis based on maximum-likelihood distances (Fig. 2). The remaining 3 unconnected networks corresponded to individuals representing C. t. townsendii (Figs. 2 and 3).
Rooting our analyses with C. t. ingens, estimates of time to most recent common ancestor from the BEAST analysis indicated that most widely distributed subspecies, C. t. townsendii, diverged approximately 52,480 years ago (95% CI = 41,190–64,140 years ago), the geographically most restricted subspecies, C. t. pallescens diverged approximately 18,630 years ago (95% CI = 12,350–22,530 years ago), and, finally, the most southerly distributed subspecies, C. t. australis, diverged approximately 12,840 years ago (95% CI = 6,415–19, 880 years ago).
Microgeographic genetic variation
Of the 4 haplotypes (E–H) detected in western Oklahoma, haplotype H only occurred in samples from the 5 northern caves (Table 1; Fig. 4). Although haplotype E was detected in bats occurring in caves within both the northern and southern gypsum deposits, that haplotype was in low frequency in northern samples but was the dominant haplotype in bats from the southern gypsum deposit (Table 1; Fig 4). Haplotype diversity (h) ranged from a low of 0.203 for bats in Washita County to a high of 0.713 for bats in Major County. In general, haplotype diversity was lowest for bats in Washita County, and nucleotide diversity (π) was low for all groupings (Table 1).
Proportions of genetic diversity attributable to variation within caves, among caves within gypsum deposits, and between the 2 disjunct gypsum deposits were 48.98%, 12.26%, and 38.76%, respectively. Pairwise ΦST revealed significant levels of differentiation between bats collected in Washita County and bats in Alabaster Caverns, Major County, and Ake's caves on the northern gypsum distribution (Table 2). The global test for nondifferentiation revealed the same pattern as results from the pairwise ΦST (Smith 2006).
Maximum-likelihood estimates of effective population sizes of females and migration rates were generally stable over multiple runs of Migrate 2.3 (Table 3). Final estimates of θ ranged from 0.0009 for Ake's Cave to 0.0128 for the northern populations in Woodward and Major counties. When migration rates were interpreted in light of the distribution of gypsum deposits, we detected high levels of migration from the northern population in Woodward and Major counties into the population at Ake's Cave and a high level of migration from our southernmost population in Washita County into the northern populations at Woodward and Major counties. Interestingly, these migration routes were highly asymmetric, and there was essentially no movement of females between Ake's Cave and Washita County or from Washita County to Ake's Cave (Table 3).
Based on genotyping 96 individuals for 5 microsatellite loci, mean HO and HE were 0.714 and 0.736 for Alabaster Caverns, 0.654 and 0.732 for Major County, 0.693 and 0.712 for Ake's Cave, and 0.665 and 0.764 for Washita County, respectively. No loci were consistently out of Hardy–Weinberg or linkage equilibrium. Hierarchical analysis of genetic variation revealed no significant genetic structuring among the disjunct gypsum deposits (ΦST = 0.09%) or among populations within gypsum deposits (ΦST = 0.00%). Thus, the AMOVA identified variation among individuals as accounting for almost 100% of the total genetic variation. That lack of genetic structuring was further corroborated based on results of pairwise ΦST and the global test for nondifferentiation (Smith 2006).
Mixture clustering at the group level for all population groupings resulted in optimal partitioning of 1 cluster (logml-value = −1,638.2). Admixture for group-level clustering also revealed 1 optimal cluster using results from mixture analysis. Mixture clustering at the individual level produced logml-values that ranged from −1,493.8 to −1,501.4, with optimal partitioning of 11–13 clusters. Results with logml of −1,493.8 (with optimal partitioning of individuals into 13 clusters) were used for admixture analysis because it produced the best goodness-of-fit for the multiple runs, or the exponential of their absolute differences in logml-values (i.e., −1,493.8 is 108 times better than logml of −1,501.4, as measured in the posterior probability—Corander and Marttinen 2005). Admixture clustering at the individual level resulted in optimal partitioning of 9 clusters, but they did not adhere strictly to any predefined population grouping; instead, population assemblages contained individuals from multiple clusters or genetically divergent groups (Smith 2006). Finally, when data were partitioned by the 4 study areas or the 2 disjunct gypsum deposits, the test for a historic bottleneck detected marginal heterozygosity excess for bats from Washita County on the southern gypsum deposit (P = 0.047).
Our phylogenetic analyses support conclusions of Piaggio and Perkins (2005) that deep phylogenetic breaks exist among populations of C. townsendii throughout their range (Figs. 2 and 3). Piaggio and Perkins (2005) suggested that western big-eared bats from central and southeastern Colorado as well as New Mexico belong to the subspecies C. t. pallescens. Although they did not include samples from western Oklahoma, Piaggio and Perkins (2005) and Handley (1959) indicated that bats in western Oklahoma, Texas, and central Mexico belong to C. t. australis (Fig. 1). Our results clearly refute placing western big-eared bats from western Oklahoma in the subspecies C. t. australis and instead suggest that they belong to C. t. pallescens (Fig. 2). Thus, our results provide additional support that of the 3 western subspecies of big eared-bats (australis, pallescens, and townsendii), C. t. pallescens is even more restricted than previously believed and deserves additional investigation (Fig. 1).
Piaggio and Perkins (2005) suggested that the 5 subspecies of C. townsendii (australis, ingens, pallescens, townsendii, and virginianus) diverged from each other during the Pleistocene epoch and that the 3 subspecies of western big-eared bats (australis, pallescens, and townsendii) probably diverged from each other about 1.0 million years ago (mya). We tested the hypotheses that the 3 subspecies of western big-eared bats diverged from each other about 1.0 mya using a coalescent approach, rooting our tree with C. t. ingens and constraining the tree to 1.0 mya. Results of coalescent analysis suggest that C. t. australis, C. t. pallescens, and C. t. townsendii last shared a common ancestor approximately 52,000 years ago. Piaggio and Perkins (2005) hypothesized that C. t. pallescens, which ranges from New Mexico to central and southeastern Colorado, and now includes western Oklahoma, probably retreated to the Chihuahuan Desert during the Pleistocene. Results of our coalescent analysis cannot refute this hypothesis. Although uncertainty exists regarding the location and number of glacial refugia as well as the patterns of glacial cycles (Demboski et al. 1999; Lessa et al. 2003), the presence of 4 unlinked haplotype networks plus several haplotypes that could not be connected (Fig. 3) suggests that several isolated populations of C. t. australis, C. t. pallescens, and C. t. townsendii existed during the last glacial maximum (18,000–20,000 years ago) and presumably expanded their ranges with glacial retreat (Pielou 1991). Thus, not only do our data support the biogeographic hypothesis of Piaggio and Perkins (2005) regarding C. t. pallescens but our results suggests that a similar pattern was exhibited by C. t. australis and C. t. townsendii.
Microgeographic variation in western Oklahoma
Understanding levels of genetic variability and partitioning of genetic variation within and among populations provides important information on levels and characteristics of gene flow and is necessary for effective management strategies (Rossiter et al. 2000; Weyandt et al. 2005; Weyandt and Van Den Bussche 2007; Worthington Wilmer et al. 1999). Castella et al. (2000) documented that the 14-km-wide Gibralter Strait is a barrier to gene flow for both males and females in the greater mouse-eared bat, a species capable of long-distance flight. Although volant, not all species of bats are capable of long-distance flight. Burland et al. (1999) documented a pattern of isolation by distance among maternity colonies of brown long-eared bats (Plecotus auritus) in Scotland that were separated by 0.1–100 km and Weyandt et al. (2005) detected genetic differentiation of maternally inherited mitochondrial haplotypes among colonies of C. t. ingens separated by <20 km.
Results of our mtDNA analysis suggest low dispersal of female C. t. pallescens. Based on hierarchical analysis of variance, 38.76% of the total genetic variation was attributable to the disjunct distribution of gypsum deposits. Distribution of mtDNA haplotypes and rates of migration between populations further suggest that C. t. pallescens does not cross this uninhabitable area in a symmetric fashion. Haplotype H was only detected in individuals on the thin continuous strip of gypsum occurring from the Kansas border south into Blaine County (Fig. 1). Similarly, haplotype E, although found in low frequency (6 of 53 individuals) in bats collected in the northern part of our study area, was the dominant haplotype in bats from Washita County (33 of 37 individuals examined). Occurrence of haplotype E in a relatively small number of individuals collected on the northern gypsum peninsula could be due to the northward dispersal of a few individuals from Washita County. In fact, Migrate 2.3 results revealed that of all the possible pairwise comparisons of gene flow among populations, significant numbers of migrating females were only occurring from the northern populations in Woodward and Major counties into Ake's Cave in the central portion of the study area and from caves in Washita County in the south into the northern populations at Woodward and Major counties. The movement of a large number of females from Washita County into Woodward and Major counties, a distance of about 100 km, suggests that females are able to transverse this area between the gypsum deposits. Piaggio and Perkins (2005) suggest that ongoing telemetry studies in Nevada indicate that C. townsendii is more than capable of moving these large distances. In fact, aerial tracking documented a pregnant C. townsendii traveling >150 km in a single night of foraging (Piaggio and Perkins 2005).
Moreover, it is known that although C. t. pallescens prefers to roost in caves, these bats also will roost in mines, rock crevices, and human structures (Barbour and Davis 1969; Handley 1959; Humphrey and Kunz 1976; Kunz and Martin 1982). Therefore, Migrate 2.3 results were not completely surprising. Lack of significant gene flow between caves in Ake's Cave and Washita County as well as the movement of individuals from Ake's Cave into Woodward and Major counties was more surprising (Table 3). The most likely explanation for the lack of gene flow in these 2 situations is the small sample size representing Ake's Cave. We were only able to collect 10 individuals from Ake's Cave. Seven of those 10 individuals possessed haplotype G, which was distributed throughout the study area, but the remaining 3 individuals possessed haplotype H, which was only found on the northern peninsula of gypsum (Fig. 4).
When considering the nuclear microsatellite data, we failed to detect significant differentiation among any of the 4 study areas. Moreover, our analyses failed to detect genetic differentiation between the 2 disjunct strips of gypsum. These results are interpreted as indicating that a sufficient number of male C. t. pallescens are frequently crossing this stretch of about 100 km of uninhabitable area and thus homogenizing genetic variation between these 2 groups of mitochondrial-defined bats. As noted above, significant differentiation of mtDNA was detected between bats collected at caves in the northern peninsula of gypsum and bats collected in the gypsum deposit in Washita County. However, we did not detect differentiation to the extent that haplotypes were restricted to a gypsum deposit, as found on opposite sides of the Gibralter Strait for greater mouse-eared bats (Castella et al. 2000). Unfortunately, dispersal dynamics of male C. t. pallescens are unknown. Although C. t. pallescens may have limited dispersal because of morphological wing constraints, Piaggio and Perkins (2005) indicate that ongoing work in Nevada using radiotelemetry data and airplane tracking reveal that C. townsendii is capable of traveling greater distances than previously documented with radiotelemetry data alone.
Even given this uncertainty regarding the dispersal ability of C. t. pallescens, it appears that the less-suitable habitat separating the disjunct gypsum deposits in western Oklahoma may be serving as a filter to gene flow for C. t. pallescens. Morphometric analysis of C. t. pallescens in Oklahoma and Texas separated by breaks in exposed gypsum deposits paralleled our results (Smith and Tumilson 2004). Individuals found in caves in gypsum deposits extending from Kansas southeast into Blaine County differed morphometrically from individuals found in caves in gypsum deposits in Washita County (Smith and Tumilson 2004). Furthermore, C. t. pallescens captured southwest of Washita County and in Texas differed in size from those sampled in Oklahoma from Washita County and north. Despite our lack of samples of C. t. pallescens in caves southwest of Washita County (most of these are on private land and we were not able to obtain access to these caves), distance between gaps of suitable habitat between caves in gypsum deposits in Washita County and other gypsum deposits to the southwest are sufficient to hypothesize further genetic differentiation of females in our study area (Fig. 1).
Most genetic studies of bats have investigated the macrogeographic scale, but crucial information (e.g., population structure, dispersal patterns, and social structure) has been found in several microgeographic-scale studies of bat species (Burland et al. 1999, 2001; Kerth et al. 2002a, 2002b; Petri et al. 1997; Rossiter et al. 2000, 2002; Weyandt et al. 2005). Low genetic diversity in conjunction with heterozygosity excess relative to the average expected value at mutation-drift equilibrium suggested that a severe bottleneck occurred in bats in Washita County. The BOTTLENECK analysis can only detect a recent bottlenecked population (2–4 Ne—Cornuet and Luikart 1996), but little historical information exists for C. t. pallescens in western Oklahoma. The Oklahoma Grotto Club documented C. t. pallescens in only 22 of 200 caves in northwestern Oklahoma (W. Caire, in litt.). C. t. pallescens is sensitive to human disturbance, particularly maternity colonies, and has a propensity to roost near cave entrances (Clark et al. 1993; Humphrey and Kunz 1976). Thus, human disturbance could effectively contribute to the declines in C. t. pallescens and to a bottleneck, particularly in isolated populations, such as the one in Washita County. These human disturbances also may explain, at least in part, the high level of migration of individuals from caves in Washita County into caves in Woodward and Major counties (Table 3). Human visitation of caves in Washita County has been noted, and most nearby residents were familiar with the caves in their county (Central Oklahoma Grotto 1984). Caves that we examined in 2003 revealed more human graffiti, trash, and destruction in Washita County than elsewhere in western Oklahoma.
Low haplotype and nucleotide diversity provided support for a recent bottleneck or a founder event by 1 of a few mtDNA lineages in populations of C. t. pallescens in Washita County. Historical demography is needed to verify a bottleneck because many factors contribute to low genetic diversity and could mask a bottleneck (e.g., founding event, disease, predation, or stochastic events—Pimm et al. 1989; Ramey et al. 2000). Other ecological or behavioral metapopulation dynamics also could negatively affect the population of bats in Washita County (e.g., mating behavior, frequent extinction and recolonization, or fragmentation—Pimm et al. 1989; Slatkin 1987). Our analyses were a crucial 1st step, which highlights the need for ongoing monitoring of C. t. pallescens in Washita County. Only long-term analysis with ≥10 polymorphic loci could provide insight into historical demography and causes of a genetic bottleneck (Cornuet and Luikart 1996; Luikart and Cornuet 1998; Luikart et al. 1998).
Although understanding ecology, behavior, and population structure of a species is essential for proper management (Moritz 1994), it also is critical to establish the correct taxonomy of the organism. Such is the case with C. townsendii in western Oklahoma. In the most recent phylogenetic study of Corynorhinus, without examining individuals from western Oklahoma, Piaggio and Perkins (2005) suggested that C. townsendii occurring in western Oklahoma belonged to C. t. australis. However, our phylogenetic analyses clearly document that C. townsendii occurring in western Oklahoma shares its closest phylogenetic affinities with members of C. t. pallescens. This is a significant finding because based on the analyses of Piaggio and Perkins (2005), C. t. pallescens has the most restricted geographic distribution. Moreover, because C. townsendii in western Oklahoma is rare and listed as a Tier II species of special concern by the Oklahoma Department of Wildlife Conservation (2005), this more restricted geographic distribution adds further support for proper management and protection of this taxon in Oklahoma.
With regard to levels of connectedness among populations of C. t. pallescens in western Oklahoma, overall, we detected low levels of gene flow of females that were highly asymmetric. Thus, despite finding high levels of gene flow for males and no population structuring, this situation could produce autonomous populations or subpopulations (Avise 2000). Based on our analyses, maternity colonies of C. t. pallescens, especially at Ake's Cave in Blaine County and the caves in Washita County, are not expected to be recolonized in the event of severe decline or loss because of the low dispersal of females (Avise 2000). Our results consistently suggested that C. t. pallescens in gypsum deposits in Washita County was in the most need of additional study and possible protection. Gypsum deposits and associated caves occur in a relatively small disjunct pocket in Washita County, resulting in minimal habitat for the small isolated population of C. t. pallescens. Although the overall population has low genetic diversity, it is not clear if this resulted from a bottleneck, low fecundity, or high mortality of juveniles. Any of these could result in inbreeding depression, fixation of deleterious alleles, and negative metapopulation effects (Cornuet and Luikart 1996; Lande 1994; Luikart et al. 1998; Mills and Smouse 1994) for this population, and stochastic events could lead to rapid extinction before the population began to suffer from negative genetic-related scenarios (Mills and Smouse 1994). Furthermore, populations of bats in Washita County appear to be impacted the most by human interactions, based on the amount of graffiti and trash in the caves (W. Caire, in litt.) and results of our study suggest high levels of gene flow from these caves into caves in the northern Woodward and Major counties without receiving gene flow from any of our other study sites.
We used the same markers used by Weyandt et al. (2005) for endangered C. t. ingens in eastern Oklahoma, which allowed direct comparisons between the 2 subspecies of big-eared bats. Based on mtDNA, each study detected 4 unique haplotypes that were not shared. Haplotypes detected among C. t. ingens revealed considerably more genetic variation than those detected in C. t. pallescens, suggesting older haplotypes in C. t. ingens. Overall, within-cave haplotype diversity was generally higher for C. t. pallescens than for C. t. ingens. Interestingly, haplotype diversity within Washita County (0.203) was most similar to that detected in C. t. ingens (0.000–0.473—Weyandt et al. 2005). Low haplotype diversity in these situations probably reflects isolation of big-eared bat subspecies in Washita County and eastern Oklahoma. Mean numbers of alleles per locus were similar in both studies, but mean HOs and HEs were higher in C. t. pallescens than in C. t. ingens. Weyandt et al. (2005) found no support for a historic bottleneck affecting C. t. ingens, but we detected a slight heterozygosity excess, reflective of a potential population bottleneck in bats residing in Washita County.
Similar to the conclusions of Weyandt et al. (2005) for C. t. ingens, long-term monitoring of populations of C. t. pallescens, with special emphasis on the geographically isolated population in Washita County, is warranted. Caves southwest of Washita County extending into Texas should be examined for colonies of C. t. pallescens. Because of the disjunct nature of gypsum deposits in this region, comparable genetic data from bats in southwestern Oklahoma and Texas will enhance conservation and our understanding of population dynamics of C. t. pallescens on the periphery of its range.
We thank the Veterans Affairs Department of Vocational Rehabilitation and Employment for funding this research and D. Jackson, Director Vocational Rehabilitation and Employment, Oklahoma City. The Oklahoma Cooperative Fish and Wildlife Research Unit (Oklahoma Department of Wildlife Conservation, Oklahoma State University, United States Geological Survey, United States Fish and Wildlife Service, and Wildlife Management Institute cooperating) provided additional funding. We thank private landowners and park managers for access to their land, especially C. and G. Inman. Finally, we thank R. Froese, K. Phelps, and R. Allison for field assistance and S. Davis and S. Weyandt for assistance in the laboratory.
Specimens examined.—For each specimen, we provide the GenBank accession number; haplotype designations; the state or country; county (United States), state (Mexico), or province (Canada); and subspecies as determined by Piaggio and Perkins (2005). Based on results of our phylogenetic analyses (Figs. 2 and 3), Corynorhinus townsendii from western Oklahoma belongs to C. t. pallescens
Table 1.—Genetic diversity of Corynorhinus townsendii pallescens from 8 caves in western Oklahoma using mitochondrial sequence data (mtDNA) and nuclear microsatellite data (nDNA). Individuals were analyzed according to the 8 cave entrances where they were captured, or grouped (combining caves) based on known maximum nightly movement of C. t. pallescens (Major County and Washita County caves) and continuous gypsum distribution (Fig. 1). For the mtDNA data, number of individuals sampled (nmt), number of individuals with each haplotype, and haplotype diversity (h) and nucleotide diversity (π) are given; for the microsatellite data, number of individuals sampled (nms), alleles per locus averaged across loci (A), and mean observed heterozygosity (HO) and expected heterozygosity (HE) are given. Values are presented for individual caves as well as the grouping of caves in Major and Washita counties
Table 2.—Pairwise Φ-statistics using sequential Bonferroni correction (Rice 1989) of Corynorhinus townsendii pallescens in western Oklahoma based on mitochondrial DNA sequence data. Bold numbers designate statistically significant differentiation (P = 0.000)
Table 3.—Migration rates and population size estimates from Markov chain Monte Carlo simulations implemented in the computer program Migrate 2.3 (Beerli and Felsenstein 2001). θ is 2Nefμ, where Nef is female effective population size and μ is the mutation rate per site per generation. Nefm is the number of migrating females per generation. Northern is the grouping of bats collected in Woodward and Major counties, whereas Ake's refers to Ake's Cave and Washita refers to the population of bats collected in Washita County