Next-generation DNA sequencing (NGS) is used to study the genome sequences of Cryptosporidium spp., but NGS is challenging when pure Cryptosporidium oocysts are limited in number or not available. Varying levels of parasites present in fecal samples, combined with the abundance of host cells, bacterial and other microbial cells, and undigested food particles, often result in fecal DNA samples with ~0.1% Cryptosporidium DNA, making genome-scale sequencing of Cryptosporidium from such samples cost-prohibitive. DNA extractions from fecal samples are, however, widely available and commonly used for polymerase chain reaction (PCR)-based diagnostics which can detect fg levels of Cryptosporidium DNA in complex DNA mixtures. Here, we describe an Illumina NGS sample preparation protocol (iNextEra) that can generate libraries from a wide range of DNA input (<1 ng to >60 ng). We then use those libraries within a modified myBaits capture hybridization protocol using CryptoCap_100K baits to enrich Cryptosporidium genomic DNA from a complex DNA background to increase the percentage of generated sequence reads that map to target Cryptosporidium reference genome sequences. Thus, iNextEra libraries and capture hybridization facilitate genome-level sequencing of this critical pathogen from widely available samples with less cost, thereby opening new opportunities to understand the complex biology of this important pathogen.

Ticks are blood-feeding arthropods with approximately 1,000 species, however, only 24 species currently have a genome assembly. These genome assemblies are important resources to advance tick biology and control of tick-associated diseases. Generating tick genome assemblies is challenging due to their small body size (low DNA input), large genome size (approximately 2.4 Gbp for hard ticks), DNA contamination (from microbiota and host bloodmeals), and abundant transposable elements. Advances in sequencing technologies and assembly software have facilitated an increasing number of tick assemblies from 2011 to 2025. We characterize and assess the 54 tick genome assemblies within public genome databases using QUAST-LG and BUSCO compleasm. Then we evaluate the impact on these tick genome assemblies of biological source material, sequencing platform, and gene and repetitive element annotation. We identify 34 high-quality assemblies from 21 species that are well-suited for a variety of downstream analyses. Future tick genome assemblies should use long-read sequencing and Hi-C scaffolding to improve genomic resources for these unique blood-feeding parasites. ### Competing Interest Statement The authors have declared no competing interest. Centers of Disease Control and Prevention Pathogen Genomics Center of Excellence, Contract No. NU50CK000626 National Science Foundation, DGE-1545433 Dutch Ministry of Health, Welfare and Sport Howard Hughes Medical Institute, Life Sciences Research Foundation

Cryptosporidium spp. are protozoan parasites that cause severe illness in vulnerable human populations. Obtaining pure Cryptosporidium DNA from clinical and environmental samples is challenging because the oocysts shed in contaminated feces are limited in quantity, difficult to purify efficiently, may derive from multiple species, and yield limited DNA (<40 fg/oocyst). Here, we develop and validate a set of 100,000 RNA baits (CryptoCap_100k) based on six human-infecting Cryptosporidium spp. ( C. cuniculus , C. hominis , C. meleagridis , C. parvum , C. tyzzeri , and C. viatorum ) to enrich Cryptosporidium spp. DNA from a wide array of samples. We demonstrate that CryptoCap_100k increases the percentage of reads mapping to target Cryptosporidium references in a wide variety of scenarios, increasing the depth and breadth of genome coverage, facilitating increased accuracy of detecting and analyzing species within a given sample, while simultaneously decreasing costs, thereby opening new opportunities to understand the complex biology of these important pathogens.

The exponential growth of molecular sequence data over the past decade has enabled the construction of numerous clade-specific phylogenies encompassing hundreds or thousands of taxa. These independent studies often include overlapping data, presenting a unique opportunity to build macrophylogenies (phylogenies sampling >1000 taxa) for entire classes across the Tree of Life. However, the inference of large trees remains constrained by logistical, computational, and methodological challenges. The Avian Tree of Life provides an ideal model for evaluating strategies to robustly infer macrophylogenies from intersecting data sets derived from smaller studies. In this study, we leveraged a comprehensive resource of sequence capture data sets to evaluate the phylogenetic accuracy and computational costs of four methodological approaches: (1) supermatrix approaches using concatenation, including the “fast” maximum likelihood (ML) methods, (2) filtering data sets to reduce heterogeneity, (3) supertree estimation based on published phylogenomic trees, and (4) a “divide-and-conquer” strategy, wherein smaller ML trees were estimated and subsequently combined using a supertree approach. Additionally, we examined the impact of these methods on divergence time estimation using a data set that includes newly vetted fossil calibrations for the Avian Tree of Life. Our findings highlight the advantages of recently developed fast tree search approaches initiated with parsimony starting trees, which offer a reasonable compromise between computational efficiency and phylogenetic accuracy, facilitating inference of macrophylogenies.

Glacial cycles operating across Beringia have repeatedly exposed large swathes of the Bering Land Bridge, intermittently isolating and reuniting North American and Eurasian taxa. In high-latitude birds, these cycles are hypothesized to have been important in driving divergence and speciation. These repeated events have resulted in multiple trans-Beringian avian sister populations of varying degrees of taxonomic depth distributed across modern Beringia. We asked how these cyclic pulses have affected the temporal distribution and number of overall divergence events across Beringia. We sequenced full mitogenomes at high depth from 39 lineage pairs of varying levels of divergence, totaling 432 individuals of seven orders, 14 families, and 49 species from both Eurasia and North America. We then used a hierarchical approximate Bayesian comparative (hABC) approach to estimate the number and distribution of divergence events between the population pairs, using subsampled datasets. Net nucleotide divergence (DA) and Jukes-Cantor distance (JC-distance) were also calculated for each pairwise comparison to estimate divergence dates between taxa, using calibrated rates appropriate for shallow avian divergence events. Average divergence times were 200,000 ya for population-level taxa (n = 16), 720,000 ya for subspecies (n = 12), and 1 Mya for species (n = 11), although we consider these dating estimates conservative because of a lack of appropriate calibration for data of this quality. We found eighteen taxon pairs to be significantly differentiated (p < 0.05) by FST or substantially differentiated by haplotype clade, bounding the number of potential overall divergence events from 1 to 18, and two subsets of the full mitogenomic dataset analyzed in MTML-msBayes strongly supported simultaneous divergence of all Beringian lineages. However, this finding of simultaneous divergence is biologically unusual given the substantial variation in divergence dates among taxa and might indicate a relatively continuous spread of vicariance events, which is difficult to distinguish from a single, simultaneous vicariance event.

ABSTRACTTicks are obligate blood-feeding parasites associated with a huge diversity of diseases globally. The hard tickIxodes ricinusis the key vector of Lyme borreliosis and tick-borne encephalitis in Western Eurasia.Ixodesticks have large and repetitive genomes that are not yet well characterized. Here we generate two high-qualityI.ricinusgenome assemblies, with haploid genome sizes of approximately 2.15 Gbp. We find transposable elements comprise at least 69% of the twoI. ricinusgenomes, amongst the highest proportions found in animals. The transposable elements in ticks are highly diverse and novel, so we constructed a repeat library for ticks using ourI.ricinusgenomes and the genome ofI.scapularis, another major tick vector of Lyme borreliosis. To understand the impact of transposable elements on tick genomes we compared their accumulation in the twoIxodessister species. We find transposable elements in these two species to be drivers of genome evolution in ticks. TheI.ricinusgenome assemblies and our tick repeat library will be valuable resources for biological insights into this important ectoparasite. Our findings highlight that further research into the impact of transposable elements on the genomes of blood-feeding parasites is required.

Sarcoptic mange is a skin disease caused by the parasitic mite Sarcoptes scabiei that affects humans and over 150 species of domestic animals and wildlife worldwide. For the American black bear (Ursus americanus), sarcoptic mange is considered a significant emerging disease and case numbers have increased dramatically over the past four decades. While molecular techniques, such as DNA sequencing, are commonly used for diagnosis of clinical cases, they remain underutilized for studying cross-species transmission and regional mite variation. We captured complete mitochondrial genomes from S. scabiei mites (n = 17) collected from infested black bears in the Eastern United States (U.S.) using custom-designed oligonucleotide baits. These full mitogenomes were then used to refine primers targeting the cytochrome c oxidase subunit I (COI) gene. We then amplified and analyzed 285 COI sequences from S. scabiei collected across multiple North American wildlife hosts to assess their utility in inferring regional diversity and host-associated genetic structure. Full mitogenome comparisons with sequences from the U.S., Japan, and Australia revealed three major global S. scabiei clades, with two circulating among North American black bears. The COI haplotype networks showed regional clustering rather than host-specific patterns, suggesting that ecological or geographic factors influence transmission dynamics. Our findings confirm the presence of multiple S. scabiei lineages in American black bears and other U.S. wildlife and support the use of the COI gene for large-scale, cost-effective screening of genetic variation and shared haplotypes across host species. The integration of mitogenome and COI data offers a scalable model for investigating S. scabiei diversity in wildlife and domestic hosts.

Cryptosporidium parvum is a significant pathogen causing gastrointestinal infections in humans and animals. It is spread through ingesting contaminated food and water. Despite its global health significance, generating a C. parvum genome sequence has been challenging for many reasons including cloning and challenging subtelomeric regions. A new, gapless, hybrid, telomere-to-telomere genome assembly was created for C. parvum IOWA II, here termed CpBGF. It reveals 8 chromosomes, a genome size of 9,259,183 bp, and resolves complex subtelomeric regions. To facilitate ease of use and consistency with the literature, the chromosomes have been oriented, and genes in this annotation have been given similar gene IDs as those used in the 2004, C. parvum IOWA II reference genome sequence. The new annotation utilized considerable RNA expression evidence including single-molecule Iso-Seq data; thus, untranslated regions, long noncoding RNAs, and antisense RNAs are annotated. The CpBGF genome assembly serves as a valuable resource for understanding the biology, pathogenesis, and transmission of C. parvum, and it facilitates the development of diagnostics, drugs, and vaccines against cryptosporidiosis.

Once considered rare in eukaryotes, polycistronic mRNA expression has been identified in kinetoplastids and, more recently, green algae, red algae, and certain fungi. This study provides comprehensive evidence supporting the existence of polycistronic mRNA expression in the apicomplexan parasite Cryptosporidium parvum. Leveraging long-read RNA-seq data from different parasite strains and using multiple long-read technologies, we demonstrate the existence of defined polycistronic transcripts containing 2-4 protein encoding genes, several validated with RT-PCR. Some polycistrons exhibit differential expression profiles, usually involving the generation of internal monocistronic transcripts at different times during development. ATAC-seq in sporozoites reveals that polycistronic transcripts usually have a single open chromatin peak at their 5-prime ends, which contains a single E2F binding site motif. Polycistronic genes do not appear enriched for either male or female exclusive genes. This study elucidates a potentially complex layer of gene regulation with distinct chromatin accessibility akin to monocistronic transcripts. This is the first report of polycistronic transcription in an apicomplexan and expands our understanding of gene expression strategies in this medically important organism.

Ticks are a major health threat to humans and other animals, through direct damage, toxicoses, and transmission of pathogens. An estimated half a million people are treated annually in the United States for Lyme disease, a disease caused by the bite of a black-legged tick (Ixodes scapularis Say, 1821) infected with the bacterial pathogen Borrelia burgdorferi. This tick species also transmits another 6 human-disease causing pathogens, for which vaccines are currently unavailable. While I. scapularis are sexually dimorphic at the adult life stage, DNA sequence differences between male and female I. scapularis that could be used as a sex-specific marker have not yet been established. Here we identify sex-specific DNA sequences for I. scapularis (male heterogametic system with XY), using whole-genome resequencing and restriction site-associated DNA sequencing. Then we identify a male-specific marker that we use as the foundation of a molecular sex identification method (duplex PCR) to differentiate the sex of an I. scapularis tick. In addition, we provide evidence that this molecular sexing method can establish the mating status of adult females that have been mated and inseminated with male-determining sperm. Our molecular tool allows the characterization of mating and sex-specific biology for I. scapularis, a major pathogen vector, which is crucial for a better understanding of their biology and controlling tick populations.

The generation and maintenance of biodiversity are driven by population divergence and speciation. We investigated divergence, gene flow, and speciation in Beringia, a region at the top of the North Pacific Ocean with a history of dramatic landscape alteration through Pleistocene glacial cycles. These cycles repeatedly split and connected the Asian and North American continents, separating and reconnecting avian populations. Glacial refugia within Beringia also isolated some populations for a time before potentially enabling them to reunite during interglacial periods. Prior work suggests gene flow plays an important role in the divergence of Beringian birds. To improve our understanding of the generation of avian diversity in Beringia, we tested models of demographic history in 11 lineages from five avian orders (Anseriformes, Gaviiformes, Charadriiformes, Piciformes and Passeriformes) using population-, subspecies- and species-level pairwise comparisons. We sequenced an average of 3710 ultraconserved element (UCE) loci from the nuclear genomes of these taxa to examine genetic differentiation and test models of divergence through diffusion analysis for demographic inference (δaδi). All of the inferred best-fit models of divergence included gene flow. Together with prior work, this corroborates that divergence with gene flow is the predominant mode of divergence and speciation in Beringian birds.

Cryptosporidium is a globally endemic parasite genus with over 40 recognized species. While C. hominis and C. parvum are responsible for most human infections, human cases involving other species have also been reported. Furthermore, there is increasing evidence of simultaneous infections with multiple species. Therefore, we devised a new means to identify various species of Cryptosporidium in mixed infections by sequencing a 431 bp amplicon of the 18S rRNA gene encompassing two variable regions. Using the DADA2 pipeline, amplicons were first identified to a genus using the SILVA 132 reference database; then Cryptosporidium amplicons to a species using a custom database. This approach demonstrated sensitivity, successfully detecting and accurately identifying as little as 0.001 ng of C. parvum DNA in a complex stool background. Notably, we differentiated mixed infections and demonstrated the ability to identify potentially novel species of Cryptosporidium both in situ and in vitro. Using this method, we identified Cryptosporidium parvum in Egyptian rabbits with three samples showing minor mixed infections. By contrast, no mixed infections were detected in Egyptian children, who were primarily infected with C. hominis. Thus, this pipeline provides a sensitive tool for Cryptosporidium species-level identification, allowing for the detection and accurate identification of minor variants and mixed infections.IMPORTANCECryptosporidium is a eukaryotic parasite and a leading global cause of waterborne diarrhea, with over 40 recognized species infecting livestock, wildlife, and people. While we have effective tools for detecting Cryptosporidium in clinical and agricultural water samples, there is still a need for a method that can efficiently identify known species as well as infections with multiple Cryptosporidium species, which are increasingly being reported. In this study, we utilized sequencing of a specific region to develop a sensitive and accurate identification workflow for Cryptosporidium species based on high-throughput sequencing. This method can distinguish between all 40 recognized species and accurately detect mixed infections. Our approach provides a sensitive and reliable means to identify Cryptosporidium species in complex clinical and agricultural samples. This has important implications for clinical diagnostics, biosurveillance, and understanding disease transmission, ultimately benefiting clinicians and produce growers.

Multiple displacement amplification (MDA) outperforms conventional PCR in long fragment and whole-genome amplification, making it attractive to couple MDA with long-read sequencing of samples with limited quantities of DNA to obtain improved genome assemblies. Here, we explore the efficacy and limits of MDA for efficient low-cost genome sequence assembly using Oxford Nanopore Technologies (ONTs) rapid library preparations and minION sequencing. We successfully generated almost complete genome sequences for all organisms examined, including Gram-positive (Staphylococcus aureus, Enterococcus faecium) and Gram-negative (Escherichia coli) prokaryotes and one challenging eukaryotic pathogen (Cryptosporidium spp) representing a broad spectrum of critical infectious disease pathogens. High-quality data from those samples were generated starting with only 0.025 ng of total DNA. Controlled sheared DNA samples exhibited a distinct pattern of size increase after MDA, which may be associated with the amplification of long, low-abundance fragments present in the assay, as well as generating concatemeric sequences during amplification. To address concatemers, we developed a computational pipeline (CADECT: Concatemer Detection Tool) to identify and remove putative concatemeric sequences. This study highlights the efficacy of MDA in generating high-quality genome assemblies from limited amounts of input DNA. Also, the CADECT pipeline effectively mitigated the impact of concatemeric sequences, enabling the assembly of contiguous sequences even in cases where the input genomic DNA was degraded. These results have significant implications for the study of organisms that are challenging to culture in vitro, such as Cryptosporidium, and for expediting critical results in clinical settings with limited quantities of available genomic DNA.

Local adaptation occurs when populations evolve traits in response to local environmental challenges. Isolated island populations often experience different selection pressures than their mainland counterparts, which enables the study of how phenotypes and genotypes respond to differing selection regimes. We studied a group of five phenotypically differentiated subspecies of song sparrow (Melospiza melodia) in Alaska that demonstrate striking body size, color, and migratory behavioral differences to examine the effects of local adaptation on phenotypes and genotypes. We examined the phenotypic attributes of these populations and used whole-genome data to determine relationships and test candidate loci for evidence of selection. Phenotypic measurements of museum specimens (n = 227) quantified the dramatic size differences among these populations, with westernmost M. m. maxima being ~1.6 times larger than easternmost M. m. rufina. Using ultraconserved elements (UCEs) and McDonald-Kreitman tests, we showed that seven candidate genes associated with bill size, circadian rhythm regulation, plumage color, and salt tolerance exhibited signs of putative positive selection. Phylogenetic analysis of UCEs identified M. m. maxima as sister to the other Alaska M. melodia subspecies. This suggests M. m. maxima colonized earliest, perhaps before the last glacial maximum, and that Alaska was later recolonized by ancestors of the remaining four subspecies.

The trunks of forest trees store massive amounts of carbon, but fungi actively and invisibly decay wood inside even seemingly healthy trees. Wood-decay fungi are responsible for the loss of stored carbon in living trees, and they make trees susceptible to snapping and uprooting in storms. We used sonic tomography to measure the prevalence and severity of decay in 1744 live trees (≥20 cm diameter) of 171 species on the 50-ha Forest Dynamics Plot on Barro Colorado Island, Panama. A median of <2% of the cross-sectional trunk area showed decay, but 15% of trees had >20% decay. Twenty percent of the combined basal area showed decay, representing a loss of approximately 1% of aboveground biomass. Larger trees more often showed internal decay, with one quarter of trees showing decay before reaching canopy height. Decay severity varied by species; 23% of species showed <2% decay while 9% of species lost over half their basal area. Rare species were more affected than locally abundant species, and species with traits associated with a fast life history were more susceptible to decay. These results suggest that hidden wood decay affects a large proportion of living tropical forest trees.

The evolutionary histories of different genomic regions typically differ from each other and from the underlying species phylogeny. This makes species tree estimation challenging. Here, we examine the performance of phylogenomic methods using a well-resolved phylogeny that nevertheless contains many difficult nodes, the species tree of living birds. We compared trees generated by maximum likelihood (ML) analysis of concatenated data, gene tree summary methods, and SVDquartets. We also conduct the first empirical test of a “new” method called METAL (Metric algorithm for Estimation of Trees based on Aggregation of Loci), which is based on evolutionary distances calculated using concatenated data. We conducted this test using a novel dataset comprising more than 4,000 ultraconserved element (UCE) loci from almost all bird families and two existing UCE and intron datasets sampled from almost all avian orders. We identified “reliable clades” very likely to be present in the true avian species tree and used them to assess method performance. ML analyses of concatenated data recovered almost all reliable clades with less data and greater robustness to missing data than other methods. METAL recovered many reliable clades, but only performed well with the largest datasets. Gene tree summary methods (weighted ASTRAL and weighted ASTRID) performed well; they required less data than METAL but more data than ML concatenation. SVDquartets exhibited the worst performance of the methods tested. In addition to the methodological insights, this study provides a novel estimate of avian phylogeny with almost 99% of the currently recognized avian families. Only one of the 181 reliable clades we examined was consistently resolved differently by ML concatenation versus other methods, suggesting that it may be possible to achieve consensus on the deep phylogeny of extant birds.

GITDs are among the most common causes of death in adult and young horses in the United States (US). Previous studies have indicated a connection between GITDs and the equine gut microbiome. However, the low taxonomic resolution of the current microbiome sequencing methods has hampered the identification of specific bacterial changes associated with GITDs in horses. Here, we have compared TEHC, a new approach for 16S rRNA gene selection and sequencing, with conventional 16S rRNA gene amplicon sequencing for the characterization of the equine fecal microbiome. Both sequencing approaches were used to determine the fecal microbiome of four adult horses and one commercial mock microbiome. Our results show that TEHC yielded significantly more operational taxonomic units (OTUs) than conventional 16S amplicon sequencing when the same number of reads were used in the analysis. This translated into a deeper and more accurate characterization of the fecal microbiome when the samples were sequenced with TEHC according to the relative abundance analysis. Alpha and beta diversity metrics corroborated these findings and demonstrated that the microbiome of the fecal samples was significantly richer when sequenced with TEHC compared to 16S amplicon sequencing. Altogether, our study suggests that the TEHC strategy provides a more extensive characterization of the fecal microbiome of horses than the current alternative based on the PCR amplification of a portion of the 16S rRNA gene.

During the COVID-19 pandemic, the detection and sequencing of SARS-CoV-2 from wastewater proved to be a valuable tool in assessing trends at the community level. Several whole genome enrichment methods have been proposed for sequencing SARS-CoV-2 from the mixed wastewater community, but there is little consensus on the most appropriate sequencing methods for variant detection or abundance estimations. Few studies have elucidated the errors associated with these methods or have established minimum sequencing requirements for correct interpretation of the results. To address these needs, we systematically assessed the efficacy of three tiled amplicon enrichment methods (Freed/Midnight, ARTIC V4, NEB VarSkip) for whole genome sequencing of SARS-CoV-2 variants using mock wastewater communities with variants at known proportions. We found the ARTIC V4 approach yielded the most accurate results for variant identification and variant abundance estimation, followed by the NEB VarSkip approach. Conversely, the NEB VarSkip method obtained the highest genomic coverage, with the ARTIC V4 method achieving the second highest coverage. Finally, we determined that the Freed/Midnight library preparation methods are not well-suited for use with short read sequencing. Based on the present results, the ARTIC V4 workflow appears to be the most robust and cost-effective approach for monitoring circulating SARS-CoV-2 variants with wastewater surveillance. IMPORTANCE This work is informative for practitioners of wastewater-based epidemiology. Here, we detail a systematic comparison of three tiled amplicon sequencing approaches for enrichment of SARS-CoV-2 variants from wastewater. Using mock communities of known variant composition, we validate the analysis methods previously published by Baaijens et al. in Genome Biology (2022) for estimating variant abundance from wastewater using an RNAseq pipeline, kallisto. We provide recommendations for minimum sequencing requirements for accurate abundance estimates of SARS-CoV-2 variants in wastewater. The sequences generated from the mock communities have been uploaded to NCBI’s Sequence Read Archive and will be useful to other practitioners seeking to validate their sequencing methods or bioinformatic pipelines. ### Competing Interest Statement The authors have declared no competing interest.

SUMMARY Public health microbiology focuses on microorganisms and infectious agents that impact human health. For years, this field has relied on culture or molecular methods to investigate complex samples of public health importance. However, with the increase in accuracy and decrease in sequencing cost over the last decade, there has been a transition to the use of next-generation sequencing in public health microbiology. Nevertheless, many available sequencing methods (e.g., shotgun metagenomics and amplicon sequencing) do not work well in complex sample types, require deep sequencing, or have inherent biases associated with them. Hybridization bait capture, also known as target enrichment, brings in solutions for such limitations. It is an increasingly popular technique to simultaneously characterize many thousands of genetic elements while reducing the amount of sequencing needed (thereby reducing the sequencing costs). Here, we summarize the concept of hybridization bait capture for public health, reviewing a total of 35 bait sets designed in six key topic areas for public health microbiology [i.e., antimicrobial resistance (AMR), bacteria, fungi, parasites, vectors, and viruses], and compare hybridization bait capture to previously relied upon methods. Furthermore, we provide an in-depth comparison of the three most popular bait sets designed for AMR by evaluating each of them against three major AMR databases: Comprehensive Antibiotic Resistance Database, Microbial Ecology Group Antimicrobial Resistance Database, and Pathogenicity Island Database. Thus, this article provides a review of hybridization bait capture for public health microbiologists.

The origin and eventual loss of biogeographic barriers can create alternating periods of allopatry and secondary contact, facilitating gene flow among distinct metapopulations and generating reticulate evolutionary histories that are not adequately described by a bifurcating evolutionary tree. One such example may exist in the two-lined salamander (Eurycea bislineata) species complex, where discordance among morphological and molecular datasets has created a "vexing taxonomic challenge." Previous phylogeographic analyses of mitochondrial DNA (mtDNA) suggested that the reorganization of Miocene paleodrainages drove vicariance and dispersal, but the inherent limitations of a single-locus dataset precluded the evaluation of subsequent gene flow. Here, we generate triple-enzyme restriction site-associated DNA sequencing (3RAD) data for > 100 individuals representing all major mtDNA lineages and use a suite of complementary methods to demonstrate that discordance among earlier datasets is best explained by a reticulate evolutionary history influenced by river drainage reorganization. Systematics of such groups should acknowledge these complex histories and relationships that are not strictly hierarchical.