Discordance between taxonomy and population genomic data: An avian example relevant to the United States Endangered Species Act

Abstract

Introduction

The Endangered Species Act (ESA) is designed to help protect species of conservation concern, and it turned 50 years old in 2023. Many authors have pointed out the practical problems with language in the ESA, particularly as it pertains to hybrids (e.g. (1–4)). These problems are exacerbated by the rapidly increasing resolution of population genomic data that allows biologists to more clearly identify and quantify signatures of effective hybridization (i.e. introgression). In sexual organisms, genomic data indicate that many related species are obvious evolutionary amalgamations produced largely by horizontal gene flow (5–7).

The U.S. Fish and Wildlife Service (USFWS) recently listed the lesser prairie-chicken (LEPC) under the ESA, identifying two distinct population segments (DPSs). The more isolated Southern DPS is now federally endangered and the Northern DPS is now threatened (Figure 1a; (8, 9)). The LEPC is one of three lekking prairie grouse species in the genus Tympanuchus, which also includes the greater prairie-chicken (GRPC) (10–12). The LEPC and GRPC diverged from one another 0.6–0.9 Mya (13, 14), and they can hybridize in regions of geographic sympatry (Figure 1a). Hybrids are apparently fertile as indicated by the direct evidence of controlled crosses that have produced F1 and F2 interspecific hybrids (15) as well as by indirect genetic evidence of introgression in wild birds (16, 17).

Fig. 1.

Geographic sampling regime and Principal Component Analysis (PCA) of whole-genome sequences. a) Sampling locations of 411 LEPC (T. pallidicinctus) collected from the Northern and Southern DPSs in the LEPC current range. Additionally, 22 GRPC (T. cupidio) were collected from sympatric (within LEPC range) and allopatric (outside of LEPC range) locations. Sampling points largely correspond to sampling lek, with 1–20 birds collected from each, but some individuals were obtained outside of leks. The black borders in the Northern DPS and Southern DPS correspond to the four different ecoregions. b) PCA of allele frequencies derived from ∼ 18 M polymorphic (maf > 0.01) sites across the nuclear genome. Principal component 1 is labeled on the Y-axis, and principal component 2 is labeled on the X-axis to better orient the points geographically. Less than 10% of the overall variation is explained by taxonomy or geography. Above the PCA, individual admixture proportions (K = 3) for mixed ancestry samples (> 0.20) identify three birds (at ends of bar plot; see Figure 2 for additional information) where genomic identities did not correspond to morphological species assignments (sample numbers included in white font). For symbols in the map and PCA, colors in the diagrammed silhouettes correspond to the phenotypic species identification of LEPC (brown) and GRPC (goldenrod). The dashed diagonal (a) and horizontal line (b) depicts the separation between the Northern and Southern DPSs.

Open in new tabDownload slide

Previous genetic studies of LEPC using a few genetic markers revealed moderate genetic differentiation and similar diversity levels in samples from different ecoregions (11, 16, 18, 19), but the recent USFWS (8) Final Rule states that “There are concerns about the implications of genetic introgression (dilution) of LEPC genes, particularly given that potential effects are poorly understood.” Here, we performed a high-resolution assessment of LEPC population structure and potential interspecific introgression that are relevant to new conservation policies. Our whole-genome sequence data from 433 individual LEPCs and GRPCs demonstrate that: (i) the LEPC has modest levels of genomic diversity compared to other avian species; (ii) genomic diversity has not substantively changed over time in either species; (iii) there is distinct but moderate differentiation between LEPC sampled from the Northern and Southern DPSs; (iv) interspecific hybridization between LEPC and GRPC has resulted in bidirectional introgression that blurs assignment of a few individuals to a species; and (v) interspecific divergence between LEPC and GRPC is very low. The boundary between these two species is still porous (11), and the current scientific nomenclature (i.e. discrete Latin binomials for each species), which is critical for effective conservation policy, poorly reflects their ongoing evolutionary dynamics.

Results

Population structure and introgression

LEPC samples were collected over more than a decade from the Southern DPS (n = 150) and the Northern DPS (n = 261); GRPC samples were collected from areas of sympatry with LEPC (n = 13) and areas of allopatry (n = 9; Figure 1a). A principal component analysis (PCA) and admixture analysis were performed on ∼ 18 M variable sites among 433 individual samples. A distance-based approach (i.e. model-free) using allele frequencies across the entire genome explained only a small proportion of the overall variability in the dataset (principal component 1 = 3.60%, Y-axis). The PCA clearly separated the Southern DPS and revealed a transect of genetic connectivity between LEPC from the Northern DPS to the allopatric GRPC (Figure 1b). Samples along principal component 2 (2.77% explained variation, X-axis) poorly discriminated LEPC and GRPC sampled north of the DPS line (Figure 1b). A virtually identical PCA was generated when using the chicken (Gallus gallus) reference genome (GCA_016700215.1), demonstrating that there was minimal ascertainment bias due to reference genome selection (FigureS1). Admixture analysis using a model-based (Hardy–Weinberg) approach identified two genomic groups (K = 2; Figure S2) that clustered individuals by collection location above or below the official DPS boundary line but did not discriminate between phenotypic species assignment in the Northern part of the range (Figure 2a). However, the next most likely number of genomic groups (K = 3) categorized individuals based upon DPS and species, revealing bidirectional interspecific introgression (63 individuals with co-ancestry proportions > 0.001; Figure 2a). The K = 3 model revealed taxonomic discrepancies between field (phenotypic) species assignments and genomic affiliations in two GRPC individuals and one LEPC. One nominal LEPC from the northern DPS (F430) was apparently misidentified or misrecorded as a GRPC. One sample (F1270) identified in the field as a putative hybrid was of GRPC (0.73) and Northern DPS LEPC ancestry (0.27).

Fig. 2.

Admixture plots of whole-genome sequence data. Individual admixture proportions for each bird, illustrating the assignment to the optimal number of clusters (K = 2) and next best model (K = 3; Figure S2), according to the method of Evanno et al. (20). Colors of individual bars correspond to genetic clusters, not the phenotypic species assignment. Admixture proportions illustrated in bar plots for K = 2 (a) and K = 3 (b) are sorted from north to south latitude (left to right). For each bar plot, individual ancestry coefficients were also continuously plotted over geographic space (maps). The vertical white solid lines in the bar plots correspond to the separation between field assigned species, with 411 LEPC (T. pallidicinctus) collected from the Northern and Southern DPSs and 22 GRPC (T. cupidio) that based upon collection local were sympatric (within LEPC range) and allopatric (outside of LEPC range). For K = 2, the admixture analysis groups the Southern DPS separately from the Northern DPS + GRPC. The second most likely number of clusters (K = 3) is mostly consistent with DPS and taxonomic integrity while identifying biparental introgression between LEPC and GRPC.

Open in new tabDownload slide

Of the remaining 410 LEPC samples, 63 (15.4%) individuals had some degree of GRPC hybrid ancestry. These hybrid ancestry individuals were collected from all three ecoregions in the Northern DPS, but primarily from the Short-Grass Prairie/Conservation Reserve Program (CRP) Mosaic Ecoregion. The 63 birds with hybrid ancestry comprised 24.1% of the putative LEPCs that we sampled in the Northern DPS. The mean proportion of Northern DPS LEPC genomes attributable to GRPC (i.e. the degree of introgression) was 4.5% among these 63 individuals. Out of 13 sympatric GRPCs, six had LEPC hybrid ancestry but two were likely misidentifications in the field. The remaining four admixed GRPCs have a mean of 42% of their genome from Northern DPS LEPC. The ABBA–BABA tests and associated D-statistics found significant but modest levels of gene flow among all examined groups, including between species (Table S1). Global weighted estimates of pairwise F_ST between LEPC collected from the Northern and Southern DPSs were low to moderate (F_ST = 0.044), demonstrating that 4% of the total variation across the genome was partitioned between the two DPSs. The greatest differentiation occurred between LEPC from the Southern DPS vs allopatric (F_ST = 0.145) and sympatric (F_ST = 0.095) GRPC. As expected with interspecific gene flow, levels of differentiation were reduced when comparing LEPC from the Northern DPS to allopatric (F_ST = 0.109) and sympatric (F_ST = 0.06) GRPC (Table 1).

Genome-wide patterns of diversity

Genome-wide mean heterozygosity (H) was 0.0031 ± 0.0005 (Figure 3a) across all 411 nominal LEPC, with reduced levels of diversity in the Southern DPS (H = 0.00295 ± 3.27e – 4) compared to the Northern DPS (H = 0.00340 ± 5.30e – 4). There were no obvious trends in temporal H among samples collected from 2008 to 2023, either between DPSs or species (Figure 3b). Overall LEPC inbreeding was low (fROH_total = 0.067 ± 0.062) and similar to Greater Sage-Grouse (Centrocercus urophasianus; fROH_total = 0.072 ± 0.06), Gunnison Sage-Grouse (Centrocercus minimus; fROH_total = 0.03 ± 0.025), and other galliform species (21). Inbreeding was statistically greater in the Southern DPS than in the Northern DPS (P-values < 2e – 16). One LEPC individual (out of 411 in total) sampled in the Northern DPS was inbred (fROH_total = 0.61). No significant differences in inbreeding were detected when comparing the two species for any fROH category (P-values = 0.056–0.83; Figure 3c).

Fig. 3.

Genomic diversity and autozygosity. 411 nominal LEPC (T. pallidicinctus) were collected from the Northern and Southern DPSs and 22 GRPC (T. cupidio) that based upon collection local were sympatric (within LEPC range) and allopatric (outside of LEPC range). a) Box plot of 433 individual heterozygosity (H) estimates derived from whole-genome resequencing, grouped by sampling region (allopatric, sympatric, Northern DPS, and Southern DPS), and colored by species. b) Box plot of individual heterozygosity estimates for each sampled region, factored by the year the sample was collected. Note that several birds were missing the sampling year in the attribute data, so were labeled as NA (not applicable). There is very little temporal variation in heterozygosity regardless of geographic or taxonomic affiliation. c) The proportion of each individual genome with runs of homozygosity (fROH), binned by length (100 kb–1 Mb; > 1 Mb; total) and grouped by sampling region. For all Panels (a–c), colors in the diagramed silhouettes correspond to LEPC (brown) and GRPC (goldenrod) samples.

Open in new tabDownload slide

A distance-based nuclear tree largely resolved each species (i.e. it clustered GRPC together within a clade) and DPS, but it revealed some GRPC interspersed among LEPC collected from the Northern DPS (Figure 4a). Mixing of putative species in the nuclear tree can be partially explained based upon the 16 samples that were misidentified in the field or showed hybrid ancestry (Figure 1b). In contrast, the maximum likelihood mtDNA tree provided no resolution as it detected no clear structure between LEPC from either DPS or between species and neither did the haplotype network (Figures 4b and S12). There was no definitive evidence of reciprocal monophyly, one common yardstick of “good” species (22), in our analyses of mtDNA genomes (Figure 4). Haplotype diversity statistics associated with the mtDNA genomes reveal reduced haplotype and nucleotide diversity for LEPC sampled from the Southern DPS, consistent with the nuclear genomic results (Table 1). Furthermore, despite the larger sample size, LEPC contained less than half the level of haplotype diversity observed in GRPC (Table 1).

Fig. 4.

Nuclear and mitochondrial phylogeographic trees. Midrooted distance tree derived from filtered nuclear genotype likelihoods (a) and a maximum likelihood tree from aligned reconstructed mitochondrial assemblies (b), with colors representing the two species and shapes corresponding to the source location (Northern DPS, Southern DPS, allopatric, and sympatric). For both trees, colors in the diagrammed silhouettes correspond to the LEPC (T. pallidicinctus) and the GRPC (T. cupidio). The red arrows in a) identify the division between LEPC collected from the Northern and Southern DPSs, and the labeled sample numbers represent the potentially misidentified or introgressed birds illustrated in the bar plot above Figure 1b. The goldenrod arrows in b) are illustrated to orient the location of several GRPC samples throughout the mtDNA tree. The nuclear tree (a) shows increased species discrimination compared to the mitochondrial tree (b), but neither exhibits a clear pattern of reciprocal monophyly that is often associated with distinct species or evolutionarily significant units (22) driven in part by bidirectional introgression between species.

Open in new tabDownload slide

Species divergence and effective population sizes

Grouse have traditionally been classified using phenotypic and behavioral characteristics, but species boundaries are ill-defined due to hybridization among species (17, 23). To focus on the greatest divergence between LEPC and GRPC, the two extremes of our sampling distribution were utilized (i.e. Southern DPS of LEPC vs allopatric GRPC). A nonoverlapping 50 kb sliding window analysis of genome-wide pairwise nucleotide diversity (ø_π) showed relatively stable and highly correlated (Pearson's r = 0.968) ø_π for both LEPC (mean = 0.00361, range = 0.00003 – 0.03952) and GRPC (0.00386, 0.00036 – 0.03683), illustrating similar genomic diversity and architecture (Figure S3). Mean ø_π among both species was 0.00366 (0.00012 – 0.03609), signifying there were not any more notable differences occurring between species relative to within species. Under neutrality, ø_π = 4N_e(µ) and equivalently, N_e = ø_π / 4µ, where N_e is the effective population size and µ is the mutation rate. Using a grouse-specific mean mutation rate of 3.28 × 10⁻⁹ (24), we estimate a long-term N_e = 276,152 in LEPC and N_e = 294,207 in GRPC.

For context, we quantified divergence and differentiation using data from the Gunnison and Greater Sage-Grouse (25). Those results (Figure S4) reveal substantially greater diversity observed among the Greater Sage-Grouse population (ø_π = 0.00211) relative to the Gunnison Sage-Grouse populations (ø_π = 0.00081). Thus, Sage-Grouse have roughly half as much nucleotide diversity as prairie-chickens. Low levels of interspecific differentiation (F_ST = 0.073 across ∼18 M sites) existed between all LEPC and all GRPC samples, whereas the two Sage-Grouse species showed ∼5 × more interspecific genomic differentiation (F_ST = 0.384).

Discussion

Conservation policies ostensibly rely on the best available science. Whole-genome resequencing studies like ours have the capacity to capture relevant patterns of genomic diversity and to clarify the underlying evolutionary processes. Genomic diversity is positively associated with evolutionary fitness, with a lack of diversity often seen in populations of conservation concern (11, 26). Our study of whole-genome sequences has revealed no obvious red flags in observed levels of diversity in the LEPC. We explicitly note that this situation could change quickly in the face of rapid and severe demographic bottlenecks due to external factors such as infectious disease, additional habitat loss, and climate change. We certainly advocate the continuation of best management practices to help sustain LEPC populations for the foreseeable future.

Neither LEPC DPS is depauperate of genomic diversity, neither is particularly inbred, temporal diversity is stable in each DPS (Figure 3), and diversity is similar across ecoregions (Figure S5), almost certainly because of the large (relative to contemporary size estimates (27)) long-term N_e. There is modest but distinct genome-wide differentiation (Table 1) that reflects both the anthropogenic habitat fragmentation that has gradually isolated each DPS since the 19th century (28) and any natural interspecific isolation that might have arisen prior to anthropogenic habitat modification. Our data reveal lower levels of heterozygosity (Figure 3a) and haplotype diversity (Table 2) in the Southern DPS compared to all other examined groups. Mean levels of genomic diversity (H = 0.003) observed in our prairie-chicken samples are similar to Greater Sage-Grouse (H = 0.002) as well as many other common avian species (29). We found limited evidence for inbreeding among LEPC samples (mean fROH_total = 6%), which were similar to those estimated in Gunnison (fROH_total = 3%) and Greater Sage-Grouse (fROH_total = 7%; (25)) but substantially lower than in a critically imperiled avian species like the Stewart Island Kākāpō (fROH_total ∼80%; (30)). The lack of appreciable inbreeding and the modest level of temporally stable diversity in our dataset (Figure 3) indicates that genomic imprints of demographic reductions and geographic isolation (12, 31) have not yet manifested themselves in prairie-chickens (though see (32)).

Hybridization and subsequent backcrossing that leads to introgression is well-documented in many avian sub(species) and can be resolved using genomic approaches (32, 33). According to the standard method for identifying the optimal K-value (20), the most likely model of admixture identified two main groups (K = 2) of prairie-chickens with relatively little differentiation among birds sampled north of the DPS line (Figure 2a). This lack of differentiation is likely driven by recent or ongoing gene flow among leks, ecoregions (Figure S6), and between putative species (Figures 1b and 2a, b). A transect of genomic variability between allopatric GRPC and LEPC collected from the Northern DPS (Figure 1a) and our phylogeographic trees (Figure 4a and especially b) reveal that individual GRPC grouped with LEPC from the Northern DPS. We cannot disregard the influence of incomplete lineage sorting in the observed topography, but the lack of clear monophyletic groups is consistent with the ABBA–BABA tests and associated D-statistics (Table S1).

Our second-best model of admixture identified three main groups (K = 3) of prairie-chickens (Figure 2b). To us, the K = 3 model makes the most biological sense because it best demarcates each species and each DPS while also quantifying bidirectional, interspecific introgression. Much of the interspecific gene flow is likely due to ongoing hybridization and introgression between species in the Short-Grass Prairie/CRP Mosaic Ecoregion of the Northern DPS (Table S1; Figure S7). Lineages are expected to show extensive genetic similarity under a model of speciation with ongoing gene flow (34), as observed between LEPC and GRPC when comparing levels of overall diversity (Figure S3), sequence divergence (0.003%; see Supporting Information), and differentiation (Table 1). Collectively, our results indicate there have been hybridization and introgression between the two species, consistent with behavioral observations (35, 36).

Our genomic analyses indicate that LEPC genes effectively flow among leks, ecoregions, and species as reflected by a lack of differentiation in the PCA (Figures 1b and S4), the admixture analyses (Figures 2 and S7), the analysis of variance (F_ST = 0.044), and the phylogeographic trees (Figure 4). We note that our mtDNA tree (Figure 4b) exhibits little genetic structure between LEPC sampled from the Northern and Southern DPSs or between species. This lack of structure reflected by the maternally inherited mtDNA tree may in part be the result of effective female-biased dispersal. Introgression between LEPC and GRPC has long been suspected (16, 17), and our genomic data confirm and quantify it (Figures 1b, 2b, S4, S7, and S9; (32)). GRPC and LEPC are estimated to have diverged from one another 0.6 – 0.9 Mya (13, 14). This relatively recent divergence has led to prairie grouse phylogenies that are fraught with uncertainty because of incomplete lineage sorting and further uncertainty exists because of secondary contact and introgression among lineages after the Pleistocene glaciers receded from the Great Plains (17, 37).

In the future, we expect that bidirectional introgression could continue or even expand in geographic scope and/or in the extent of effective hybridization. Environmental niche models indicate that over the last 130,000 years, there has been substantial overlap in potential habitat for each species (11). Furthermore, the area of sympatry harbors the largest contemporary population of LEPC (27, 38). If LEPCs expand habitat occupancy northward through management efforts like the Conservation Reserve Program (39), there may be increasing opportunities for future hybridization and introgression. Effective gene flow within and between species should continue to help maintain existing genomic variation and temper demographic bottlenecks that can otherwise linger and lead to genomic erosion (40).

Many closely related taxa (like the LEPC and GRPC) are formally recognized as two distinct species, each subject to different conservation strategies and formal legal frameworks. Unfortunately, the static Latin names of biological entities often do not accurately reflect nature. Incipient species can diverge despite genetic exchange via hybridization and introgression (5), and horizontal gene transfer between any two sexual species may have important policy ramifications. Specifically, the inability of whole-genome sequence data to consistently and clearly delineate individuals into species could raise questions about whether an individual is legally protected, or whether a population segment is in fact “discrete” as defined under the ESA (26). However, this depends entirely on how complex and ongoing evolutionary dynamics (such as speciation, hybridization, and introgression) are interpreted in the context of static conservation policies that vary at the state, federal, or international levels.

Rapid advancements in bioinformatics algorithms and next-generation sequencing technology (41) and the concomitant cost reduction are enabling unprecedented insight into evolutionary processes such as speciation and introgressive hybridization. While this is an exciting time for evolutionary biology, it also poses significant challenges to conservation, policy, and law (42). For example, there is the growing realization that individual eukaryotic organisms are often genomic composites of different species (including humans; (43)), so how does one legally enforce a policy when one cannot consistently define a species? The US ESA, which defines species as “any subspecies of fish or wildlife or plants, and any DPS of any species of vertebrate fish or wildlife which interbreeds when mature,” is now over 50 years old and was crafted decades before whole-genome sequencing. If current management practices are to represent the best available science and technology, a consensus is needed on how to best classify conservation units (e.g. which species concept to use). This is especially important when considering the management of populations that exhibit signals of recent and/or ongoing introgressive hybridization. Until legislation better reflects the emerging view of biological reality, species will remain the most important taxonomic unit for the conservation of biodiversity (44)—even when their boundaries are porous and ill-defined. The US ESA is a prime example whereby the lack of a consistent “hybrid policy” means that listing decisions are often handled on a case-by-case basis subject to little formal guidance (4). Increasingly, datasets such as ours indicate that hybridization and introgression are common in sexual organisms, and we think the dense information content contained in whole-genome sequences has the potential to be leveraged into more realistic, consistent, and objective policies where there is the political will to do so.

Materials and methods

Comprehensive methodological details and supplemental results (Tables and Figures S1-S12) are provided in the Supporting Information. Briefly, whole genomes from 433 prairie-chickens were sequenced and aligned to the LEPC reference genome (29). Genotype likelihoods were then leveraged to determine individual heterozygosity, genome-wide pairwise nucleotide diversity, ABBA/BABA tests of introgression, and genomic differentiation using Angsd (45). Runs of homozygosity were identified using Bcftools (46), and population stratification was assessed using Pcangsd (47) and Ngsadmix (48). A midpoint rooted nuclear genomic tree was created using Ngsdist (49), and a mtDNA tree was created using Iqtree (50).

Acknowledgments

The authors thank Dan Sullins, Samantha Robinson, John Kraft, Jonathan Lautenbach, Joseph Lautenbach, Reid Plumb, and Carleigh M. Stein for collecting/sharing biological material used in this study, Samarth Mathur for constructive feedback, Natalie Allen for assistance with estimating sequence divergence, and Megan Vhay for creation of Figure 1a. Bioinformatic processing was conducted on computing clusters housed at Purdue University's Rosen Center for Advanced Computing and Pittsburgh Supercomputing Center, which was made possible using ACCESS allocations (https://allocations.access-ci.org/).

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

This work was funded by the Western Association of Fish and Wildlife Agencies. J.A.D. was supported in part by the National Institute for Food and Agriculture. Funding for this work was derived from industry and government sources, but neither had any role in any aspect of study design, execution, interpretation, or writing/editing. Thus, the authors declare there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US government.

Author Contributions

A.N.B., A.J.M., J.Y.J., and J.A.D. conceived the ideas and designed methodology. D.H., K.J.B., K.A.F., A.G., and B.G. collected and/or processed the samples. A.N.B., A.J.M., and J.Y.J. analyzed the data. A.N.B. and J.A.D. led the writing of the manuscript. Z.E.L. and J.A.D. acquired funding and provided mentorship. All authors contributed critically to the drafts and gave final approval for publication.

Data Availability

Paired-end Illumina shotgun sequence data from 481 prairie-chicken samples has been deposited in the NCBI Sequence Read Archive under Bioproject PRJNA986511. All code, attribute data, and scripts have been uploaded to a publicly accessible repository (https://github.com/Andrew-N-Black/LEPC-popgen).

References

Author notes