Introduction
Two non-sister species of woodrats (Neotoma floridana and N. micropus; Edwards and Bradley 2002; Longhofer and Bradley 2006) occur parapatrically from the Gulf of Mexico to southeastern Colorado (Figure 1; Hall 1981). Although both species can be found in a variety of habitats, N. floridana typically occupies more mesic riparian habitats (Wiley 1980), whereas N. micropus generally exploits more shrub-like, xeric habitats (Braun and Mares 1989). The distributions of these two species are separated by a few kilometers at several localities, and by less than one kilometer at others (Spencer 1968; Birney 1973, 1976; Stangl et al. 1992; Schmidly 2004; Mauldin et al. 2014, 2021). Based on the results of morphologic, allozymic, karyotypic, and genetic data, previous studies determined hybridization occurred along the North Canadian River in Major County, Oklahoma (Spencer 1968; Birney 1973, 1976; Mauldin et al. 2014).
Recently, Mauldin et al. (2021) examined genotypic variation in individuals collected from Major and Woodward counties and reported that hybridization was intermittent with potentially transient contact zones in this region, as evidence of genetic introgression was present at 11 of 12 sampled localities. Despite this apparent widespread evidence of genetic introgression, only two localities contained mitochondrial DNA (mtDNA) haplotypes of both species and individuals with highly admixed nuclear genomes (Mauldin et al. 2021). Additionally, two temporal sampling events (separated by 22 years) from the hybrid zone indicated ongoing and potentially ephemeral hybridization is occurring between the two species in western Oklahoma (Mauldin et al. 2021). Similarity of hybrid zone characteristics (i. e., location of the zone, frequency of hybrids detected, directionality of hybridization, level of population substructure detected, etc.) in both datasets indicated short term stability of the hybrid zone; however, expanded geographic sampling detected varying levels of nuclear admixture at 10 of 11 peripheral localities. Presence of individuals with N. floridana mtDNA haplotypes and N. micropus nuclear genomes at two localities west of the known area of sympatry suggested the location of the parapatric boundary between these taxa may: 1) be larger than the hybrid zone examined by Mauldin et al. (2021) or 2) there may be multiple sites of active hybridization (Mauldin et al. 2021).
Evidence of intermittent hybridization in Major County, Oklahoma (Mauldin et al. 2014, 2021), has lent support to the possibility that additional areas of hybridization may exist throughout the area of parapatry (Spencer 1968; Birney 1973). A second potential area of contact, along the south bank of the Red River (Locality 20, Figure 1) was sampled at intervals over several years (Stangl et al. 1992). Although no morphological evidence of hybridization was reported, Stangl et al. (1992) collected N. floridana and N. micropus within 100 m of each other, thereby establishing the possibility that the two species were in contact. Superficially, this region is similar to that of the known hybrid zone in Major County (Spencer 1968; Birney 1973; Mauldin et al. 2014, 2021), as the Red River bisects the parapatric border of these species, and riparian habitat typically exploited by N. floridana interdigitates with sage brush and sand dunes, more commonly inhabited by N. micropus. In addition to current areas of hybridization, detection of admixed individuals at localities peripheral to the current parapatric boundary could provide insight into the stability of the distributions of these species, and the effect dynamic distributions may have on hybridization in this system.
Given the potential ephemeral nature of the previously studied hybrid zone, along with the long parapatric border shared by these species, Mauldin et al. (2021) advocated for further taxonomic sampling along the border of parapatry. They suggested further study was need to determine if 1) additional areas of hybridization exist and 2) evidence of dynamic species distributions could be substantiated. Therefore, the goal of this study was to examine potential areas of sympatry for evidence of hybridization, and to inspect areas peripheral to the parapatric border for evidence of genetic introgression. To this end, multiple objectives were addressed: 1) collect and genotype individuals from localities along and at varying distances from the proposed parapatric border, 2) examine localities for presence or absence of evidence of genetic introgression, 3) determine the maximum recorded distance of hybrid individuals from the current estimated border of parapatry, and 4) utilize microsatellite data to examine population substructure and level of genetic introgression in areas of sympatry.
Materials and methods
Samples. State collecting permits, as well as permission of property owners or appropriate state agencies (e. g., Kansas Department of Wildlife and Parks, Oklahoma Department of Wildlife Conservation, Texas Parks and Wildlife) were received prior to any collection efforts. One hundred and forty woodrats were collected from 21 localities throughout Kansas, Oklahoma, and Texas (Figure 1) between July 2009 and May 2012. Spatial distribution of individual capture sites (middens) were identified with UTM coordinates. Most woodrats were collected with Sherman live-traps (Sherman live-trap Co. Tallahassee, Florida), others were collected with Havahart ® live-traps (Woodstream Corporation, Lititz, Pennsylvania, USA), and some were captured by hand after excavation of middens (nests) to ensure all occupants were collected. Individuals and embryos of sufficient size to ensure extraction of embryonic DNA were given a unique identification number (TK number), sexed, measured, and sacrificed. Individual woodrats were assigned putative species identifications based on morphologic characteristics (Hall 1981; Schmidly 2004), however, given previous results and the inability to distinguish hybrids based solely on morphology, a formal morphological identification based on pelage color was not considered in hybrid identification. Animal care and use guidelines conformed to those proposed by the American Society of Mammalogists (Animal Care and Use Committee 1998) and were approved by the Texas Tech University Institutional Animal Care and Use Committee (IACUC protocol 11009-03). In cases where females and their offspring were captured in the same midden, mother and offspring were cross-referenced; similarly, pregnant females were cross-referenced to embryos. Blood and tissue samples (heart, kidney, liver, lung, muscle, and spleen) were obtained and tissues were immediately frozen in liquid nitrogen, and subsequently archived at the Natural Science Research Laboratory (NSRL) at the Museum of Texas Tech University. Voucher specimens (skulls, postcranial skeletons, and skins) were prepared and deposited in the NSRL (Appendix 1). Additionally, liver samples of four woodrats (indicated by prefix TJM in Tables A1 and A2) were obtained from the lab of Ivan Castro-Arellano at Texas State University.
DNA Isolation. Total genomic DNA (nuclear and mitochondrial) was isolated from each individual using approximately 0.1 g of liver and the Qiagen DNeasy kit (Qiagen Inc.; Valencia, California, USA). In some cases, entire embryos were required to isolate sufficient DNA. DNA samples were stored at -20° C for subsequent analyses.
Genotype Analyses. All genotype analyses followed the protocol outlined in detail by Mauldin et al. (2014, 2021). Eight N. floridana and seven N. micropus collected a minimum distance of 125 km from the parapatric border, and previously utilized by Mauldin et al. (2014, 2021) were included as reference samples (Figure 1; Appendix 1). Three loci were examined, two autosomal loci (intron two of the vertebrate alcohol dehydrogenase gene (Adh1-12) and intron seven of the beta-fibrinogen gene (Fgb-I7)) and one mitochondrial DNA locus (Cytochrome-b, Cytb). Additionally, individuals from Locality 20 were genotyped for six microsatellite loci (Nma01, Nma04, Nma05, Nma06, Nma10, and Nma11) developed by Castleberry et al. (2000) to detect genetic structure within the population.
Adh1-12 Assay. A banding pattern unique to N. floridana was produced using the restriction enzyme NsiI (ATGCA/T) with a fragment of the Adh1-12 region of either 566 bp or 390 bp that had been amplified using PCR methods modified from Amman et al. (2006) and Longhofer and Bradley (2006) using one of the following primer pairs: ExonII-F and 2340-II (566 bp product) or 350F and 2340-II (390 bp product; Amman et al. 2006). Restriction digests were conducted following manufacturer’s methods and are outlined in Mauldin et al. (2014).
Fgb-I7 Assay. Mauldin et al. (2014) reported that although no restriction enzyme was diagnostic, three diagnostic nucleotide substitutions were identified (positions 428, 497, and 493). Therefore, sequence data was collected on a 609-610 bp fragment amplified using PCR primers Fgb-I7L-Rattus and Fgb-I7U-Rattus from Wickliffe et al. (2003) and following PCR methods modified from Prychitko and Moore (2000) as outlined in Carroll and Bradley (2005). Sequence data has been deposited in GenBank (Appendix 1). Though previous studies have utilized this as a diagnostic marker (Mauldin et al. 2014, 2021), it is possible unsorted polymorphisms not detected in reference samples may exist.
Cytb Assay. The entire Cytb gene was amplified using two PCR primers (LGL765 forward-Bickham et al. 1995 and LGL 766 reverse-Bickham et al. 2004) and conditions outlined in Edwards and Bradley (2002). The restriction enzyme [BsaI (GGTCTC(N)1/)] produced a cut that was unique to N. floridana following methods outlined by the manufacturer and reported by Mauldin et al. (2014, 2021).
Microsatellite Assay. The six microsatellite loci developed by Castleberry et al. (2000) and utilized by Mauldin et al. (2014, 2021) were amplified and analyzed for all individuals collected at Locality 20 (Figure 1) as described by Haynie et al. (2007). Alleles were scored using GeneMapper software (version 4.0; Applied Biosystems Inc.).
Data Analysis. Based on the results of molecular assays outlined above, each individual was scored as either N. micropus or N. floridana for the mitochondrial genome, and as homozygous N. micropus, heterozygous, or homozygous N. floridana for the Adh1-I2 and Fgb-I7 markers. GenAlEx (version 6.5; Peakall and Smouse 2012) was utilized to identify presence of duplicate genotypes, test microsatellite loci for deviation from Hardy-Weinberg equilibrium (HWE) expectations, and format data for use in Structure (v2.3.4; Pritchard et al. 2000) as codominant nuclear markers with only adults and subadults being included in the analyses. Mitochondrial data were not analyzed in Structure or NewHybrids (v.1.1Beta3; Anderson and Thompson 2002) but were included in result plots to aid in identification of hybrid individuals and examine any potential bias present in directionality of introgression.
Based on preliminary results (nuclear introgression and geographic proximity of both mtDNA haplotypes), complete nuclear genotypes (Adh1-I2, Fgb-I7, and six microsatellite loci) for all individuals collected at Locality 20 were analyzed in Structure to examine population structure and quantify potential admixture between the two species. Structure runs utilized the admixture model with independent allele frequency option, a burnin of 500,000, run length of 1,000,000 iterations, and examined values of K (clusters) from 1-5. Two separate parameter sets were run, one assigned reference individuals to a priori populations using the popflag designation (parameter set A), whereas the other did not (parameter set B). Neither dataset used prior population assignment information for study samples. Structure result files were uploaded to Structure Harvester (Earl and vonHoldt 2012) to determine the value of K which best fit the data using the Evanno method (Evanno et al. 2005).
Results of the Structure run with the smallest variance value from parameter set A (K = 2) were used to generate a plot for examination of admixture between the two species. Furthermore, individuals from Locality 20 were analyzed in NewHybrids to determine the posterior probability values (PPVs) of individuals belonging to one of six classifications (pure parental N. floridana, pure parental N. micropus, F1, F2, backcross to N. floridana, backcross to N. micropus) based on admixture of nuclear genomes with no prior allele frequency data, Uniform priors, a burnin of 100,000, and 1,000,000 sweeps after burnin. Structure and NewHybrids output files were visualized using Excel 2010 (Microsoft Corporation, Redmond, Washington, USA). Assignment to hybrid classifications followed the protocol outlined by Mauldin et al. (2014, 2021).
Electronic species distributions of N. floridana and N. micropus (Patterson et al. 2007) generated by digitizing previously published range maps (i. e., Hall 1981) were used to approximate the location of the parapatric border. Distance of each sampling locality to the closest point along the approximated parapatric boundary was then measured with the use of ArcGIS Software (ESRI, Redlands, California, USA), based on UTM coordinates of localities. Distances were measured to each distributional boundary (N. floridana and N. micropus) along the same vector, and the two distances were averaged for the final estimate. Additionally, samples from Locality 20 were collected from two nonadjacent parcels of private property; however, given the proximity of localities (all samples collected within ~2.5 km), samples from both properties were consolidated into a single locality for simplicity. However, these localities are examined both jointly (Locality 20) and independently (Localities 20a and 20b) to better examine patterns of interspecific genetic introgression at multiple scales.
Randomization tests of goodness-of-fit utilized 20,000 iterations and were conducted with Excel 2010 (Microsoft Corporation, Redmond, Washington, USA) following methods described by McDonald (2009) to determine if the following proportions deviated significantly from an equal contribution: 1) proportion mtDNA haplotypes of each species at Localities 20, 20a, and 20b, 2) proportion of Adh1-I2 and Fgb-I7 alleles detected at localities within each presumed species distribution, and 3) proportion of Adh1-I2 and Fgb-I7 alleles detected east and west of the proposed center of the hybrid zone in Major County, Oklahoma (data from Mauldin et al. 2014). The proportion test within the statistical package R (Team 2008) was utilized to compare the following proportions: N. floridana mtDNA haplotypes detected at Localities 20a and 20b, hybrid individuals with introgression detected at the Adh1-I2 locus within the distributions of N. micropus and N. floridana, respectively, and hybrid individuals with introgression detected at the Fgb-I7 locus within the distributions of N. micropus and N. floridana, respectively.
Results
Results of molecular assays are available in Appendix 2. Evidence of mixed ancestry was detected in 55 of 140 (39 %) sampled individuals, at 10 of 21 (48 %) localities (Figure 1). A high percentage of individuals with mixed ancestry (>50 %) was recorded at three localities (4, 9, and 20). Genetic introgression was detected at both nuclear loci at only two localities (9 and 20), whereas only Locality 20 contained mitochondrial DNA (mtDNA) haplotypes of both species. Given that Locality 20 is the only site at which both mtDNA haplotypes were detected, and the possibility that it represents an area of current or recent contact (Stangl et al. 1992), individuals from Locality 20 were not included in examination of differential detection of admixture between loci. Of the 22 woodrats with mixed ancestry collected within the distributional limits of N. micropus, evidence of nuclear admixture was detected at the Adh1-I2 locus in one animal, and at the Fgb-I7 locus in 21 individuals (P < 0.0001). Admixture was detected only at the Adh1-I2 locus for all 21 admixed individuals identified within the distribution of N. floridana (P < 0.0001). No individuals were heterozygous at both loci. A similar bias was identified through use of randomization test of goodness-of-fit that examined data from the Major County hybrid zone, as detection of admixture in individuals collected west of the proposed center of the zone was significantly biased towards Fgb-I7 locus (P = 0.012), and detection of admixture east of the zone was biased, although not significantly, to the Adh1-I2 locus (P = 0.073). The furthest distance from the parapatric border at which nuclear admixture was detected was approximately 150 km within the species distribution of N. micropus (Locality 6; Figure 1). Additional distance data for localities and individuals is available in Table 1.
Category | Mean | Median | Minimum | Maximum |
---|---|---|---|---|
All localities | 42 | 27 | 4 | 152 |
hybrid localities | 45 | 29 | 12 | 152 |
‘pure’ localities | 40 | 27 | 4 | 152 |
all individuals | 34 | 26 | 4 | 152 |
hybrid individuals | 27 | 26 | 12 | 152 |
‘pure’ individuals | 38 | 26 | 4 | 152 |
Examination of microsatellite data with GenAlEx identified no duplicate genotypes. The following markers was determined to deviate significantly from HWE expectations within the sampled population, Nma05 in Locality 20b (P = 0.030). Results of Structure and Structure Harvester analyses determined K = 2 as the most appropriate number of clusters for both parameter sets. Results of Structure analyses detected no genetic introgression or population substructure at Locality 20 (Figure 2). Results of NewHybrids analyses of samples from Locality 20 identified only one sample as less than 90 % probability of belonging to the classification of ‘pure’ N. micropus (Figure 3; TK179266 = 87.76 %; mean N. micropus PPV = 96.56 %, median N. micropus PPV = 98.77 %). Spatial distribution of mtDNA haplotypes within Locality 20 is depicted in Figure 4. Results of randomization test of goodness-of-fit determined the proportion of mtDNA haplotypes of each species present at Locality 20 did not vary significantly from the null model of equal contribution (P = 0.118), nor did Locality 20a (P = 0.690); however, Locality 20b was significantly biased towards N. floridana mtDNA, with no N. micropus mtDNA haplotypes detected (P = 0.004). Results of analyses of the proportion of N. floridana mtDNA haplotypes at Locality 20b were significantly higher than that of Locality 20a (P = 0.024). Results of analyses of the proportion of hybrids with introgression detected at the Adh1-I2 (P < 0.0001) and Fgb-I7 (P < 0.0001) loci varied significantly depending upon the distribution from which they were collected.
Discussion
A high proportion of sampled woodrats (39 %) were determined to be of mixed ancestry, including individuals from 10 of 21 (48 %) sampled localities throughout Texas, Oklahoma, and Kansas (Appendix 2). Given the small number of molecular markers examined, these values likely underestimate the true number of genetically admixed individuals and localities at which they are found. These results suggest some degree of hybridization has occurred, or currently occurs, at multiple localities along the parapatric border. Additionally, at three localities (4, 9, and 20) greater than 50 % of examined individuals exhibited some level of genetic introgression; although no genetic admixture was detected at three localities (13, 14, and 19) of similar or lesser distances to the border. Furthermore, the geographic distance between some sampled hybrids and the putative location of the parapatric boundary is substantial (e. g., Locality 6: ~150 km).
Finally, the locus at which admixture was detected was dependent upon the species distribution from which the samples were collected. For individuals collected within the distribution of N. micropus, genetic introgression was detected most frequently at the Fgb-I7 locus; however, for individuals collected within the distribution of N. floridana, introgression was detected only at the Adh1-I2 marker. The statistically significant difference in the locus at which exotic alleles were detected within each species distribution suggests that selection favors inclusion of foreign DNA sequences at different loci based upon the genomic background of the organism (i. e., predominantly N. floridana or N. micropus nuclear genomes). Examination of the Major County hybrid zone data generated by Mauldin et al. (2014) identified a similar bias at a smaller geographic scale.
Results of Structure analyses failed to detect nuclear admixture at Locality 20 and estimated that nuclear genomes of sampled woodrats were predominantly (>99 %) composed of N. micropus alleles, although the majority contained N. floridana mtDNA haplotypes. Examination of results of NewHybrids analyses in combination with Adh1-I2, Fgb-I7, and Cytb data determined all individuals from Locality 20 were either backcrosses to N. micropus (12) or putatively pure N. micropus (3), with no N. floridana parental types detected. The statistically significant change in proportion of mtDNA haplotypes present between Localities 20a and 20b, combined with paucity of N. floridana parental types, and nuclear genomes of all individuals composed primarily of N. micropus alleles suggests the location of the hybrid zone has shifted from the approximate location of Locality 20a to some location east of the sampled area. The easternmost sample (TK 179251) was collected ~100 m east of the I-44 Bridge reported to be an area of contact (Stangl et al. 1992) and appeared to be N. micropus (pelage and nuclear genome). Therefore, it is possible that the shift in area of sympatry began prior to the study by Stangl et al. (1992), at which time it had reached the I-44 Bridge, and has subsequently continued east along the Red River, with the mtDNA haplotype of N. floridana occurring throughout its now displaced range.
Similar cytonuclear discordance, although smaller and directionally reversed, was reported between the positions of the mitochondrial and nuclear boundaries between these species in Major County, Oklahoma, suggesting that the areas of hybridization are somewhat transient as distributional borders of these species shift over generations (Mauldin et al. 2021). Given the large variation in degree of genetic introgression detected over relatively small distances and the small proportion of individuals detected with highly admixed nuclear genomes (F1 and F2-like individuals) reported along the North Canadian River (Mauldin et al. 2021), distance to the putative area of sympatry cannot be estimated with any certainty. However, it is worth noting that an individual identified as a putatively pure N. floridana was collected from Locality 19 (~10 km east of Locality 20) in Oklahoma.
Various methodologies, including morphologic, karyotypic, allozymic, and genotypic data have been used to examine hybridization at various geographic scales within this system (Spencer 1968; Birney 1973, 1976; Mauldin et al. 2014, 2021). Examination of these studies and the research presented herein has determined the following characteristics are demonstrative of hybridization between these species: a high percentage of genetically admixed individuals, no evidence of reduced fertility in hybrid individuals, a paucity of F1 and F2-like genotypes, significant linkage disequilibrium, limited population structure, differential genetic introgression of nuclear loci, and varying levels of hybrid zone ephemerality. Examination of these characteristics in the framework of mechanistic models of hybrid zone maintenance and criteria set forth by Endler (1977) and Moore (1977), as summarized by Van Den Bussche et al. (1993), indicate that either the hybrid equilibrium model (wherein hybrids and parental types are equally fit) or the hybrid superiority model (in which hybrids have higher fitness within an ecotone or certain set of environmental conditions) might be responsible for maintaining hybridization between these species. Additional data concerning hybrid fitness, selection pressures, and possible correlations to environmental conditions are needed to distinguish between these models relative to the dynamics responsible for maintaining hybrid zones between these species.
Historical distribution changes of N. floridana have been documented by Quaternary fossil records (Richards 2013), ‘recent fossil’ remains dating to the late Holocene (~1,450 years before present - Eshelman 1971; Richards 2013), and temporal sampling in the 19th century (Cope 1872; Blatchley 1897). Additionally, a study examining distributions of N. micropus and N. albigula in southern New Mexico, determined interspecific competition led to displacement of N. micropus over a portion of the study area (Wright 1973). Given the evidence presented herein, the ephemeral nature of distributional boundaries, and the documented occurrence of interspecific displacement within the genus Neotoma, it is possible that the evidence of introgression detected at peripheral localities is a result of some combination of distributional shifts and dispersal of alleles over generations. Subsequently, the differential detection of alleles at distinct nuclear loci might be the result of disparity in persistence of certain alleles within populations, the rate at which those same alleles disperse over generations, or some combination thereof. Additionally, the possibility of unsorted polymorphisms may exist, potentially impacting the nuclear introgression calculations; however as all but one of the molecularly identified hybrids exhibited cytonuclear discordance, this would not change the overall results or the classification for most animals examined.
In conclusion, nuclear introgression was detected at multiple localities throughout a large portion of the parapatric border including sites near Burkburnett, Texas, Seiling, Oklahoma, as well as Great Bend and Syracuse, Kansas, among others (see black dots in Figure 1). Additionally, this introgression appears to be variable with regard to prevalence of admixture detected at separate nuclear markers dependent upon the genomic background of the organism, as the N. micropus genome appears to tolerate N. floridana alleles at the Fgb-I7 locus better than at the Adh1-I2 locus, and N. floridana genome is more commonly infiltrated with N. micropus alleles at the Adh1-I2 locus than the Fgb-I7 locus. The presence of cytonuclear discordance at Locality 20, and similar evidence reported in Major and Woodward counties (Mauldin et al. 2021) provide evidence of nuclear genome displacement, likely caused by distributional shifts. Although introgression appears common throughout the parapatric border, the differential introgression of alleles and paucity of individuals determined to have highly admixed nuclear genomes, suggest hybridization does not pose a major threat to the gene pools of either species.