Dispersal and gene flow are significant components of metapopulation structure and dynamics. We examined gene flow in the pygmy rabbit (Brachylagus idahoensis), a sagebrush specialist that occurs in small, isolated populations in the Great Basin region and is believed to have limited dispersal abilities. We examined genetic diversity, gene flow, and population genetic structure, and investigated the effects of putative barriers to gene flow in 8 sample locations within 2 study areas in Idaho during 2002–2006. Potential barriers included secondary roads, highways, creeks, and agricultural pastures. Distance between sample locations ranged from 0.5 to 32 km. We expected that gene flow would be limited and that moderate levels of genetic structure would be apparent among sample locations separated by barriers and distances over several kilometers. However, analyses based on 15 microsatellite loci from 239 rabbits revealed low to moderate levels of genetic differentiation among sample locations. Bayesian cluster analyses suggested that sample locations separated by <14 km comprised 1 genetic cluster with high levels of gene flow among sample locations within clusters. We documented no significant evidence of a sex bias in dispersal based on sex-specific FST among sample locations, and relatedness (r), mean assignment index (mAI), and variance of assignment index (vAI) within sample locations; however, spatial autocorrelation suggested that juvenile females might disperse farther than males. These results indicate that dispersal capabilities of pygmy rabbits greatly exceed previous estimates, and that creeks and roads were not effective barriers to gene flow in our study areas.
Dispersal influences patterns of gene flow and population genetic structure and contributes to metapopulation dynamics (for reviews, see Chepko-Sade and Halpin 1987; Clobert et al. 2001; Hanski and Gilpin 1997). Dispersal and gene flow also help maintain genetic diversity, allowing populations to respond to changing environments. Loss and fragmentation of habitat can lead to reduced gene flow among small, fragmented populations and a subsequent decline in genetic diversity within remaining populations (Frankham et al. 2002; Garner et al. 2005). Because genetic diversity is correlated with population fitness (Reed and Frankham 2003), understanding patterns of gene flow and dispersal is necessary for conservation of wildlife in fragmented landscapes.
Despite the importance of dispersal to maintaining genetic diversity and metapopulation dynamics, knowledge of dispersal patterns and capabilities is incomplete for many mammalian species. Mammals can exhibit both natal and breeding dispersal. Natal dispersal, the permanent movement from the natal area to the area where an individual 1st breeds, typically occurs more commonly and over longer distances than breeding dispersal, which refers to movement between successive breeding events (Greenwood 1980). Logistical constraints of field techniques traditionally used to study natal dispersal (e.g., radiotelemetry and mark–recapture) often limit collection of data and lead to biased estimates of juvenile movements (Koenig et al. 1996). Genetic techniques based on markers with different modes of inheritance are increasingly being used to infer patterns of dispersal and gene flow from population genetic structure or the identification of immigrants (Favre et al. 1997; Goudet et al. 2002; Mossman and Waser 1999; Paetkau et al. 1995; Prugnolle and de Meeûs 2002; Rannala and Mountain 1997). The growing suite of genetic tools can advance understanding of dispersal behavior and gene flow in mammals and improve our ability to conserve species occupying fragmented landscapes. Such techniques are especially useful for species for which ecology, demographics, or behavior limit the use of traditional field methods.
The pygmy rabbit (Brachylagus idahoensis) is a burrowing, sagebrush specialist distributed patchily, both historically and currently, throughout the Great Basin and adjacent intermountain areas (Green and Flinders 1980; Orr 1940; Weiss and Verts 1984). The pygmy rabbit is the smallest rabbit in North America (425–525 g) and is restricted to areas with tall, dense clumps of big sagebrush (Artemisia tridentata spp.) and soils deep enough to support burrow systems (Green 1978; Orr 1940; Wilde 1978). Because of their small size and specific habitat requirements, pygmy rabbits were believed previously to be relatively sedentary (Janson 1946; Wilde 1978). However, recent radiotelemetry studies of natal dispersal (Estes-Zumpf and Rachlow 2009) and adult home-range size and spacing behavior (Crawford 2008; Sanchez and Rachlow 2008) indicate that both adult and juvenile pygmy rabbits are capable of moving distances up to 10–12 km.
Despite the ability of pygmy rabbits to move farther and more frequently than previously estimated, our understanding of connectivity among populations is incomplete. Loss, fragmentation, and degradation of sagebrush steppe has led to increased concern about the effects of landscape change on this species (United States Fish and Wildlife Service 2003, 2008), especially following the rapid decline and subsequent extirpation of pygmy rabbits in the Columbia Basin of eastern Washington (United States Fish and Wildlife Service 2003). Currently, no information exists about fine-scale genetic structure or about gene flow among populations across natural or anthropogenic landscape features. However, continued loss and fragmentation of sagebrush habitat (Knick and Rotenberry 1997) makes information about connectivity among populations of pygmy rabbits critical to conservation planning (Rachlow and Svancara 2006).
We used genetic methods to evaluate patterns of natal dispersal and gene flow within and among sample locations in 2 study areas in Idaho. Sample sites were separated by different distances and putative barriers to gene flow (perennial creeks, rural highways, and agricultural expanses). We examined population genetic structure and gene flow within and among sample locations. We hypothesized that rural highways, perennial creeks, and agricultural expanses would impede gene flow, resulting in evidence of population substructure within each study area. We also examined evidence for sex-biased natal dispersal in pygmy rabbits. We hypothesized that pygmy rabbits, like many other mammals, would exhibit male-biased natal dispersal resulting in increased population genetic structure and isolation by distance among females relative to males. Because pygmy rabbits occur in small, fragmented populations, an understanding of how natural and anthropogenic features influence dispersal and gene flow among populations will contribute to conservation planning for this species in a changing landscape.
Materials and Methods
Study area
We collected genetic samples during 2002–2006 from 2 study areas, 1 in the Lemhi Valley in east-central Idaho, and 1 in the Camas Prairie in central Idaho. Both study areas are dominated by sagebrush-steppe vegetation. Land is managed by private, state, and federal (Bureau of Land Management and United States Forest Service) landowners. Public lands are managed for multiple uses, including seasonal cattle grazing. Grazing and agricultural production of crops occur on surrounding private lands.
The Lemhi Valley study area is located near Leadore, Idaho, and is bounded by the Beaverhead Mountains on the east and by the Lemhi Range on the west. Elevation ranges from 1,900 to 2,100 m. The study area is bisected by 2 rural highways, secondary paved and gravel roads, and several perennial streams (Fig. 1a). Dominant shrub cover on the study area includes big sagebrush (A. tridentata spp.) and green rabbitbrush (Chrysothamnus viscidiflorus). Density of shrubs varies across the study area (Sanchez and Rachlow 2008); however, much of the area is characterized by a mounded microtopography (mima mounds) with taller, denser shrub cover on mounds and lower, sparser shrub cover between mounds (Tullis 1995). We sampled pygmy rabbits from 5 sites delineated within the Lemhi Valley study area: Smokey Cub (SC), Cedar Gulch (CG), Rocky Canyon (RC), Hawley Creek (HC), and Warm Springs (WS). Sample locations were separated by up to 13 km and various putative barriers to movement by pygmy rabbits (roads, creeks, and agricultural expanses; Fig. 1a).
Fig. 1
Distribution of sample locations and putative barriers to gene flow for pygmy rabbits (Brachylagus idahoensis) in a) the Lemhi Valley study area in east-central Idaho and b) the Camas Prairie study area in central Idaho.
The Camas Prairie study area is located between Mormon and Magic reservoirs south of the Soldier Mountains (Fig. 1b). Elevation ranges from 1,500 to 1,750 m. The study area is characterized by relatively continuous big sagebrush and is bisected by a rural highway and several 2-track roads and perennial streams. Distribution of rabbits was patchy across the study area. We collected tissue samples from pygmy rabbits at 3 sample locations separated by 9–32 km (Fig. 1b): Mormon Reservoir (MOR), Johnson Hill (JH), and Magic Reservoir (MAG).
Data collection
We collected genetic samples from 209 pygmy rabbits from the Lemhi Valley during 2002–2006 and 40 rabbits from the Camas Prairie during 2004–2005. We trapped rabbits with Tomahawk live traps (14 × 14 × 40.5 cm; Tomahawk Live Trap Co., Tomahawk, Wisconsin) set at burrow entrances down which rabbits were observed to escape. Traps were checked approximately every 15 min. We restrained rabbits with a cloth handling cone (Koprowski 2002); recorded sex, age, standard morphometric measurements, and location coordinates; and used a 2-mm biopsy punch to collect ear tissue for genetic analyses. Tissue was stored in 100% ethanol until extraction. Techniques were approved by the University of Idaho Animal Care and Use Committee (protocol 2003-50) and are in compliance with guidelines published by the American Society of Mammalogists for use of wild mammals in research (Gannon et al. 2007).
The DNA extraction and microsatellite analyses were conducted at the Laboratory for Conservation and Ecological Genetics (University of Idaho, Moscow, Idaho). We extracted DNA from tissue samples following the Qiagen tissue protocol (Qiagen Inc., Valencia, California). We amplified 9 microsatellite loci (A2, A10, A121, A124, A133, D103, D118, D121, and D126) developed for the pygmy rabbit (Estes-Zumpf et al. 2008) and 7 microsatellite loci (Sol08 and Sol30 [Rico et al. 1994]; Sat05, Sat07, Sat08, and Sat12 [Mougel et al. 1997]; and Sol44 [Surridge et al. 1997]) developed for the European rabbit (Oryctolagus cuniculus) using polymerase chain reaction conditions described by Estes-Zumpf (2008) and Estes-Zumpf et al. (2008). We visualized microsatellite genotypes with an Applied Biosystems 3130 sequencer (Applied Biosystems, Foster City, California) and scored alleles using GeneMapper version 3.7 (Applied Biosystems). We tested for deviations from Hardy–Weinberg and linkage equilibrium at each locus within and across sample locations at each study area using GENEPOP version 3.1 (Raymond and Rousset 1995). Significance levels were adjusted for multiple comparisons using the sequential Bonferroni correction (Rice 1989). We estimated expected heterozygosity (HE) and observed heterozygosity (HO) using GENEPOP version 3.1, and allelic richness (AR) using FSTAT version 2.9.3.2 (Goudet 1995, 2001). To reduce bias resulting from samples of closely related individuals while acknowledging the naturally small population sizes and clustered distribution of pygmy rabbits, we removed all juveniles likely to be littermates. Criteria for identifying littermates included juvenile rabbits of similar body mass for which relatedness coefficients (r) suggested full siblings. Theoretically, the r-value for full sibling is 0.5; however, factors such as the amount of polymorphism at loci and the extent of the genome examined can influence r-values. Therefore, we considered juvenile rabbits of similar body mass with r ≥ 0.3 to be full siblings. Pairwise relatedness values among individuals were calculated separately for each sample location using Kinship version 1.3.1 (Goodnight et al. 2004). We randomly selected and regenotyped samples from 26 rabbits (10%), and we calculated the mean error rate per locus following Pompanon et al. (2005).
Mitochondrial DNA (mtDNA) from 75 rabbits from the CG, RC, and WS sample locations in the Lemhi Valley and 33 rabbits from the 3 sample locations in the Camas Prairie was sequenced at the Molecular Genetics Laboratory, Washington Department of Fish and Wildlife (Olympia, Washington). We sequenced a 509-base pair fragment of the mtDNA control region in both directions using primers H16498 and L15997 (Ward et al. 1991). Control region sequences were amplified in 15-µl reactions consisting of 2.5 mM of MgCl2, 0.2 µM of each primer, 0.5 U/µl of AmplitaqGold® DNA Polymerase (Applied Biosystems), and 1 µl of template DNA. Samples were amplified in a Tetrad 2 thermal cycler (MJ Research, Inc., Waltham, Massachusetts) by an initial denaturation (95°C, 10 min), 20 cycles of denaturation (95°C, 30 s), annealing (51°C, 30 s, decreasing by 0.2°C each cycle), and extension (72°C, 60 s), followed by 15 cycles of denaturation (95°C, 30 s), annealing (48°C, 30 s), and extension (72°C, 60 s), and a final extension at 72°C for 7 min. We used the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) for sequencing reactions. Sequences were analyzed using an ABI 3130xl sequencer with associated Sequencing Analysis software (Applied Biosystems) and aligned using Sequencher version 4.5 (Gene Codes Corporation, Ann Arbor, Michigan). We used TCS version 1.21 (Clement et al. 2000) to determine frequencies of unique haplotypes. Although nuclear DNA occasionally can be amplified along with or instead of the target mtDNA region, we detected no ambiguous or unusual sequences that would indicate this problem. mtDNA sequences amplified in this study also aligned well with mtDNA sequences from GenBank for other lagomorph species. Representatives of all haplotypes were submitted to GenBank (accession numbers GQ166913–GQ166921, GQ166924, GQ166939, and GQ166940). We estimated haplotype (h) and nucleotide (π) diversity using Arlequin 2.0 (Schneider et al. 2000).
Genetic structure
We examined gene flow within each study area by looking for evidence of population substructure among sample locations. We calculated pairwise FST-values for both microsatellite and mtDNA data among sample locations within each study area at the P < 0.05 level based on 10,000 randomizations in Arlequin 2.0 (Schneider et al. 2000). We also tested the hypothesis that mtDNA haplotypes were distributed randomly among sample locations using an exact test for population differentiation in Arlequin 2.0. We used microsatellite data to test whether allele frequencies differed between pairs of sample locations in each study area based on 1,000 randomizations in GENEPOP version 3.1. We determined the number of loci (of 15) showing significant allelic differentiation (P < 0.05) for each pair of sample locations. For the Lemhi Valley study area, we tested for evidence of isolation by distance (Wright 1943) using Mantel tests (Mantel 1967) of both linearized FST and Nei's standard genetic distance (DS—Nei 1978) versus log geographic distance using program zt (Bonnet and Van de Peer 2002). Nei's standard genetic distance was calculated using SPAGeDi, version 1.2 (Hardy and Vekemans 2002).
We also used 3 Bayesian clustering analyses (2 nonspatial and 1 spatial) to determine the number of populations represented by our sample locations in each study area. Bayesian clustering analyses group individuals into genetic clusters that minimize Hardy–Weinberg and linkage disequilibrium. We inferred the optimal number of genetic clusters (K) from 10 independent runs of each value of K for K = 1–10 for the Lemhi Valley and for K = 1–6 in the Camas Prairie. We 1st used the admixture model in STRUCTURE, version 2.2 (Pritchard et al. 2000) with a burn-in period of 50,000 and 500,000 Markov chain Monte Carlo repetitions. The optimum K-value was determined using criteria outlined in Pritchard et al. (2000) and by calculating ΔK (Evanno et al. 2005). We also used STRUCTURAMA (Huelsenbeck and Andolfatto 2007) to infer K, assuming that the number of populations and the expected prior number of populations were random variables with a gamma distribution of shape = 1 and scale = 1. We ran 100,000 Markov chain Monte Carlo repetitions with a sample frequency of 100 and a print frequency of 25. Lastly, we used the spatial option in BAPS 4.0 (Corander et al. 2008), which uses individual sample locations when estimating K and considers that genetically distinct clusters of individuals are likely to be separated spatially. We ran 10 repetitions for each value of K for K = 1–10 for the Lemhi Valley and for K = 1–6 in the Camas Prairie.
Sex-biased dispersal
We used several genetic analyses to determine if pygmy rabbits exhibited sex-biased natal dispersal. Analyses were conducted both among and within the CG, RC, and WS sample locations at the Lemhi Valley study area, where we sampled most extensively. To increase our power to detect genetic signatures of sex-biased natal dispersal at the CG and RC sample locations, we confined our analyses to juvenile rabbits sampled postdispersal in 2004 (after mid-October, or males > 380 g and females > 450 g) as recommended by Prugnolle and de Meeûs (2002). Because pygmy rabbits have low survival rates (Crawford 2008; Estes-Zumpf and Rachlow 2009; Price et al., in press; Sanchez 2007), we also included adults caught in early 2005, which were likely postdispersal juveniles from 2004. Because of low sample sizes at WS, however, we included both adults and postdispersal juveniles sampled in 2005 and 2006. At the study-area scale we tested for sex-specific differences in FST and average pairwise relatedness (r) computed separately for males and females among sample locations using the biased dispersal option in FSTAT. Because FSTAT calculates relatedness as r = 2FST/(1 + FIT), tests based on relatedness should perform similarly to tests based on FST (Goudet et al. 2002). Differences between sexes were tested with 10,000 randomizations following the randomization procedure of Goudet et al. (2002). The sex that disperses more frequently is expected to show less population substructure and be less closely related to neighboring rabbits; thus, the dispersing sex should have lower global FST and relatedness values than the more philopatric sex (Goudet et al. 2002; Mossman and Waser 1999; Prugnolle and de Meeûs 2002).
At the Camas Prairie study area sampling covered a broader spatial extent but was less extensive than sampling in the Lemhi Valley. Thus, we looked for evidence of sex-biased dispersal between sample locations using ancestry (q) values from program STRUCTURE to identify migrants. Individuals sampled in 1 cluster that had >70% ancestry from another cluster were considered migrants. We then compared the number of male and female migrants between clusters.
We also tested for evidence of sex-biased dispersal at a finer scale among individuals within the CG, RC, and WS sample locations. We compared sex-specific average pairwise relatedness (r) values, mean assignment index (mAI), and variance of assignment index (vAI) estimated under the biased dispersal option in FSTAT and tested with 10,000 randomizations. An assignment index can be used to determine the probability that a genotype originated in the population in which it was sampled (Paetkau et al. 1995; Waser and Strobeck 1998). If sex-biased dispersal occurs, the more philopatric sex should have a greater mAI than the dispersing sex (Favre et al. 1997; Prugnolle and de Meeûs 2002). Similarly, the variance around the assignment index is expected to be larger for the dispersing sex (Favre et al. 1997; Mossman and Waser 1999; Prugnolle and de Meeûs 2002). Because the more philopatric sex is expected to show stronger evidence of isolation by distance resulting from restricted movement and subsequent gene flow, we tested for sex-specific isolation by distance among individuals within sample locations. We tested whether the relationship between pairwise relatedness and Rousset's genetic distance between individuals (â—Rousset 2000) and geographic distance differed between sexes. We also examined multilocus spatial autocorrelation (Peakall et al. 2003; Smouse and Peakall 1999) among individuals of each sex from CG, RC, and WS using GENALEX version 6 (Peakall and Smouse 2006) because spatial genetic structure is influenced by dispersal behavior. Autocorrelation coefficients (r) range from 1 to −1 with 0 indicating that autocorrelation does not differ from random. Positive autocorrelation reflects nonrandom spatial genetic structure, the extent of which can be inferred from the 1st x-intercept where r = 0. We tested the significance of correlation coefficients against the null hypothesis of random structure using 1,000 permutations. We estimated 95% confidence intervals (95% CIs) around r-values for each sex by bootstrapping (1,000 bootstraps). We compared spatial genetic structure due to isolation by distance in each sex by examining spatial autocorrelation within discrete distance classes chosen to minimize variation in number of pairwise comparisons and size of 95% CIs. We compared the overall extent of nonrandom genetic structure (i.e., genetic “patch size”—Sokal and Wartenberg 1983) for each sex using sequential autocorrelation analyses at increasing distance classes.
Results
Genetic diversity
We obtained genotypes for 249 individuals from 8 sample locations (Table 1). We removed 7 juveniles from our analyses to avoid possible repeated samples from the same litter. We also removed 3 samples determined to be duplicates of previously sampled rabbits. Two loci showed significant deviations from Hardy–Weinberg equilibrium (P < 0.05 after Bonferroni adjustment). Sat05 deviated from Hardy-Weinberg equilibrium at 1 sample location in the Lemhi Valley (CG) and was excluded from all analyses within that study area. Similarly, A10 deviated from Hardy–Weinberg equilibrium at 1 of the sample locations in the Camas Prairie (JH) and was excluded from analyses within the Camas Prairie study area. We found no evidence of physical linkage among loci. Linkage disequilibrium detected among loci was not consistent across populations, thus all loci were treated as independent markers. Mean genotyping error rate for all microsatellite loci was 0.7%.
Table 1
Indices of genetic diversity for pygmy rabbits (Brachylagus idahoensis) sampled from 5 locations in the Lemhi Valley (SC = Smokey Cub, CG = Cedar Gulch, RC = Rocky Canyon, HC = Hawley Creek, WS = Warm Springs) and 3 locations in the Camas Prairie (MOR = Mormon Reservoir, JH = Johnson Hill, MAG = Magic Reservoir), Idaho, during 2002–2006 based on 15 nuclear microsatellite loci. Indices are averaged across all loci at each sample location. n = number of individuals sampled; AR = allelic richness (subsampled for the smallest number of individuals in each study: n = 8 for Lemhi Valley sites; n = 12 for Camas Prairie sites); HE = expected heterozygosity; HO = observed heterozygosity.
We documented moderate levels of heterozygosity within each study area. Average HE for the Lemhi Valley and Camas Prairie was 0.748 and 0.694 (Table 1). Allelic richness ranged from 4.8 to 5.6 in the Lemhi Valley and from 4.3 to 5.6 in the Camas Prairie (Table 1). Overall, the highest genetic diversity in the Lemhi Valley and the Camas Prairie study areas was at the HC and the JH sample locations. Using data from the mitochondrial control region, we identified 12 unique haplotypes across the Lemhi Valley and Camas Prairie (Table 2). We documented 9 haplotypes in the Lemhi Valley and 6 in the Camas Prairie. Haplotype diversity in the Lemhi Valley and Camas Prairie was greatest at the RC and JH sample locations (Table 2).
Table 2
Haplotype frequencies and diversity (h) for the mitochondrial DNA control region for pygmy rabbits (Brachylagus idahoensis) sampled from 3 locations in the Lemhi Valley (CG = Cedar Gulch, RC = Rocky Canyon, WS = Warm Springs) and 3 locations in the Camas Prairie (MAG = Magic Reservoir, JH = Johnson Hill, MOR = Mormon Reservoir), Idaho, during 2002–2006.
Genetic structure
We documented relatively low levels of genetic structure based on nuclear microsatellite data among our 5 sample locations in the Lemhi Valley (n = 201 individuals). FST-values ranged from 0.006 to 0.054. All but 2 comparisons between sample locations deviated significantly from panmixia (FST = 0); however, most FST-values indicated only low levels of genetic differentiation (Table 3a). Number of loci showing significant genetic differentiation ranged from 1 to 9 and was greatest between CG and WS. We found no evidence of isolation by distance using either FST (Mantel coefficient [rM] = 0.042, P = 0.208) or DS (rM = 0.058, P = 0.208). Distribution of mtDNA haplotypes (Table 2) and an exact test for population differentiation based on haplotype frequencies (Table 4a) among CG, RC, and WS suggested some restriction in gene flow between WS and locations on the east side of the valley, but no restriction between CG and RC. FST-values based on mtDNA control region sequences indicated a similar trend (Table 4a).
Table 3
Population differentiation among sample locations of pygmy rabbits (Brachylagus idahoensis) in the a) Lemhi Valley (n = 201 individuals; SC = Smokey Cub, CG = Cedar Gulch, RC = Rocky Canyon, HC = Hawley Creek, WS = Warm Springs) and b) Camas Prairie (n = 38 individuals; MAG = Magic Reservoir, JH = Johnson Hill, MOR = Mormon Reservoir), Idaho, based on 15 microsatellite loci. Below diagonal: pairwise FST-values. Significant deviation from 0 (panmixia) at P < 0.05 (*) and P < 0.001 (**) was assessed using 10,000 randomizations in Arlequin 2.0. Above diagonal: number of loci (of 15) showing significant genetic differentiation based on 1,000 randomizations in GENEPOP 3.1.
Table 4
Genetic diversity and differentiation of pygmy rabbits (Brachylagus idahoensis) based on mtDNA control region sequences from a) 3 sample locations in the Lemhi Valley (n = 75 individuals; CG = Cedar Gulch, RC = Rocky Canyon, WS = Warm Springs) and b) 3 sample locations in the Camas Prairie (n = 33 individuals; MAG = Magic Reservoir, JH = Johnson Hill, MOR = Mormon Reservoir), Idaho. Below diagonal: pairwise FST-values. Significant deviation from 0 (panmixia) at P < 0.05 (*) and P < 0.001 (**) was assessed based on 10,000 randomizations in Arlequin 2.0. Above diagonal: P-values for an exact test of population differentiation with 10,000 randomizations in Arlequin 2.0.
Results from all 3 Bayesian clustering programs suggested that the 5 sample locations in the Lemhi Valley comprise 1 genetic cluster with sufficient gene flow to prevent substantial genetic differentiation among all locations. Using the method by Pritchard et al. (2000), the optimum K-value revealed by STRUCTURE was K = 1 (Fig. 2). ΔK (Evanno et al. 2005) cannot be calculated when the highest likelihood is at K = 1, and, therefore, the Evanno et al. (2005) approach was not applied to this study area. Programs BAPS (spatial option) and STRUCTURAMA also identified K = 1 as the optimal number of genetic clusters with a posterior probability of 1.0 and 0.89.
Fig. 2
Results from program STRUCTURE analysis of pygmy rabbits (Brachylagus idahoensis; n = 201) from 5 sample locations in the Lemhi Valley, Idaho. Plot displays mean log-likelihood values (LnP(D)) for 10 independent runs for each value of K for K = 1–10. The highest value was at K = 1, indicating that the 5 sample locations likely form 1 population.
In contrast to the Lemhi Valley, we documented substantial evidence of population substructure among rabbits (n = 38) in the Camas Prairie. FST-values ranged from 0.059 to 0.138, and all were significantly different from 0 (Table 3b). Number of loci showing significant genetic differentiation ranged from 8 to 12 and was greatest between MAG and JH. Distribution of mtDNA control region haplotypes (Table 2) and an exact test for population differentiation based on mtDNA (Table 4b) indicated a restriction in gene flow among all pairs of populations. FST-values based on control region sequences indicated potential restriction in gene flow between MAG and the other locations but no restriction between JH and MOR (Table 4b). However, differentiation of MAG from other sample locations might be influenced by the lower genetic diversity documented there (Tables 1 and 2). Program STRUCTURE revealed 2 genetic clusters (Fig. 3). Results from BAPS (spatial) revealed a 0.94 posterior probability of 3 genetic clusters. However, one cluster contained only 4 individuals from JH. All other samples from JH and MOR comprised the 2nd genetic cluster, whereas all individuals from MAG comprised the 3rd. STRUCTURAMA indicated a 0.58 probability of 2 clusters and a 0.34 probability of 3 clusters. When we set K = 2, rabbits from JH and MOR again formed 1 cluster, whereas all MAG rabbits formed a separate cluster. When we used K = 3, the first 2 groups each contained rabbits from both JH and MOR, whereas the 3rd cluster contained only rabbits from MAG.
Fig. 3
Results from program STRUCTURE analysis of pygmy rabbits (Brachylagus idahoensis; n = 38) from 3 sample locations in the Camas Prairie, Idaho. Plot displays mean log-likelihood (LnP(D)) and ΔK-values for 10 independent runs of each value of K for K = 1–6, indicating that the 3 sample locations likely comprise 2 populations.
Sex-biased dispersal
Although pairwise FST and relatedness values were slightly higher for females (Table 5), we found no significant evidence of sex-biased dispersal among sample locations in the Lemhi Valley. Similarly, no evidence of a sex bias in dispersal was detected when we examined ancestry values from a STRUCTURE analysis assuming 2 genetic clusters for rabbits in the Camas Prairie. All MOR and JH rabbits were assigned to the same cluster and had ancestry values greater than 93%. All rabbits sampled in MAG were assigned to the 2nd cluster, with only 1 juvenile male (230 g) having an ancestry value below 70% (68%).
Table 5
Sex-specific indices of genetic differentiation (FST) and relatedness (r) among 3 locations sampled for pygmy rabbits (Brachylagus idahoensis) in the Lemhi Valley, Idaho (CG = Cedar Gulch, RC = Rocky Canyon, WS = Warm Springs). Samples were analyzed for 48 rabbits from CG (males [M] = 24, females [F] = 24), 43 rabbits from RC (M = 18, F = 25), and 27 rabbits from WS (M = 14, F = 13). P-values were based on 10,000 randomizations performed in FSTAT.
Analyses at the individual level revealed that although females were assigned to their sample location more often than males, the mAI values for females (0.008) and males (−0.009) were not significantly different (P = 0.981). Similarly, the variance in the corrected assignment index was slightly higher for males (13.35) than females (12.31); however, the difference was not significant (P = 0.783). We also found no evidence of isolation by distance among individuals of either sex within study areas. Neither pairwise relatedness nor Rousset's â decreased significantly with geographic distance for either sex (Table 6), except for pairwise relatedness among males in CG, but these results were inconclusive (P = 0.055). Spatial autocorrelation analyses at discrete distance classes revealed fine-scale structuring at shorter distances for males (<1,100 m), and females showed a trend for lower levels of genetic structure extended over longer distances (<6,500 m; Fig. 4). However, spatial genetic structure in both sexes did not differ significantly from random except among males separated by 400–800 m. The overall extent of genetic correlation for increasing distance classes also was similar between sexes. Positive genetic structure extended up to 6–7 km for both males and females (Fig. 5a and 5b).
Fig. 4
Mean spatial genetic correlation (±95% CI) at discrete distance classes for male (M; n = 56) and female (F; n = 62) pygmy rabbits (Brachylagus idahoensis) across 3 sample locations in the Lemhi Valley, Idaho. Dashed lines representing upper and lower 95% CLs for a random distribution of male (black) and female (gray) pygmy rabbits were based on 1,000 permutations. Number of pairwise comparisons for males and females within each distance class is given at the top of the graph.
Fig. 5
Mean (±95% CI) spatial genetic correlation, r, for sequential autocorrelation analyses at increasing distance classes for a) male (n = 56) and b) female (n = 62) pygmy rabbits (Brachylagus idahoensis) in the Lemhi Valley, Idaho. Upper (U) and lower (L) 95% CIs for a random distribution were based on 1,000 permutations.
Table 6
Results from sex-specific isolation-by-distance analyses of relatedness (r) and Rousset's genetic distance (â) against geographic distance among individuals within 3 sample locations in the Lemhi Valley, Idaho (CG = Cedar Gulch, RC = Rocky Canyon, WS = Warm Springs). Mantel coefficients (rM) and P-values were generated in program zt (Bonnet and Van de Peer 2002).
Discussion
Levels of genetic diversity in pygmy rabbits in our study areas were relatively high. Examination of our data suggested that gene flow is occurring, or recently has occurred, among sample locations or between sample locations and neighboring areas occupied by rabbits. Heterozygosity documented in this study was comparable to values identified in the lagomorphs Lepus americanus (Burton and Krebs 2003) and O. cuniculus (Richardson et al. 2002; Webb et al. 1995). The only other study of genetic diversity in pygmy rabbits, largely based on museum specimens, documented similar levels of heterozygosity in rabbits from Idaho and Montana, but significantly lower levels of heterozygosity in the endangered Columbia Basin pygmy rabbit in Washington (Warheit 2001). We added species-specific loci (Estes-Zumpf et al. 2008) to those used by Warheit (2001); however, ascertainment bias detected in other mammalian species (Garner et al. 2005) was not apparent in our study.
Contrary to expectations, we did not document strong evidence of genetic substructure, based on nuclear microsatellites, among populations of pygmy rabbits within our study areas. Sample locations within approximately 13 km of each other in each study area exhibited sufficient gene flow to constitute single populations. This distance is close to maximum natal dispersal distances documented for pygmy rabbits (Estes-Zumpf and Rachlow 2009). The only location separated from others by a larger distance (approximately 30 km) at the Camas Prairie study area was identified as a separate genetic cluster. Although levels of genetic differentiation among populations differed significantly from panmixia, actual values were low to moderate (Tables 3a and 3b).
Lack of strong population substructure within study areas indicated that perennial streams and roads did not act as substantial barriers to gene flow. Genetic methods have been used to detect significant impacts of roads on species with both longer (Ursus arctos [Proctor et al. 2005] and Canis latrans and Lynx rufus [Riley et al. 2006]) and shorter (Myodes glareolus—Gerlach and Musolf 2000) generation times than pygmy rabbits. Highway 28 near Leadore in the Lemhi Valley was 1st paved in approximately 1940, and significant restriction in gene flow caused by this highway and other roads would likely have been detected if it occurred. We documented gene flow across rural 2-lane highways at both study areas and across several perennial creeks in the Lemhi study area. Furthermore, radiotelemetered animals in the Lemhi Valley were documented crossing creeks (Estes-Zumpf and Rachlow 2009; Sanchez 2007). One interesting exception was the SC sampling location, which was separated from all other locations in the Lemhi study area by both a highway and a creek. Although SC was <1 km from the closest sample location (CG), it was more differentiated from other locations and exhibited the lowest genetic diversity. The SC sampling location was relatively small and peripheral to the larger expanse of occupied habitat in the Lemhi study area, which likely affected patterns of genetic variation. Although the highway and creek also might limit dispersal between SC and other sites, they do not cause complete isolation. For example, 1 rabbit originally captured on CG was later recaptured on SC during a concurrent mark–recapture study (Estes-Zumpf 2008).
In contrast to roads and streams, agricultural expanses appeared to act as filters if not complete barriers to gene flow. Although the WS sampling location on the west side of the Lemhi Valley was roughly equidistant from all sample locations on the east side, FST-values based on nuclear microsatellites indicated that the amount of population substructure between WS and other locations appeared to increase as the amount of agriculture between locations increased (Table 3a; Fig. 1a). This pattern is consistent with either a lower rate of successful dispersal through agriculture or dispersal around agricultural fields through surrounding sagebrush habitat. Although we did not have DNA from historic rabbit samples that predate road and agricultural development in the valley for comparison, the valley historically contained relatively continuous sagebrush cover. Thus, observed patterns in genetic structure, when present, are likely due to recent landscape changes rather than historic patterns in landcover. These results suggest that agriculture might increase isolation among populations of pygmy rabbits and that more work is needed to quantify how different land uses affect permeability of the landscape for this species.
We detected little evidence of sex-biased natal dispersal in pygmy rabbits either within or among sample locations at different spatial extents. Although different indices perform better under different scenarios (Goudet et al. 2002; Mossman and Waser 1999; Prugnolle and de Meeûs 2002), Goudet et al. (2002) determined that FST had the most power to detect sex-biased dispersal under the widest range of sampling schemes and dispersal rates. Although FST, mAI, and vAI indices revealed a consistent but slight trend toward male-biased dispersal, differences between sexes were not significant. Similarly, a radiotelemetry study of natal dispersal in pygmy rabbits documented no significant difference in dispersal rates between sexes, although the trends were consistent with this study (dispersal rate for males = 90%, for females = 80%—Estes-Zumpf and Rachlow 2009).
Despite the lack of a sex bias in dispersal frequency, radiotelemetry of juvenile rabbits detected evidence of a female bias in dispersal distance, with females moving 3 times farther than males, on average (Estes-Zumpf and Rachlow 2009). This bias in dispersal distance has potential consequences for population genetic structure. However, once dispersing juveniles settle and breed, genetic signature associated with natal dispersal is likely to be lost or diluted (Prugnolle and de Meeûs 2002). Smouse and Peakall (1999) suggested that multivariate spatial autocorrelation analyses might be more sensitive to discrepancies in genetic structure between sexes than other indices. We used 2 spatial autocorrelation analyses that explore differing aspects of population genetic structure. One analysis examined isolation by distance among individuals across discrete distance classes (Smouse and Peakall 1999) and should reflect patterns of natal dispersal. The 2nd analysis examined the overall genetic patch size (Sokal and Wartenberg 1983), which reflects a combination of effective natal dispersal (juveniles that survived and reproduced) and the extent of genetic exchange among breeding adults.
Inferences about dispersal behavior of pygmy rabbits based on the 2 types of spatial autocorrelation analyses that we used suggested that patterns of sex-specific genetic structure in this species are complex, likely exhibiting both temporal and spatial variability. Results from the 1st analysis (Fig. 4) were consistent with a female bias in dispersal distance because correlation between genetic distance and geographic distance was greater among males than females. However, this trend was not apparent in the 2nd analysis over increasing distance classes (Fig. 5), likely because of extensive movement by breeding adults, especially by males (Burak 2006; Crawford 2008; Sanchez and Rachlow 2008), and high mortality rates of juvenile rabbits (69–89%—Estes-Zumpf and Rachlow 2009; Price et al., in press). Differences in frequencies of maternally inherited mtDNA haplotypes among most sample locations suggested a possible restriction in female gene flow at distances ≥ 7 km. Furthermore, FST and relatedness values among sample locations tended to be higher for females, suggesting greater clustering of genetically similar individuals in females than males. Because genetic patch size was similar among postdispersal juveniles of both sexes, and sample locations within approximately 13 km of each other were considered the same genetic cluster in this study, gene flow among sample locations might be mediated largely by movements of breeding males.
This research provides the 1st data on fine-scale population genetic structure for this mammal of conservation concern. Gene flow occurred across greater spatial extents than previously expected, and although roads in this study did not present significant barriers to gene flow, they were characterized by low traffic volume and a lack of human development. It is unknown how increased traffic volume and associated changes in land use might affect movement by this species. Other studies have documented differences in response of animals to roads with high and low traffic volumes (Clevenger et al. 2003; Gagnon et al. 2007). Vegetation structure potentially affects movement patterns of pygmy rabbits (Sanchez and Rachlow 2008), and therefore agricultural practices also might affect dispersal and gene flow in this species. Examination of our data suggested reduced gene flow across agricultural expanses, but more information is needed to assess the permeability of different landcover types. An understanding of how landscape features influence movement and gene flow in pygmy rabbits is important for predicting and mitigating the effects of ongoing land-use change in the sagebrush-steppe ecosystem.
Acknowledgments
We thank C. Scofield, J. Witham, D. Sanchez, V. Guyer, B. Waterbury, and H. Roberts for their logistical support and assistance in the field. C. Anderson, K. Gebhardt, S. Robinson, J. Adams, N. Balkenhol, and the entire Waits Lab provided considerable help with genetic analyses. We also thank the Washington Department of Fish and Wildlife, especially Cheryl Dean, for assistance with mtDNA analyses. Comments from several colleagues and an anonymous reviewer improved the manuscript. We are grateful to the United States Forest Service Leadore Ranger District, especially D. Hilliard, M. Steck, and R. Riebe, for providing logistical support and housing while in the field. We sincerely thank the Bureau of Land Management, Idaho Department of Fish and Game, Montana Department of Fish, Wildlife and Parks, the United States Forest Service, the United States Fish and Wildlife Service, and the University of Idaho for their financial and logistical support of this project. This publication is Contribution 1042 of the Idaho Forest, Wildlife, and Range Experiment Station, University of Idaho, Moscow.
