Comparative Genetic Structure of Sympatric Leporids in Southern Illinois

Many leporid species have overlapping ranges, but ecological differences may make for profoundly different population structure in the same area. In southern Illinois, swamp rabbits (Sylvilagus aquaticus) and eastern cottontails (S. floridanus) co-exist, but swamp rabbits are habitat specialists associated with floodplain forests and eastern cottontails exhibit more general habitat preferences. Additionally, swamp rabbits exist at the northern edge of their range in southern Illinois, whereas eastern cottontails are well within the core of their range. To determine whether these differences resulted in differences in genetic structure, we assessed the population differentiation of these 2 sympatric species in southern Illinois using polymorphic microsatellite markers. We employed a combination of sampling techniques including tissue sampling from trapped rabbits and collection of fecal pellets from swamp rabbit latrine logs. Swamp rabbits appeared to be differentiated within 50 km of Illinois’ Cache River watershed (K = 4 populations) suggesting that local populations are relatively isolated, whereas eastern cottontails exhibited no discernable genetic structure (K = 1). Such a result confirms the expected relationship of greater genetic structure in a habitat specialist (i.e., swamp rabbits) versus a habitat generalist (i.e., eastern cottontails). Unlike eastern cottontails, our results suggest that swamp rabbits experience low genetic connectivity in southern Illinois likely due to restricted dispersal, linear distribution of habitat, and smaller effective population size as mediated by habitat fragmentation.

Landscape features shape the demographic and genetic connectivity of wildlife populations (Holderegger and Wagner 2008;Storfer et al. 2010) and predicting population connectivity is of central importance to species conservation and landscape planning (Hilty et al. 2006). Barriers and expanses of highly inhospitable habitat such as freeways and mountain ranges may impact co-distributed species similarly by limiting gene flow (Riley et al. 2006;Crawford et al. 2007;Delaney et al. 2010). However, the population genetic structure of any 1 species in a given area can also depend on finer-scale responses to the landscape dictated by the ecological and behavioral characteristics of that species (Crawford et al. 2007). In other words, similar co-distributed species may have different patterns of connectivity in a common landscape due to slight but important differences in their natural history (Kraaijeveld-Smit et al. 2007;DiLeo et al. 2010).
In addition to those landscape factors influencing species movement, demographic effects dictating effective population size (N e ) can also influence genetic structure and diversity (Schwartz et al. 2003;Eckert et al. 2008). At species' range limits, increased variability in the persistence of populations generates the expectation that species are less abundant and occur in a patchier arrangement (Andrewartha and Birch 1954;Sagarin and Gaines 2002;Sexton et al. 2009). In turn, the bottlenecks and range expansions of populations at species' margins can lead to smaller N e , greater genetic drift, lower diversity, and greater genetic structure (Nei et al. 1975;Eckert et al. 2008). Conversely, the same geographic area may host species at the center of their range, where a larger N e may result in less genetic structure and greater diversity than for species at their range limits.
The diversity of ecological factors affecting genetic structure makes population connectivity difficult to infer using knowledge gained from species with superficial ecological similarities or co-distributed, related species. Nevertheless, to predict population viability, evolutionary potential, and a host of biological outcomes, population connectivity and the potential for genetic differentiation are often inferred based on whether a species may be classified as a specialist or generalist (Jonsen and Fahrig 1997;Swihart et al. 2003). Comparisons of sympatric mammals show that specialist species with narrower habitat requirements can exhibit greater genetic structure than generalist species (Kraaijeveld-Smit et al. 2007;Devine 2012). This can be due to a number of factors that result in decreased population connectivity such as larger distances between habitat patches and restricted dispersal corridors (Fahrig and Merriam 1994;Rosenberg et al. 1997;Hanski 1998;Swihart et al. 2003). However, co-distributed leporids have not always conformed to this generalization (Thimmayya and Buskirk 2012).
Swamp rabbits (Sylvilagus aquaticus) and eastern cottontails (S. floridanus) are leporid species that co-occur as legally harvested game species in southern Illinois. Considered habitat specialists, swamp rabbits occur in bottomland hardwood forest habitats and adjacent recently afforested areas (Scharine et al. 2011) in most of the southeastern United States (Chapman and Feldhamer 1981;Allen 1985) and their occupancy is influenced by distance to wetlands (Terrel 1972;Scharine et al. 2011). Habitat at the northernmost portion of swamp rabbit range (e.g., in southern Illinois and southern Indiana) is patchily distributed (Kjolhaug et al. 1987;Barbour et al. 2001;Roy Nielsen et al. 2008). A species with similar habitat preferences to the swamp rabbit (Allen 1985), the marsh rabbit (S. palustris), was reluctant to travel through open areas (Forys and Humphrey 1996). If similar, swamp rabbits may be effectively isolated in remnant bottomland hardwood forest due to a disinclination to cross inhospitable cover types. Such isolation could result in genetic differentiation among habitat patches due to drift over time (Frankham et al. 2002) and, thus, we predicted that swamp rabbits exhibit significant genetic structure among subpopulations.
Eastern cottontails, in contrast, are the most widely distributed member of the genus ranging from Canada through Mexico (Chapman and Litvaitis 2003). Early successional habitats and other forms that provide shelter and resting cover are important for eastern cottontails (Chapman et al. 1980;Litvaitis 2001). Their preferred cover depends on location and availability, and thus selected habitat varies among regions and seasonally (Chapman et al. 1980). In Illinois, eastern cottontails exhibit more general habitat associations than swamp rabbits (Scharine et al. 2011). The eastern cottontail is less reliant on dense shrub cover than New England cottontails (S. transitionalis) and brush rabbits (S. bachmani) potentially because of adaptations that allow more efficient predator detection in open habitats (Smith and Litvaitis 2000;Chapman and Litvaitis 2003). Swamp rabbits are not more closely associated with dense shrub cover than eastern cottontails in Illinois, possibly due to use of open water for predator avoidance (Scharine et al. 2011).
An understanding of functional landscape connectivity can provide information crucial to maintaining the viability of both rabbit populations and mammalian biodiversity in southern Illinois. Many studies have successfully applied molecular methods to illuminate genetic structure in North American rabbits (Scribner and Chesser 1993;Crouse et al. 2009;Estes-Zumpf et al. 2010;Thimmayya and Buskirk 2012;Fenderson et al. 2014). However, no studies have addressed the genetic structure of eastern cottontails and swamp rabbits with highly polymorphic microsatellites. In southern Illinois, genetic structure and diversity of these species may differ because of their ecological requirements (i.e., generalist versus specialist), the distribution of their habitats, and demographic factors related to location within their range (i.e., core versus periphery). We hypothesized that swamp rabbits would have greater genetic structure and lower genetic diversity than eastern cottontails in southern Illinois. To determine the degree of population genetic connectivity and genetic diversity of 2 ecologically similar rabbit species in a common landscape, we conducted a population genetic analysis of microsatellite markers.

Materials and Methods
We obtained genetic samples from swamp rabbits and eastern cottontails in the bottomland hardwood forests of the Cache River and adjacent watersheds in southern Illinois ( Fig. 1; approximately 37°N, 89°W). We also collected fecal pellets of swamp rabbits from latrine logs (Zollner et al. 2000;Schauber et al. 2008). Sites were selected based on their potential to be occupied by swamp rabbits and were identified using >20 years of research conducted in this area (Kjolhaug et al. 1987;Porath 1997;Woolf and Barbour 2002;Rubert 2007;Scharine et al. 2009). Bottomlands contained slight variations in relief that produced ridges, sloughs, and swamps. These features varied considerably in soil composition, drainage conditions, and tree species (Hosner and Minkler 1963). Agricultural development in the 1970s significantly altered land cover and hydrology of the Cache River watershed and extensive restoration to hardwood tree species has since occurred (Kruse and Groninger 2003 (Scharine et al. 2011). Traps were baited with a quartered apple and covered with burlap and brush to provide greater security for captured rabbits. Trapping effort was focused on reforested sites (Scharine et al. 2011) and areas where both species occurred (Crawford 2014). Swamp rabbit ears were also received from hunter harvests during 2007 and 2010. Animals were captured and handled in accordance with guidelines of the American Society of Mammalogists (Sikes et al. 2011), and methods were approved by the Institutional Animal Care and Use Committee at Southern Illinois University Carbondale (protocol 06-035). A small section of tissue (1 mm) was excised from the outer edge of the ear of captured animals (Halanych and Robinson 1997) using a biopsy punch. This piece of tissue was stored in an individual 2-ml screw-cap or 15-ml centrifuge tube and frozen at −20°C within approximately 2 h of collection. For longer term storage (> 1 month) tissue was frozen at −80°C. Swamp rabbit tissue samples were obtained from unique trapped animals (n = 114), hunter harvests (n = 13), and a roadkill (n = 1) directly adjacent to a trapping site (Fig. 1). One pair of trapped individuals (the 1st trapped in February 2007, the 2nd in January 2010) from the same site produced a matching genotype; only 1 of the genotypes was included in further analyses. Eastern cottontail samples were obtained from trapped individuals (n = 79) and road-killed individuals adjacent to trapping sites (n = 2).
Fecal pellet collection.-We collected pellets from swamp rabbit latrine sites only. Swamp rabbits use natural (Zollner et al. 2000) and artificial latrine logs (Schauber et al. 2008) for repeated defecation. Artificial latrine logs were placed on the study area during November-December 2006 (Schauber et al. 2008). These logs were revisited during November-March 2009-2010. Additionally, natural latrine logs were located from November-March 2009-2010 outside of trapped areas and in areas with artificial latrine logs. We searched all bottomland forested habitat at each site ) and collected samples from all latrines at all occupied sites. Both natural and artificial latrine logs were cleared of all visible fecal pellets and marked with flagging. Latrines were revisited ≤ 2 weeks later and new pellets were collected, placed on ice during transport, and stored at −80°C until processed. Pellets ≤ 10 cm apart were considered the same sample.
Tissue genotyping.-We extracted rabbit DNA with the Wizard Genomic DNA Purification Kit (Promega, Madison, Wisconsin) according to their Mouse Tail Protocol. We screened 23 microsatellite primer pairs published for the European rabbit (Oryctolagus cuniculus) and the pygmy rabbit (Brachylagus idahoensis- Rico et al. 1994;Mougel et al. 1997;Surridge et al. 1997;Estes-Zumpf et al. 2008) and 23 markers developed in conjunction with this study (Berkman et al. 2009) for 10 swamp rabbit tissue samples. We identified 16 polymorphic microsatellite loci that could be amplified and scored consistently in swamp rabbits and 15 loci in eastern cottontails. Reactions included 10 µl of a PCR mixture containing ThermoScientific (Grand Island, New York) PCR Master Mix, 0.2 µM of each primer (unlabeled and dye-labeled), and approximately 20 ng of genomic DNA. Reactions were subjected to 35 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 30 s for all primers except sat8, A133, D121, sfl013, and A121, which were subjected to an annealing temperature of 54°C. Multiplex combinations included sat16 with D118, sat7 with sfl011, A133 with D121, and A121 with sfl013. All other loci were amplified singly. Products were visualized with Rox-400 size standard (Gel Company, San Francisco, California) on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Grand Island, New York).
To detect scoring errors and null alleles in groups in which no substructure was expected (Chapuis and Estoup 2007), we analyzed 2 groups of rabbit genotypes from tissues sampled in very close proximity ( Fig. 1) to test for conformation to Sandusky Hickory Horseshoe Lake ) for the 2 species (light gray = eastern cottontail; dark gray = swamp rabbit) in the western hemisphere (right inset) and in Illinois (left inset). Illinois is outlined in the right inset and the extent of the study area map is given in a dark black box on the left inset. Location names are given in regular text and watercourse names are given in italicized text. Sample numbers are indicated on the map beside sample locations marked as diamonds (swamp rabbits; regular text number) or squares (eastern cottontails; italicized number). The designation "P" indicates that the genotype originated from a fecal pellet; all other genotypes were obtained from tissue samples. Light gray lines are watercourses.
Hardy-Weinberg equilibrium (HWE-Guo and Thompson 1992). For swamp rabbits, we analyzed the Horseshoe Lake group (n = 16) and the Hickory group (n = 41; Fig. 1). For eastern cottontails, we analyzed the Sandusky group (n = 21) and the Hickory group (n = 26). We calculated the probability of deviation from HWE across groups and tested for linkage disequilibrium for all pairs of loci within groups using a Bonferroni correction of α = 0.05 (Rice 1989).
Fecal pellet genotyping.-All pre-PCR procedures were conducted on a bench dedicated to noninvasive samples only, with dedicated equipment and materials. We used aerosol-resistant pipet tips and carried negative controls through the process from extraction to fragment analysis. We extracted DNA from fecal pellets with the Zymo Research Fecal DNA Mini Prep Kit (Zymo Research, Irvine, California). We cut each pellet into a minimum of 8 approximately equal-sized pieces and divided the material into 2 tubes for extraction. We used the manufacturer's protocol but modified the volume of lysis solution to 900 µl to compensate for the absorbency of the fecal material. Genomic DNA from the same pellet was combined at the end of the process to provide a sufficient template volume and concentration. An additional pellet from the same sample was extracted when the initial genomic DNA extraction ran out. Analysis of additional pellets from the same sample always produced concordant results.
PCR conditions consisted of a 40-µl PCR mixture of 0.4 µM of each primer, 0.2 mg/ml of bovine serum albumin, 1.25 units of AmpliTaq Gold DNA Polymerase, 10× PCR Gold Buffer (150 mM Tris-HCl, 500 mM KCl), 2.0 mM MgCl 2 , 200 µM each dNTP, and 4 µl of genomic DNA extract. Cycling conditions and fragment visualization were identical to those described above except that 45 cycles were performed. Eight tissue samples were genotyped with the fecal pellet protocol in addition to the tissue genotyping protocol to ensure consistency of allele scores.
Fecal samples (n = 86) from swamp rabbits were genotyped at a subset of reliable loci, specifically sfl006, sfl014, and sat16. Samples that produced a product for the subset of loci (n = 46) were genotyped for the full set of 16 loci. Samples that subsequently failed for the majority of PCR reactions were discarded. The remaining samples (n = 21) produced genotypes at ≥7 loci. To obtain accurate genotypes from fecal pellet samples, we employed the comparative approach of Frantz et al. (2003). This approach is suited to situations in which multiple pellets are expected to originate from 1 individual (Frantz et al. 2003) and can be less costly and less error prone (Hansen et al. 2008) than the more commonly used multiple tubes approach of Taberlet et al. (1996). Consensus genotypes were determined from a minimum of 2 PCRs for every locus. If a heterozygous genotype was obtained in both PCRs, we considered the consensus score achieved. If only 1 allele was observed in the first 2 PCRs or 2 different alleles were observed singly in each PCR, a 3rd PCR score was obtained. If only 1 allele was observed in all 3 PCRs, the individual was scored as a homozygote at the locus. If another allele was observed, up to 5 additional PCRs were performed until each allele was observed twice. The maximum number of PCRs was less than the 7 recommended by Frantz et al. (2003) due to resource constraints and the large number of loci examined. After performing pairwise comparisons of genotypes, 1 pair remained unresolved at 1 locus for which 1 individual produced only 1 allele in 2 PCRs. An additional 4 amplifications were performed and the 2nd allele was observed twice resolving the pairwise difference.
Consensus scores for samples from the same site were compared with one another. If the consensus genotypes matched, they were considered to have originated from the same individual and grouped. If ≥3 alleles at any locus were observed in a pairwise comparison, the 2 genotypes were considered to have originated from different individuals. All data from grouped consensus genotypes were pooled to produce a final genotype for analysis. Consequently, a missing genotype at a locus from 1 sample was replaced if a matching sample contained a consensus genotype at that locus. Error rates for pellets were estimated by comparing consensus genotypes to the PCR replicates used to obtain them.
The probability of identity, the probability that 2 siblings shared the same genotype (P IDSIB ), was calculated with the program GenAlEx v. 6.2 (Peakall and Smouse 2006). We calculated P IDSIB at each locus and cumulatively for all loci. Not all pellet samples produced genotypes for the full suite of loci, so for these individuals, P IDSIB was calculated for the available loci.
Genetic diversity and N e .-We assessed differences in diversity among eastern cottontails and swamp rabbits by comparing the allelic richness of loci that were typed for both species (n = 10). Allelic richness was calculated using HPRare (Kalinowski 2005) and was standardized to the minimum number of scored alleles (n = 134) in the sample. We compared jackknife estimates of the mean allelic richness and the 95% confidence intervals (CI). We calculated the inbreeding effective population size with the linkage disequilibrium method implemented in NeEstimator v. 2 (Do et al. 2014). We used only alleles with frequency > 0.05, which resulted in 1,245 independent comparisons for swamp rabbits and 1,152 for eastern cottontails. We compared estimates of N e and their 95% jackknife confidence intervals.
Genetic differentiation.-We used Bayesian clustering approaches to examine evidence for distinct subpopulations of swamp rabbits and eastern cottontails. We employed alternative analyses because different statistical models make different assumptions that may or may not be violated by the unknown demographic history of the population under study (Francois and Durand 2010). Two Bayesian clustering methods were used: 1 aspatial model (STRUCTURE) and 1 spatial model (TESS). We chose approaches that offered admixture models (Falush et al. 2003;Chen et al. 2007) because these models performed well with simulations of weakly diverged parental populations admixed along a spatial gradient (Francois and Durand 2010). Such a dynamic was expected given our relatively continuous sampling distribution (Fig. 1). STRUCTURE (Pritchard et al. 2000) attempts to find groups that minimize HWE and linkage disequilibrium and does not assume that all source populations have been thoroughly sampled. Ten STRUCTURE simulations with K ranging from 1 to 15 were performed using 300,000 Markov chain Monte Carlo steps with 300,000 burnin steps assuming correlated allele frequencies. We examined the 1st modal peak of the likelihood of the data (lnP(D)) to determine the best choice of K (Pritchard et al. 2000). TESS takes into account location of individuals when determining the most likely number of interbreeding groups and jointly estimates admixture proportions and optimal number of clusters by allowing admixture proportions to vary in space. The TESS algorithm was run with 10,000 burn-in sweeps and a running period of 50,000 sweeps. We ran the program 10 times for each K MAX between 2 and 15 with the CAR model of admixture (Durand et al. 2009b). To determine the most likely number of interbreeding groups with TESS, we averaged the deviance information criterion (DIC) for replicate runs and examined plots of the average DIC for increasing K MAX to determine when DIC decreased and reached a plateau (Durand et al. 2009a). Additionally, we examined genetic admixture plots for individual rabbits and looked for stabilization, which is the point where membership in a particular genetic cluster does not change with an increase in K MAX (Durand et al. 2009a). The resulting genetic admixture values for each individual resulting from multiple Bayesian simulations were summarized with CLUMPP v. 1.1.2 (Jakobsson and Rosenberg 2007).
Using the clusters informed by the aspatial and spatial Bayesian analyses, we grouped swamp rabbit samples into subpopulations. Locations (Fig. 1) were assigned to a swamp rabbit subpopulation based on the average proportion of membership in a particular genetic cluster (Q) of the individual genotypes from that location. Average Q-values from the 2 Bayesian clustering methods and an initial run of BAPS (Corander et al. 2008) were examined. For swamp rabbit subpopulations designated using the clustering results, we calculated global and pairwise F-statistics (Weir and Cockerham 1984) with GenAlEx (Peakall and Smouse 2006) and their significance was determined with 999 bootstrap replicates. One individual from the Ballard location in Kentucky ( Fig. 1) was excluded from this analysis because we could not confidently assign it to a subpopulation based on location. We divided eastern cottontail samples into subpopulations that matched the subpopulation boundaries that we designated for swamp rabbits in order to compare random differentiation among groups as measured by F ST . These groups were not intended to reflect true subpopulations of eastern cottontails. Only 4 eastern cottontails were in the southwest subpopulation (Fig. 1), so we grouped them with the south-central subpopulation and calculated global and pairwise F-statistics similarly for only 3 groups of eastern cottontails.
Spatial autocorrelation.-We examined the decrease in genetic similarity over distance for both species to assess finescale genetic structure and differences in dispersal patterns. The calculation of the spatial autocorrelation coefficient (r) was performed among pairs of individuals using GenAlEx v. 6.2 (Peakall and Smouse 2006) and 95% CI were determined with 999 bootstrap replicates. We chose 1 km as the lag distance because this is the maximum distance that most individual cottontails move (Shields 1960;Chapman and Trethewey 1972;Scribner and Warren 1990).

Results
The number of alleles per locus for swamp rabbits ranged from 2-10 and averaged 5.1. For eastern cottontails, 1 locus could not be amplified (sol28). The remaining loci ranged from 2-16 and averaged 9.1 alleles per locus. Locus sfl008 approached significance for deviation from HWE in the Horseshoe Lake group of swamp rabbits and did not contain sufficient polymorphism in the Hickory groups so we removed it from further analyses (Supporting Information S1). Four loci (sfl006, sat7, sfl008, and A121) deviated from HWE in eastern cottontails following Bonferroni correction and 1 (sat 8) was monomorphic in the Sandusky group (Supporting Information S1). We removed these 5 loci from further analyses of eastern cottontail samples for a total of 10 microsatellite markers and retained all loci but sfl008 for the swamp rabbit analysis for a total of 15 microsatellite markers (Supporting Information S2). No linkage disequilibrium was indicated between any loci.
We identified 13 unique swamp rabbits from genotypes recovered from fecal pellets (Fig. 1). These genotypes were included in population-level analyses. The P IDSIB for a full suite of loci was 2.9 × 10 -5 . For those pellets that produced genotypes identical to another pellet sample with ≥ 7 and < 15 loci genotyped, the maximum P IDSIB was 0.016. This gave a minimum confidence of 98% that the pellets originated from the same individual (Mills et al. 2000). The false allele rate for pellets ranged from 0-5% and averaged 2%. Allelic dropout errors per locus ranged from 16-42% and averaged 26%, comparable to other fecal DNA studies (Frantz et al. 2003).
For swamp rabbits, STRUCTURE produced the highest lnP(D) at K = 5; however, the trend of increasing lnP(D) plateaued at K = 4 with little increase for K = 5 (Fig. 2). For the swamp rabbit data set, the DIC produced by the algorithm of TESS appeared to level after K MAX = 4 (Fig. 2). Admixture plots also stabilized at K = 4. Assignment probabilities (Q) of swamp rabbits to genetic groups calculated with STRUCTURE showed a pattern related to location (Fig. 3A) including southwestern, south-central, north-central, and northeastern groups. All locations appeared to contain some individuals with mixed membership, usually from the neighboring cluster, but a few exceptions did exist. Results from TESS showed a similar pattern, but most samples had Q-values that indicated a stronger correspondence to geographic location (Fig. 3B). However, this was not the case with the north-central cluster in which Q-values obtained with TESS were generally lower (< 0.5) than individuals in the southeast, south-central, and northeast clusters (> 0.8). The northeast genetic cluster delineated by TESS showed a similar pattern to that obtained with STRUCTURE, and indicated that swamp rabbits from Bay Creek were more genetically similar to samples in the main channel of the Cache River (i.e., the north-central cluster) than to those directly west in the northeast cluster (Figs. 3 and 4).
For eastern cottontails, the highest lnP(D) was found at K = 1 (Fig. 2). For eastern cottontails, TESS does not run simulations for K MAX = 1, so we examined admixture plots for K MAX = 2 after which DIC increased. Three eastern cottontails had assignment probabilities to a 2nd cluster > 0.5. This pattern did not correspond to admixture results from STRUCTURE when K = 2. Because evidence for 2 genetic groups was not supported with STRUCTURE simulations and only weakly supported with TESS simulations, we regarded K = 1 as the best solution for the number of genetic groups for eastern cottontails.
Given the results of Bayesian clustering, locations were grouped into 4 subpopulations (Table 1, Fig. 4; Supporting Information S2). For swamp rabbits, global F ST was 0.067 (P < 0.001) and global F IS was 0.039 (P = 0.009). Pairwise F ST values were all significantly different from 0 (P < 0.001) and ranged from 0.045-0.096 (Table 2). The largest pairwise F ST value was observed between the southwest and northcentral groups ( Table 2). The smallest F ST values occurred between neighboring groups: south-central and southwest, and south-central and north-central. For eastern cottontails, global F ST was 0.028 (P < 0.001) and global F IS was 0.035 (P = 0.023). Pairwise F ST ranged from 0.018-0.033 and all values were significantly different from 0 (P < 0.01).
The spatial autocorrelation coefficient decreased more sharply over a shorter distance for eastern cottontails than for swamp rabbits (Fig. 5). Autocorrelation decreased to a level not significantly different from 0 at 4 km for eastern cottontails. Autocorrelation coefficients for eastern cottontails were significantly less than for swamp rabbits at the 1-, 3-, and 4-km lag distances (Fig. 5). For swamp rabbits, autocorrelation was significantly greater than 0 (i.e., 95% CI did not overlap 0) until 7 km (r = −0.16).

Discussion
The swamp rabbit and eastern cottontail conform to the expectation that habitat specialists have greater genetic structure than co-occurring habitat generalists. We found 4 genetic groups of swamp rabbits within approximately 50 km along the Cache River and nearby watersheds in southern Illinois compared with evidence for only 1 genetic group of eastern cottontails in the same area. In general, swamp rabbit groups showed a moderate level (0.05 < F ST < 0.10) of genetic differentiation. Genetic groups were distributed among sections of the Cache River: the southwest, south-central, north-central, and northeast. Eastern cottontails that were delineated into similar groups for comparison had a lower level of genetic differentiation (F ST < 0.05).
Lower N e and genetic diversity of swamp rabbits may have been partly due to demographic effects at the range margin. However, factors that could have contributed to the fine-scale genetic structure we observed in our spatial autocorrelation analysis included patch isolation, restricted dispersal corridors, and habitat shape. Habitat patches for swamp rabbits are farther apart than for eastern cottontails in Illinois (Illinois Natural History Survey 2014) and more distant than the maximum movements (<1 km) of most individuals (Shields 1960;Chapman and Trethewey 1972;Scribner and Warren 1990). Thus, lack of connectivity due to patch isolation (Fig. 4) may have contributed to the observed genetic structure of swamp rabbits. Habitat use and, likely, dispersal of swamp rabbits are restricted to riparian corridors (Scharine et al. 2011), whereas eastern cottontails can readily occupy and move through early successional upland habitats, row crop agriculture, and even urban areas (Chapman et al. 1980;Mankin and Warner 1999;Bond et al. 2002). Consequently, eastern cottontails likely have more potential dispersal routes as well as shorter distance between habitat patches. Finally, greater genetic structure may be due to the linear nature of swamp rabbit habitat. Under a stepping stone model, greater genetic differentiation occurs in a 1-dimensional scenario (i.e., a linearly distributed population) than in a 2-dimensional scenario (i.e., a population on a plane- Kimura and Weiss 1964). Eastern cottontails occupy more of our study area (Scharine et al. 2011) and, due to their generalist nature, they can contact more home ranges in all directions leading to less inbreeding, greater N e , less genetic drift, and less genetic structure than swamp rabbits. Thus, the  (Pritchard et al. 2000) is shown in the top graph and the deviance information criterion (DIC) from TESS (Chen et al. 2007) is shown on the bottom. Analyses were performed with swamp rabbit (Sylvilagus aquaticus; diamond markers) and eastern cottontail (S. floridanus; square markers) samples obtained from southern Illinois, 2004Illinois, -2011 difference in habitat dimensionality among the 2 species likely contributed to the difference in genetic structure and, in turn, to the differences in genetic diversity and N e .
Several studies have examined population genetic structure of leporids in a geographical context (Litvaitis et al. 1997;Crouse et al. 2009;Estes-Zumpf et al. 2010) but, until recently, few have provided evidence that anthropogenic land use may be responsible for genetic structure. Although habitat fragmentation through land use changes was hypothesized to cause genetic structure in the pygmy rabbit, strong evidence of  genetic structure among populations was lacking (Estes-Zumpf et al. 2010). Populations separated by < 30 km were interbreeding despite intervening roads and perennial streams. European brown hares (Lepus europeus) also showed no evidence of genetic structure within a large heterogeneous study area separated by highways (Fickel et al. 2005). However, due to the recent availability of microsatellite markers for North American rabbits for use in landscape genetic analysis (Estes-Zumpf et al. 2008;Berkman et al. 2009), evidence is mounting that indicates anthropogenic land use changes can alter the genetic structure of imperiled rabbits (Fenderson et al. 2014). Our study demonstrates genetic structure in a leporid (i.e., swamp rabbits) may be attributed to such anthropogenic changes including alterations in hydrology that impact the spatial structure of wetland habitat. In particular, anthropogenic activities have altered or destroyed the bottomland hardwood forests swamp rabbits prefer, threatening their persistence at the edge of their range (Terrel 1972;Korte and Fredrickson 1977;Whitaker and Abrell 1986;Kjolhaug et al. 1987;Scharine et al. 2009). Much of the remaining habitat in Illinois is highly fragmented or too small to maintain a viable swamp rabbit population (Woolf and Barbour 2002). Furthermore, loss of early successional habitat in the eastern United States due to forest maturation may continue to reduce available habitat for swamp rabbits and eastern cottontails (Litvaitis 2001;Scharine et al. 2009). We considered factors that may have produced the genetic patterns we observed. Expanses of unsuitable land cover may have restricted gene flow and isolated some groups such as the northeast group (Fig. 4), but such expanses did not occur between the southern groups. We found evidence of swamp rabbit gene flow between the Bay Creek and Cache River watersheds, specifically among the separate groups of the north-central subpopulation (Fig. 4). Bay Creek samples were expected to be most similar to the northeastern subpopulation due to proximity, but this was not the case. The relationship of the Bay Creek swamp rabbits to the north-central subpopulation of the Cache River may reflect historical habitat configuration and the route of water flow through the Cache River Valley prior to 1950 (Gough 2005 ; Fig. 4). In the early 20th century, the Ohio River and Bay Creek floodwaters flowed west into the Cache River. After the installation of the Bay Creek relief channel and a levee near Reevesville, Illinois in 1952, this connection was cut and the subsequent flood protection allowed the conversion of wetlands to agricultural use (Demissie et al. 1990). Unfortunately, we were unable to recover genotypes from an occupied area intervening the 2 north-central groups so limited sampling in this area precludes any strong conclusions.
The relatively restricted dispersal of swamp rabbits apparent in our spatial autocorrelation analysis means that swamp rabbits likely mate more often with neighbors, which can produce isolation-by-distance (Wright 1943). When sampling across a landscape is clustered rather than random and isolation-bydistance exists, the algorithm of STRUCTURE may identify discrete clusters when no barriers to gene flow exist (Serre  below diagonal) and 10 microsatellite loci in 81 eastern cottontails (S. floridanus; above diagonal) from southern Illinois, during 2004-2011. All values were significantly different from 0 (P < 0.01) as determined by 999 random permutations. Only 4 eastern cottontails were assigned to the southwest population making an inadequate sample size, thus they were assigned to the south-central subpopulation.  Schwartz and McKelvey 2009). However, models that take into account spatial location of samples can provide robust estimates of genetic structure despite the sampling scheme (Frantz et al. 2009). Additionally, STRUCTURE and TESS have been shown to perform well under simulations of genetic clines where true population subdivision exists (Chen et al. 2007;Francois and Durand 2010). Concurrence between the results of TESS and STRUCTURE support our interpretations. Significant genetic divergence of 4 groups within 1 watershed < 1,600 km 2 in area represents a very low degree of population connectivity for swamp rabbits and has implications for population persistence in southern Illinois if habitat connectivity decreases further. Loss of habitat leading to the extinction of 1 genetic cluster would decrease genetic diversity of swamp rabbits inhabiting Illinois. In turn, decreased connectivity between the remaining groups would likely increase genetic differentiation and impact metapopulation viability as a whole (Hilty et al. 2006). Furthermore, very little habitat for swamp rabbits occurs between the Cache River watershed and occupied areas further north (Rubert 2007). Genetic divergence of swamp rabbits at sites north of the study area and in other parts of their range margin is likely (Lesica and Allendorf 1995).
Our study reinforces the relationship between relative habitat specialization and population connectivity (Jonsen and Fahrig 1997;Swihart et al. 2003) and extends the evidence to co-occurring leporids. We have highlighted the need to account for the degree of habitat specialization and placement within a species range rather than using superficial similarities among organisms to make inferences about population connectivity and gene flow. Finally, due to the potential for linearly distributed habitat to exacerbate genetic drift (Kimura and Weiss 1964), we express the need to account for habitat shape when considering restoration goals of improving gene flow and genetic diversity for swamp rabbits.