Sustaining genetic diversity in an imperiled pelagophilic fish despite genetic drift

Sustaining genetic diversity in an imperiled pelagophilic fish despite genetic drift | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Sustaining genetic diversity in an imperiled pelagophilic fish despite genetic drift Guilherme Caeiro-Dias, Alexander Cameron, Thomas Turner, Megan Osborne This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7863499/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Patterns of genetic variation within species may differ in space and time, particularly in the face of rapid anthropogenic changes. Genomic time-series population-level data can reveal patterns of allele frequency change and can help understanding on-going evolutionary processes. This study used a time-series spanning 6 years to evaluate genomic variability and trends in effective population size, and to assess the influence of genetic drift on Macrhybopsis tetranema genomic diversity. Temporal trends were evaluated from samples collected from 2015 to 2020. Results from genomic variation and variance effective population size across time suggests M. tetranema upstream (New Mexico) are affected by genetic drift, but genetic diversity did not change significantly across the time-series. Absolute allele frequency differences across time-series revealed that drift is a result of many loci with moderate allele frequency changes, particularly after 2018, rather than few loci with high frequency changes. Also, linkage disequilibrium effective population size showed fluctuations and a general increasing trend. Overall, the results presented here demonstrated that genetic drift is an overarching evolutionary force driving trends in genetic variation in M. tetranema . Together, genetic data reported here and ecological data from recent studies, suggest that M. tetranema in the upstream portion of the remnant distribution experiences genetic drift likely due of downstream-biased gene flow, while upstream movement of adults can explain low relatedness and maintenance of genetic diversity upstream in the face of genetic drift. genomic time-series neutral genetic diversity peppered chub Macrhybopsis tetranema South Canadian River Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Introduction Contemporary patterns of intraspecific genetic variation are the result of demographic and evolutionary processes. Identifying specific processes generating these patterns is critically important to ecological, evolutionary and conservation studies. In particular, contemporary demographic and microevolutionary processes may impact recent patterns of genetic variation (Messer et al. 2016 ). Time-series population-level genetic data can reveal patterns of allele frequency change over contemporary time scales and can be important for understanding recent evolutionary processes, including short term maintenance (Osborne et al. 2012 ) and loss (Osborne et al. 2023a ) of genetic diversity in the face of wild population augmentation, differences in genetic declines in distinct subpopulations in a metapopulation (Mathieu-Bégné et al. 2019 ), or genetic drift after isolation (Demandt 2010 ). Although genome-wide time-series data are slowly accumulating, they are still uncommon in the literature. Nonetheless, such studies have uncovered evidence of rapid adaptation (Bergland et al. 2014 ), climate-mediated hybrid zone movement (Ryan et al. 2018 ), maintenance of high genetic diversity due to gene flow despite genetic drift (Gompert et al. 2021 ), and genomic signatures of genetic drift and increased inbreeding in the face of massive population declines despite population supplementation efforts (Osborne et al. 2023b ). Patterns of genetic variation within species may also differ in space. In lotic freshwater ecosystems with dendritic connectivity (e.g., rivers, streams) several empirical studies found an increase of genetic diversity downstream (Alp et al. 2012 ; Torterotot et al. 2014 ; Paz-Vinas et al. 2018 ). While upstream-directed colonization can replenish diversity at upstream sites, the distribution of genetic diversity across the riverscape also depends on species’ intrinsic traits such as fecundity, reproductive strategy, and migratory behavior (Morrissey and de Kerckhove 2009 ). Anthropogenic changes can further impact patterns of genetic diversity in space and time in riverine systems. Environmental alterations induced by anthropogenic activities such as habitat changes, climate change, and/or introduction of non-native species, together with stochastic events have disrupted freshwater fish population connectivity, altered population dynamics (Kuczynski et al. 2018 ; Chavez et al. 2024 ; Franklin et al. 2024 ), and caused substantial changes in abundance (Marschall and Crowder 1996 ; Yackulic et al. 2022 ), which in turn can affect levels of genetic diversity (DiBattista 2008 ) and its distribution across the riverscape (Blanchet et al. 2020 ). Across the Great Plains of North America, fishes belonging to the pelagic broadcast spawning reproductive guild have declined in range and abundance (Perkin et al. 2015 ). Pelagophils release neutrally buoyant eggs that are fertilized once released (Platania and Altenbach 1998 ) and drift while developing (Bottrell et al. 1964 ). Macrhybopsis tetranema (peppered chub), a pelagic spawning minnow, was historically distributed throughout the Arkansas River basin in New Mexico (NM), Kansas, Texas (TX), and Oklahoma (Eisenhour 1999 ). This species is now extirpated from the majority of its historical range. Until about a decade ago, three isolated populations remained but a drought cycle in the late 1980’s and early 1990’s resulted in decline and eventual extirpation (last seen in 2011) from the Cimarron River (U.S. Fish and Wildlife Service 2018 ) and further drought between 2011 and 2013 resulted in extirpation of the Ninnescah and Arkansas river population in Kansas (Perkin et al. 2015 ; Pennock et al. 2017 ). Today a single extant population persists in a 218 km stretch of the South Canadian River between Ute Lake (NM) and Lake Meredith (TX) (Bonner and Wilde 2000 ; Pennock et al. 2017 ). As such, M. tetranema was listed as federally endangered (U.S. Fish and Wildlife Service 2022). Macrhybopsis tetranema evolved life-history traits that facilitated its persistence in rivers with marked seasonality and a highly variable flow regime. Traits include small body size, short generation time (1–3 years; Wilde and Durham 2008 ), and non-adhesive semi-buoyant eggs characteristic of the pelagophilic reproductive guild of fishes (Balon 1975 , 1981 ). Because eggs drift while developing, persistence of upstream populations must depend on the retention of eggs and larvae in local nursery habitats or upstream dispersal. Upstream movement was observed in post-larval life stages of several pelagophilic species (Archdeacon et al. 2018 ; Platania et al. 2020 ) including young-of-the-year (Stuart and Sharpe 2020 ). River fragmentation by human-made barriers (e.g., impoundments, dams) disrupt upstream dispersal. These barriers and the habitat changes they cause, contribute to declines of pelagophilic minnows including M. tetranema (Luttrell et al. 1999 ; Dudley and Platania 2007 ; Perkin and Gido 2011 ; Perkin et al. 2019 ; Archdeacon et al. 2020 ). River fragmentation impacts habitat availability and quality, reduces stream discharge, increases frequency of stream dewatering, alters flow periodicity, and increases stream channelization (Hoagstrom and Turner 2015). Together these factors disrupt the reproductive cycle, decrease available nursery habitat, and enhance probabilities of recruitment failure in pelagic broadcast spawners (Dudley and Platania 2007 ; Archdeacon et al. 2020 )d tetranema is no exception (Wilde and Durham 2008 ; Perkin et al. 2019 ). For a short-lived minnow species, recruitment failure contributes to population decline and potential reduction of genetic diversity. The combination of fragmented habitats, limited geographic distribution, short lifespan and recruitment failure increase the chance of local population extirpation (Perkin et al. 2015 ; Pennock et al. 2017 ) and, in the worst case, can lead to the extinction of this species. Using microsatellite and mtDNA data collected between 2015 and 2018, Osborne et al. ( 2021 ) found that the remnant M. tetranema population maintained genetic diversity across the time-series and that effective population size (N e ) estimates showed an increasing trend consistent with the demographic trajectory. We re-evaluated putatively neutral genetic variability and effective size (N e ) trends using genome-wide time-series data from 2015 to 2020. Specifically, we (1) evaluated the magnitude of genetic drift from temporal sampling, (2) tested if there was a significant loss of genetic diversity across the extended time-series, and (3) assessed temporal trends of N e . Results from genome-wide data provided a more comprehensive understanding of recent changes of neutral genetic diversity in M. tetranema but also raised new questions about population dynamics in the remnant range in the South Canadian River. Material and Methods Sampling Temporal sampling for reduced-representation sequencing included archived DNA isolates from 2015 (n = 34), 2016 (n = 31), 2017 (n = 21), and 2018 (n = 9), and new samples from 2019 (n = 41) and 2020 (n = 43) collected from multiple localities across the South Canadian River in New Mexico (NM) from Ute Late to NM-TX border (Fig. 1 ). This region is referend in this study as “upstream”. Because it was not possible to obtain recent collections from Texas, analysis from downstream of the distribution range was restricted to relatively small number of archived samples collected in 2017 (n = 10) from a single site (Fig. 1 ). For the draft genome sequencing, a M. tetranema individual was collected from the South Canadian River, NM (approximate geographic coordinates: 35.389376, − 103.355331) in August 2021. Species identity was determined by visual examination of external morphology by experienced personnel from US Fish and Wildlife Service. The fish was euthanized with an overdose of MS222. Remaining tissue and the DNA isolate were deposited in the Museum of Southwestern Biology (MSB Cat no. 117041 for M. tetranema ; https://msb.unm.edu/ ). A complete mitochondrial genome has previously been reported for this specimen (Osborne et al. 2023b ). Draft genome sequencing and assembly High-molecular-weight genomic DNA was isolated from muscle and fin tissue using genomic-tips 20/G (Qiagen®) according to the manufacturer’s directions. Prior to sequencing, DNA quality was assessed using Qubit HS assays (Invitrogen®, Thermo Fisher Scientific) and using pulse-field electrophoresis on a 1% agarose gel. A genomic library was prepared using PacBio HiFi SMRTbell® and the whole genome was sequenced using 1 SMRT cell on a PacBio Sequel II platform at the UC Davis Genomics Core facility. To estimate M. tetranema genome size, Jellyfish v. 2.2.7 (Marçais and Kingsford 2011 ) was used to obtain the k-mers count which was then used as input to GenomeScope v. 2.0 online tool ( http://genomescope.org/genomescope2.0/ ; Ranallo-Benavidez et al. 2020 ). The genome was assembled from PacBio reads with the HiFiAsm assembler v.0.18.1-r466 (Cheng et al. 2021 ) using the default parameters. Genome assembly was improved by scaffolding the contigs assembled by HiFiAsm to scaffold-level assemblies available for other leuciscid fishes using the homology-based scaffolding software RagTag (Alonge et al. 2022 ). Briefly, assembled contigs from M. tetranema were mapped to primary genomic assemblies for Spikedace and Loach Minnow ( Meda fulgida and Tiaroga cobitis ; Alexandre et al., 2023 ) and the Eurasian Minnow ( Phoxinus phoxinus , Nunn et al., 2024 ). The default aligner minimap2 (Li 2018 ) was used and mapping quality score was increased to 20. Resulting assembly gap path (AGP) files were then merged with the edge weight set to 2 (e.g., scaffold joins had to be supported by at least 2 AGPs). The online platform CoGe ( https://genomevolution.org/coge/ ) was used to obtain summary statistics from both contig and scaffold-level assemblies. NextRAD sequencing, variant call, quality filtering and microhaplotype data Isolation of DNA from tissue samples (2019 and 2020 collections) was performed using E.Z.N.A. tissue DNA (Omega Bio-Tek Inc.) or Zymo Quick-DNA (Zymo Research Corp.) kits according to the procedures outlined by the manufacturer. Isolations of DNA from previous genetic monitoring (2015–2018; Osborne et al. 2021 ) were purified using Zymo DNA clean and concentrator kits (Zymo Research Corp.) to remove phenol which can inhibit DNA sequencing. All samples were treated with RNase A following DNA extraction. Each isolate was evaluated for the presence of high molecular weight DNA and the absence of RNA contamination by electrophoresis on a 1.2% agarose gel and by quantification using Qubit HS assays. One-hundred and ninety samples with sufficient high-quality DNA were sent to SNPsaurus, LLC (University of Oregon) for library preparation and sequencing. Genomic DNA was converted into a Nextera-tagmented reductively-amplified DNA (nextRAD) library and sequenced according to Russello et al. ( 2015 ). Genomic DNA was initially fragmented with Nextera DNA Flex reagent (Illumina®, Inc), which also ligates short adapter sequences to the ends of the fragments. The Nextera reaction was scaled for fragmenting 30 ng of genomic DNA, although 60 ng of genomic DNA was used for input to compensate for degraded DNA and to increase fragment sizes in samples. Fragmented DNA was then amplified for 27 cycles at 74 degrees, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence GTGTAGAGCC. Only fragments starting with a sequence that can be hybridized by the selective sequence of the primer will be efficiently amplified. All samples were pooled for the nextRAD library and sequenced with 150 base pair (bp) single-end reads on one Illumina Hi-Seq 4000 lane. Raw sequence reads were received from the sequencing provider, trimmed for adapter sequences, and demultiplexed by individual and by lane. Raw reads were further trimmed to remove low quality bases with Trimmomatic v. 0.39 (Bolger et al. 2014 ). Bases on both extremes of each read were removed if quality was bellow 20 or if ambiguous (N). Each read was then scanned with a 5-base wide sliding window, cutting when the average quality per base drops below 10. After trimming, reads were discarded if smaller than 60 bases. Retained reads were mapped against the scaffold level draft genome with Bowtie v. 2.4.2 (Langmead and Salzberg 2012 ) using the ‘local alignment’ and default ‘very sensitive’ options. Successfully aligned reads were filtered with SAMtools v. 1.16 (Li et al. 2009 ; Danecek et al. 2021 ) to remove reads with mapping quality lower than 20. Before variant calling, we used Picard tools v. 2.26.2 (Broad Institute 2019) to add read group (RG) flags to bam files. Genetic variants were identified using FreeBayes v. 1.3.6 (Garrison and Marth 2012 ). FreeBayes uses base quality scores to estimate a probability for each allele. We kept a maximum of 10 raw variants from each alignment with higher probabilities and at least a base quality of five. To remove erroneous or potentially erroneous variants, we used VCFtools v. 0.1.16 (Danecek et al. 2011 ) to filter out variants with mean depth of coverage lower than 20 and higher than 200, minor allele count less than three, minor allele frequency lower than 2%, genotype depth of coverage lower than five, and with genotype quality lower than 20. Multi-nucleotide states were decomposed into single variants with vcflib ( https://github.com/ekg/vcflib ) and VCFtools was used to filter out nucleotide insertions and deletions and to retain only the bi-allelic SNPs. The dataset was then filtered by missing data, keeping SNPs present in at least 80% of samples and removing individuals with more than 30% missing data. After this step, SNPs were filtered using the bash script dDocent_filters ( https://github.com/jpuritz/dDocent/blob/master/scripts/dDocent_filters ) that uses vcflib and VCFtools to filter loci based on allelic balance at heterozygous genotypes, strand representation, and quality vs depth of coverage. First, loci were removed if at heterozygous positions, the alternative allele had a coverage lower than 20% or higher than 80% compared with the reference allele as reads with alleles from heterozygous positions are expected to have similar frequencies in the same individual. Alleles with frequencies smaller than 0.01 and higher that 0.99 were not removed to account for fixed alleles. Additionally, if the quality sum of the reference or alternative allele was zero, the locus was removed. This removes positions with spurious heterozygous genotype calls. Then loci with the ratio between the mean mapping quality of the alternative and reference allele lower than 0.25 or higher than 1.75 were removed, because loci from the same genomic location should not have large discrepancy between mapping qualities of two alleles. Furthermore, loci with quality scores less than half of the total depth were excluded because excessive depth inflates FreeBayes quality scores. Of the remaining loci, the average depth and standard deviation across all individuals was calculated. Loci with depth greater than the average depth plus one standard deviation were removed if the quality score was less than two times the depth. Finally, this script removed loci with a mean depth across individuals greater than two times the mode (98) that corresponded approximately to the 95th percentile of mean depth. Subsequently potential erroneous SNPs were filtered based on Hardy-Weinberg equilibrium (HWE) expectations with the pearl script filter_hwe_by_pop.pl ( https://github.com/jpuritz/dDocent/blob/master/scripts/filter_hwe_by_pop.pl ). Typically, errors would have a low p-value and would be present in many populations; SNPs present in more that 50% of the populations (here each year was considered a ‘population’) and with an HWE p-value lower than 0.001 were removed. We filtered out potential incorrectly assembled paralogous loci that exhibited a large variation in read depth across all individuals. Standard deviation was estimated with package stats implemented in R v. 4.2.1 (R Core Team 2022) and read depth with VCFtools. Additional filtering based on missing data per locus (keeping loci present in 80% of individuals) was then reapplied. Remaining SNPs were used to identify haplotypes within loci (referred to as microhaplotypes). Haplotyping SNPs within a locus also eliminates possible paralogous loci while neutralizing physical linkage without losing data (Willis et al. 2017 ). This was performed with the rad_haplotyper.pl pearl script ( https://github.com/chollenbeck/rad_haplotyper ) excluding microhaplotypes if considered paralogs at least in five individuals and if missing from more than 30% of individuals. Retained loci were tested for deviations from HWE and for linkage disequilibrium (LD) considering individuals captured in each year as a single ‘population’. Departures from HWE were accessed using a chi-square test on microhaplotype data with R package pegas v. 1.0 (Paradis 2010 ) and using the Bonferroni correction for multiple comparisons implemented in the R package rcompanion v. 2.4.0 (Mangiafico 2021 ), as implemented in the R function multi_HWE_tests ( https://github.com/gcaeirodias/multi HWE tests; Caeiro-Dias et al. 2024 ). Estimations of LD were performed on SNP data using the SNP of each microhaplotype with higher minimum allele frequency. If a SNP was found in LD, then the entire locus was removed. Tests for LD were performed using the chi-square test implemented in the R package GUSLD v. 1.0.1 (Bilton et al. 2018 ) and the Bonferroni correction to account for multiple simultaneous tests as implemented in the R pipeline significantLD ( https://github.com/gcaeirodias/significantLD ; Caeiro-Dias et al. 2024 ). If loci were found in multiple significant LD pairs, the loci that appeared in the highest number of comparisons were discarded to keep the maximum number of loci possible. In the remainder of instances, one locus from each pair were discarded randomly. Loci were considered as deviating from HWE and to be in LD if tests were significant across the six temporal samples (p-value < 0.05). The resulting dataset should represent a robust genome-wide neutral SNP dataset. After filtering, samples from one collection (TX 2017) had increased missing data when compared to other collections (see subsection NextRAD sequencing, variant call, quality filtering and microhaplotype data of Results). Missing data can bias some results. As such, two datasets were created from the filtered set of markers and individuals described before. One dataset included all the microhaplotypes, but the TX 2017 collection was removed (NM dataset). For the other dataset, we removed the loci with more than 20% missing data in TX 2017 collection while retaining all the individuals after the filtering described above (NM-TX dataset). Unless otherwise specified, the analyses described below were performed with both datasets. Temporal genetic variation and evaluation of genetic drift Temporal genetic variation was first visualized using a discriminant analysis of principal components (DAPC), which summarizes genotypes in principal components to construct linear functions that maximize among group variation while minimizing within group variation. Analysis was performed using the R package adegenet v. 1.3–1 (Jombart 2008 ; Jombart and Ahmed 2011 ). Prior to DAPC we replaced missing data within each temporal sample using the Breiman's regression random forest algorithm (Breiman 2001 ) implemented in R package randomForest v. 4.6–14 (Liaw and Wiener 2002 ). Values of missing data were predicted from 500 independently constructed regression trees and 50 bootstrap iterations with default bootstrap sample size. This was preferred over the default “mean method” (i.e., missing genotypes are replaced by the average estimated across the data set) implemented in adegenet to ensure that we did not artificially increase similarity of allele frequencies across temporal samples. An initial DAPC was performed using temporal samples as groups, allele frequencies centered but not scaled, retaining all PCA and DA axes, and keeping other options as default. The a-score method was used to select the optimal number of principal components to retain for the final DAPC, using the maximum number of PCs, all DAs and the other options were set as default. The final DAPC was performed using the optimal number of PCs, two DAs and keeping the other default options. We assessed genetic differentiation over time by estimating pairwise F ST between all temporal collections and p-values using 1,000 bootstrap iterations over loci as implemented in GenoDive v. 3.06 (Meirmans 2020 ). Significant values of F ST between consecutive years or those separated by a few years, reflect differences in allele frequencies and can be a sign of genetic drift. Next, we directly quantified allele frequency changes in each locus across time by calculating absolute allele frequency differences (AFD) between consecutive years and identified outliers of frequency changes in each pair of years. This analysis was conducted to evaluate the proportion of the loci with inflated changes in allele frequency. When using the NM-TX dataset we did not include TX 2017 because there were no consecutive years to estimate allele frequency differences for this location. Allele frequencies were estimated with GenoDive and for each locus, allele frequency differences were averaged across alleles. Temporal genetic diversity, inbreeding, and relatedness For each temporal collection, we estimated standard genetic diversity and inbreeding metrics. We used R package diverRsity v. 1.9.90 (Keenan et al. 2013 ) to estimate expected heterozygosity (HS; Nei 1987 ), observed heterozygosity (H O ), allelic richness (A R ), and inbreeding coefficient (F IS ). All metrics were estimated for each population and by locus within each population. If changes to genomic diversity occurred across time, we expect to see significant changes in the distribution of values of those metrics across temporal collections. To evaluate this, we tested if the distribution of values for each metric estimated with both datasets were similar across populations using non-parametric Kruskal-Wallis tests, because diversity measures did not follow a normal distribution and/or variance across populations was not homogeneous (data not shown). When tests were significant, we performed pairwise Wilcoxon tests between populations using a Bonferroni correction for multiple comparisons. All tests were performed with the built-in R package stats. Maximum likelihood estimation (MLE) of identity-by-descent (IBD) for each pair of individuals was estimated with the R package SNPrelate v. 1.32.2 (Zheng et al. 2012 ) using the Expectation-Maximization algorithm. Because SNPrelate does not take loci with more than two alleles as input, we selected one SNP at random for each locus. SNPs with more than 5% missing data were excluded as it can have an impact on relatedness analyses. Resulting MLE-IBD coefficients (k0 and k1) were used to construct k0,k1 plots (e.g., Galván-Femenía et al. 2017 ), where k0 and k1 are the probabilities that two individuals share zero or one IBD alleles, respectively. This analysis was performed to evaluate if the overall relatedness changed overtime. Temporal trends in effective population size Linkage disequilibrium effective size (LD N e ) was estimated from both datasets for every year using the method from Hill ( 1981 ) implemented in NeEstimator v. 2.1 (Do et al. 2014 ) excluding microhaplotypes with allele frequencies lower than 2%. This metric reflects the effective number of parents that produced the progeny from which the sample was drawn (see Waples 2005 ). To estimate the 95% confidence intervals (CIs) the jackknife approach may be more suitable for data sets with large numbers of loci (Jones et al. 2016 ) as loci may not be entirely independent (Gilbert and Whitlock 2015 ). However, datasets in this study were filtered to avoid linkage and indeed linkage among loci was limited as r 2 per population with more than 20 samples was always lower than 0.06. Hence, a parametric approach was used to estimate the 95% CIs. For consecutive temporal samples, variance effective population size (N eV ) and 95% CIs were estimated using the temporal method of Nei and Tajima ( 1981 ) as implemented in NeEstimator v. 2.1, using the Plan I option (sampling with replacement), excluding microhaplotyes with allele frequencies lower than 2%, as highly polymorphic loci with many rare alleles can result in biased estimates of N eV (Hedrick 1999; Turner et al. 2001). The parametric approach was used to estimate the 95% CIs. To estimate N eV we also used both datasets, but TX 2017 was excluded from NM-TX dataset as we only had a single temporal sample from TX. Results Draft genome sequencing A total of 2,244,201 raw PacBio reads were generated, containing 29.72 Giga bases (Gb) of sequence data. Genome size was estimated from the k-mers count to be 1.04 Gb. Contig assembly resulted in 960 sequences, ranging from 6,390 to 38,396,035 bp (mean = 1,331,581, mode = 30,561) with a total length of 1,259,117,741 bp (1.26 Gb) and an N50 of 18,070,155 bp. After scaffolding the number of sequences decreased to 923, ranging from 6,390 to 48,550,500 bp (mean = 1,364,162, mode = 29,197) with a total length of 1,259,121,441 bp and an N50 of 24,795,794 bp. This corresponds to a 23.6x genome coverage if we assume that the total assembly size is similar to the true genome size, which is expected to be about 1.2 Gb based on knowledge on genomes of other related species (e.g., Alexandre et al. 2023 ). NextRAD sequencing, variant call, quality filtering and microhaplotype data After trimming, read alignment to the draft genome yielded an average of 3.5 million (M) reads per individual (minimum = 0.1 M; maximum = 5 M). FreeBayes identified 1.6 M raw variants (including SNPs, multi-nucleotide polymorphisms, indels and other complex variants) across the 190 individuals. After filtering, a total of 2,804 loci containing 6,725 SNPs across 187 individuals with a maximum of 30% missing data. Average depth per locus and per individual was 46.1 (ranging from 20.4 to 98.5 and 13.4 to 79.8, respectively). Missing data estimates across all loci was on average 0.14 for TX 2017 while for other collections ranged from 0.02 to 0.06 (Figure S1 , Supplementary material). From this dataset two others were generated. The NM dataset included all 2,804 loci and 179 individuals after TX 2017 was excluded. The NM-TX dataset was restricted to 1,908 loci with less than 20% missing data and all 187 individuals. Dataset details are presented in Table 1 . Table 1 – Number of loci and samples from each temporal collection retained in each analyzed dataset, after bioinformatics filtering. Year 2015 2016 2017 2018 2019 2020 Dataset NM (2,804 loci) Upstream 34 31 21 9 41 43 NM-TX (1,908 loci) Upstream 34 31 21 9 41 43 Downstream - - 8 - - - Patterns of temporal genomic variation Regardless of the dataset, a considerable number of samples overlapped across the DAPC space (Fig. 2 ), however the DAPC performed with the NM-TX dataset (1,908 loci) shows less variability among NM samples (Figs. 1 a and 1 b) than the DAPC performed with the NM dataset (2,804 loci; Figs. 1 c and 1 d). Using all 2,804 loci, NM samples largely overlapped from 2015 until 2018, but most of the samples showed a shift in genetic variation in 2019 and a further shift in 2020. The estimated pairwise F ST values between all possible comparisons and using fewer loci (NM-TX dataset) were very low (≤0.003) and mostly not statistically different from zero (Fig. 3 a). The exceptions were three comparisons involving mostly NM 2016 and NM 2020. When using the NM dataset, the magnitude of the F ST values was similar (≤0.002) but eight comparisons were significant, including the consecutive years 2015–2016, 2016–2017, and 2019–2020. Five comparisons involved NM 2020, and four included NM 2016 (Fig. 3 b). Regardless of the dataset used, NM 2016 and NM 2020 had more significant F ST differences when compared with other temporal/spatial samples. However, it is worth mentioning that the significance of F ST is impacted by low sample size, such as in TX 2017 and NM 2018. Mean AFD between consecutive years using both datasets showed overall small changes (ranging from 0.04 to 0.07 in both cases), but some loci exhibited substantial fluctuations (Fig. 4 ). Overall, more loci and bigger changes occurred between 2017–2018 and 2018–2019 (Fig. 4 and Table 2 ), but those values might be biased due to small sample size in 2018. However, excluding those comparisons, frequency changes of 0.1 or higher were noted in 79 to 173 loci with NM-TX dataset and 117 to 267 loci with NM dataset (Table 2 ). Moreover, 76 to 86 loci were outliers in terms of AFD between consecutive years when using the NM-TX dataset and 105 to 130 when using the NM dataset (Fig. 4 ). Most of the outlier loci were shared between datasets (Table 2 ). Moreover, regardless of the dataset, only two loci located in different contigs showed AFD of 0.1 or higher across all consecutive years. However, allele frequency changes do not show any particular trend (Figure S2, Supplementary material). Table 2 – Number of loci with absolute allele frequency differences (AFD) between consecutive years distributed in five classes of 0.1. The highest absolute allele frequency difference between two consecutive years was 0.41. Dataset AFD Class 2015–2016 2016–2017 2017–2018 2018–2019 2019–2020 NM (2,804 loci) [0,0.1[ 2605 (92.9%) 2537 (90.48%) 2236 (79.74%) 2296 (81.9%) 2687 (95.83%) [0.1,0.2[ 192 (6.85%) 251 (8.95%) 489 (17.44%) 452 (16.1%) 116 (4.13%) [0.2,0.3[ 6 (0.21%) 16 (0.57%) 66 (2.35%) 48 (1.7%) 1 (0.04%) [0.3,0.4[ 1 (0.04%) 0 12 (0.43%) 8 (0.3%) 0 [0.4,0.5[ 0 0 1 (0.04%) 0 0 Outliers 105 (3.74%) 105 (3.74%) 130 (4.64%) 126 (4.49%) 124 (4.42%) NM-TX (1,908 loci) [0,0.1[ 1763 (92.4%) 1735 (91%) 1536 (80.5%) 1578 (82.7%) 1829 (95.9%) [0.1,0.2[ 141 (7.4%) 165 (8.6%) 324 (16.98%) 285 (14.94%) 79 (4.1%) [0.2,0.3[ 4 (0.2%) 8 (0.4%) 38 (2%) 38 (1.99%) 0 [0.3,0.4[ 0 0 10 (0.52%) 7 (0.37%) 0 [0.4,0.5[ 0 0 0 0 0 Outliers 76 (3.98%) 64 (3.35%) 92 (4.82%) 88 (4.61%) 79 (4.14%) Outliers common to both datasets 76 63 86 83 79 Temporal genetic diversity, inbreeding, and relatedness Genetic diversity metrics (A R , H O , and H S ) estimated with both datasets were very similar across time and generally not significantly different between temporal collections (Fig. 5 a–c and 4 e–g). However, Kruskal-Wallis tests were significant for A R estimated with both dataset (NM-TX dataset p-value = 7.87x10 − 18 , NM dataset p-value = 1.2x10 − 11 ; Fig. 5 a and 5 d) and for H O estimated with the NM dataset (p-value = 0.02; Fig. 5 e). Pairwise Wilcoxon tests revealed that A R values were significantly different when comparing TX 2017 or NM 2018 to other populations (Tables S1 and S2, Supplementary material) and revealed H O as significantly different between NM 2018 and 2020 (Tables S3, Supplementary material). In both cases, genetic diversity metrics were significant when one of the populations had smaller sample size suggesting that sample size was a possible source of bias. Although the average inbreeding coefficients estimated with both datasets were close to zero (Fig. 5 d and 5 h), most pairwise comparisons were significantly different (Tables S4 and S5, Supplementary material). In most cases differences between collections were subtle (Fig. 5 d and 5 h), but an obvious lower genome-wide F IS was found in TX 2017 (Fig. 5 d), which had the lowest Wilcoxon test p-values with all other temporal collections (except when compared to NM 2018; Tables S4, Supplementary material). These results can be indicative of an excess of heterozygotes in many loci. Relatedness analyses were performed with a single SNP per locus from the NM-TX and NM datasets. Removing SNPs with missing data higher than 5% resulted in 1,647 and 2,136 SNPs for the NM-TX and NM datasets, respectively. Results with both datasets were very similar (Fig. 6 ). In general, sampled individuals were unrelated to each other (theoretical expectations for totally unrelated individuals are k0 ≈ 1 and k1 ≈ 0). A few pairs exhibit a probability of shared alleles similar to the expectations for first cousins (k0 = 0.75, k1 = 0.25) but values of k0 > 0.75 and k1 < 0.25 may also include more distant relatives. Two individuals exhibited probabilities of shared alleles between what is expected for first cousins and second-degree relationships (k0 = 0.5, k1 = 0.5; half-siblings, avuncular, grandchild-grandparent). Among those few pairs of potential relatives (from second degree to distant relatives), highest relatedness was between an individual from NM 2017 and an individual from NM 2020 (k0 = 0.57, k1 = 0.43 with NM-TX dataset; k0 = 0.62, k1 = 0.34 with NM dataset;) but all the other pairs with k0 ≤ 0.8 were from NM 2015 and or NM 2016. In addition, two other pairs of individuals from NM 2019 exhibited k0 ≈ 0 and k1 ≈ 0.12, meaning that the probability of sharing zero alleles by descent was close to 0% and the probability of sharing one allele was about 12%. Because this analysis used biallelic data, the probability of sharing two alleles by descent was around 88%, which means that those two pairs were either identical twins or, most likely, resulted from matings between related individuals. Temporal trends in effective population size Estimates of LD N e and 95% CIs in NM were consistent between both datasets (Fig. 7 a and Table S6, Supplementary material). This estimator suggests that the population in New Mexico increased from 2014 (LD N e = 2,255, estimated from 2015 sample) to 2015 (LD N e = 3,666, estimated from 2016 sample). Whether the trend was maintained until 2017 is not clear because point estimates and CIs of LD N e for 2016 and 2017 parental generations were infinite. Infinite values mean that LD N e values were very high or alternatively the power to estimate LD N e was low due to small sample sizes. In 2018, LD N e declined drastically (estimated from the 2019 sample) followed by a notable increase in 2019 (estimated from the 2020 sample). Estimates of N eV in NM obtained with both datasets were very small (< 351; Fig. 7 b and Table S7, Supplementary material). Only the point estimate obtained with NM-TX dataset between 2017 and 2018 was infinite (this is also likely to because of the small 2018 sample size). Estimates obtained with NM-TX dataset were similar across time-series while the NM dataset suggested a small increase in 2017–2018. Nevertheless, the confidence intervals obtained from both datasets largely overlapped, suggesting that N eV never exceeded 300, except between 2017 and 2018 (NM dataset N eV = 351; CI = [196, 1,374]). However, the sample size for those two years were the smallest and these values should be interpreted with caution. Discussion In this study, temporal trends of genomic variation, genetic diversity metrics and effective population size in M. tetranema were evaluated from samples collected from 2015 to 2020. Results from genomic variation and N e across time suggests that M. tetranema upstream (NM) is affected by genetic drift, but most genetic diversity metrics did not change significantly across the time-series. Moreover, relatedness metrics suggest that the large majority of the breeding populations upstream are unrelated. Evidence of upstream genetic drift from temporal genomic variation and effective population size The DAPC results using different datasets exhibited some differences. The NM dataset (2,804 loci) samples from 2019 and 2020 showed a higher dispersion and were shifted compared to all other samples while the NM-TX dataset (1908 loci) detected more subtle shifts. Notably, the dataset with more loci (NM dataset) revealed some changes in genomic variation not detected with the reduced dataset (NM-TX dataset). Consistent with these results, most pairwise F ST estimates were significantly different from zero when analyzing the NM dataset but not when using the NM-TX dataset. However, outlier AFD loci across time (i.e., loci with higher differences) were essentially the same in both datasets. This result suggests that the differences observed in the DAPC and F ST between datasets is likely a result of many loci with moderate AFD values present in the NM dataset, particularly after 2018, rather than few outlier loci (i.e., higher frequency changes). Moreover, regardless of the magnitude of AFD no particular trend was identified, i.e., changes were random. The moderate to large values of AFD detected in our datasets are expected when the population experiences genetic drift (Gompert et al. 2021 ). Across time-series, N eV estimates were similar and consistently low, regardless of the dataset used, suggesting substantial changes in allele frequencies between consecutive years. These results were very similar to those obtained with microsatellite data (Osborne et al. 2021 ). Low N eV estimates are a sign that genetic drift is currently a strong evolutionary force (Nei and Tajima 1981 ; Jorde and Ryman 1995 ; Gompert et al. 2021 ; Osborne et al. 2023b ). Unlike N eV , estimates of LD N e in New Mexico (upstream) showed a small increase in N e from 2014 to 2015. These results are also consistent with microsatellite data that showed a progressive increase in LD N e between 2014 and 2017 in NM (Osborne et al. 2021 ). However, the values reported here are higher and with finite confidence intervals (e.g., LD N e [NM dataset 2014] = 2,255.3, CI = [1,917.5; 2,736.5]; LD N e [microsatellites 2014] = 1,386, CI = [265; ∞]) due to the larger number of loci. Improved precision is an advantage offered by genomic time-series data when compared to microsatellite temporal data (Osborne et al. 2023b ). However, it is worth noting that when analyzing a high number of loci and relatively small number of individuals, estimates of N e based on linkage disequilibrium may be inflated (Waples and Do 2010 ; Ragsdale and Gravel 2020 ). Following the gradual increase, LD N e suffered an accentuated decline in 2018 (estimated from 2019 sample) and in the following year it was the highest recorded across the entire time-series. While the general increasing trend of LD N e described here is consistent with relative abundance data (catch-per-unit-effort, CPUE) obtained between 2012 and 2019 for the same population (Osborne et al. 2021 ), the LD N e decrease in 2018 was contrary to the increased CPUE of M. tetranema over that period. This suggests that despite increasing abundance, in some years, reproductive success may be limited to a small fraction of the population. For example, this can happen if the population is not entirely panmictic, contrary to what has generally been assumed (e.g., Osborne et al. 2021 ). Also, the LD N e decrease in 2018 is in line with the biggest changes detected with DAPC after 2018. Moreover, LD N e fluctuations are expected in pelagophilic fish species inhabiting desert rivers due to the highly variable environmental conditions (Osborne et al. 2012 , 2021 , 2023b ). Overall, the results presented here demonstrate that genetic drift is an overarching evolutionary force driving trends in genetic variation in M. tetranema . Its predominance is explained by ecological and demographic factors. For example, a regional drought between 2011 and 2013 resulted in extirpation of the population from the Ninnescah and Arkansas rivers in Kansas (Perkin et al. 2015 ; Pennock et al. 2017 ) and led to decreases in M. tetranema abundance and distribution in the South Canadian River with lowest abundance occurring in 2012 (Pennock et al. 2017 ; Osborne et al. 2021 ). Low abundance likely increased the effect of genetic drift in the population. Yet, population recovery after 2012 did not reduce the influence of genetic drift as suggested by the results presented here. In fact, some analysis suggest that genetic drift might have increased in more recent years despite increased abundance. Genetic drift and differential reproductive success (i.e., non-panmictic population) in the upstream portion of the species’ range maybe driven by downstream biased gene flow due to egg and larval drift, but this prediction lacks explicit testing in the current study. Maintenance of genetic diversity in face of genetic drift Despite signs of genetic drift, significant changes were not observed in most genetic diversity metrics between temporal collections. Exceptions were A R and H O in a single instance when sample size of one of the compared collections was small. Estimates of H S and generally H O were more robust to bias introduced by small sample sizes while A R was more sensitive as expected (Leberg 2002 ). Regarding F IS , many distributions from both datasets are statistically different, but in general follow similar patterns with F IS close to zero for most of the loci suggesting that most M. tetranema genome does not experience excesses of heterozygosity or homozygosity across time. Despite some significant tests likely due to small sample sizes, the results presented here suggest that genetic diversity in NM across time remained similar. Genetic drift and maintenance of genetic diversity in upstream might be related with movement ecology in this species. Previous field observations reported high numbers of M. tetranema larvae and young-of-the-year downstream while upstream population is dominated by adults (S. Davenport and J. Hatt, personal communication). While there’s a lack of robust data, such observation might be related with downstream drift of fertilized eggs while developing (Bottrell et al. 1964 ) if eggs and larvae retention in upstream nursery habitats is limited. Upstream spawning migration by larger, sexually mature M. tetranema prior and during reproductive season was previously suggested (Bonner 2000 ), and a recent mark-recapture study found that adults exhibit upstream-biased movement in the South Canadian River (Steffensmeier et al. 2024 ). Together, genetic and ecological data suggest that M. tetranema upstream experiences genetic drift likely due of downstream-biased gene flow, while upstream movement of adults may explain low relatedness and maintenance of genetic diversity upstream in the face of genetic drift. Conservation implications Although the results presented in this study do not suggest immediate detrimental effects on genomic diversity, abrupt demographic changes, particularly in human-mediated environmental changes is a real threat in fish with similar life history traits to M. tetranema (e.g., Osborne et al. 2023a ). Moreover, the maintenance of genetic diversity and N e upstream might be linked to upstream movement of adults. In fragmented habitats, another pelagophilic species experienced genetic diversity and N e declines upstream with persistence sustained by population augmentation (Osborne et al. In Review). Maintaining connectivity in the South Canadian River where M. tetranema persists should be a priority. Future research Temporal sampling upstream and downstream with increased sample size is needed to test the hypotheses related with asymmetric gene flow and distribution of genomic diversity across the riverscape. This would allow more robust comparisons of genetic variability and the effect of genetic drift across time and space and would also allow estimation of the magnitude and direction of movement between up and downstream sites. Moreover, the transition from microsatellite-based genetic monitoring to SNP-based assessments will be an important component of future conservation and management efforts. This requires the development of a panel of markers that can be consistently monitored across time. A SNP-based panel can be used for both long-term genetic monitoring and to obtain genetic data to test various hypotheses. Genotyping-in-Thousands by sequencing (GT-seq, Campbell et al. 2015 ) is an efficient and cost-effective method of targeted SNP genotyping that uses multiplexed PCR amplicon sequencing and facilitates simultaneous amplification of hundreds of targeted genetic loci in hundreds to thousands of samples. The SNP and microhaplotype data developed in this study offer an excellent opportunity to identify and select a set of loci that can be used to achieve those research goals. Declarations Acknowledgments We sincerely thank Joanna Hatt, Andrew Monie, John Caldwell, Eliza Gilbert (New Mexico Department of Game and Fish), David Kitcheyan, Stephen Davenport, and Daniel Fenner (U.S. Fish and Wildlife Service) for samples collection. We thank Emily DeArmon and Alexandra Snyder (UNM Museum of Southwestern Biology Division of Fishes) for expert curatorial services. High performance computing was conducted at the UNM Center for Advanced Research Computing that is supported, in part, by the National Science Foundation. This study was funded by the New Mexico Department of Game and Fish. Founded in 1889, the University of New Mexico sits on the traditional homelands of the Pueblo of Sandia. The original peoples of New Mexico – Pueblo, Navajo, and Apache – since time immemorial, have deep connections to the land and have made significant contributions to the broader community statewide. We honor the land itself and those who remain stewards of this land throughout the generations and also acknowledge our committed relationship to Indigenous peoples. We gratefully recognize our history. Funding This study was funded by the New Mexico Department of Game and Fish. Competing Interests The authors declare no conflicts of interest. Author Contributions G.C.D. was responsible for project design, bioinformatics, data analysis, and writing the manuscript. A.C. was responsible for bioinformatics and critical revision of the manuscript, T.T. was responsible for critical revision of the manuscript. M.O. was responsible for project conception, design and coordination, laboratory work, critical revision of the manuscript, and obtaining funding. Data Availability Individual genotype data from microhaplotypes are available from authors upon request. Raw sequence reads from nextRAD-seq (190 individuals), are deposited in the NCBI Sequence Read Archive (SRA), with the BioProject accession number XXXXXX (to be added upon acceptance). References Alexandre NM, Cameron AC, Tian D, et al (2023) Chromosome-level reference genomes of two imperiled desert fishes: spikedace ( Meda fulgida ) and loach minnow ( Tiaroga cobitis ). G3: Genes, Genomes, Genetics 13:jkad157 Alonge M, Lebeigle L, Kirsche M, et al (2022) Automated assembly scaffolding using RagTag elevates a new tomato system for high-throughput genome editing. Genome biology 23:258 Alp M, Keller I, Westram AM, Robinson CT (2012) How river structure and biological traits influence gene flow: a population genetic study of two stream invertebrates with differing dispersal abilities. Freshwater Biology 57:969–981 Archdeacon TP, Davenport SR, Grant JD, Henry EB (2018) Mass upstream dispersal of pelagic-broadcast spawning cyprinids in the Rio Grande and Pecos River, New Mexico. Western North American Naturalist 78:100–105 Archdeacon TP, Diver-Franssen TA, Bertrand NG, Grant JD (2020) Drought results in recruitment failure of Rio Grande silvery minnow ( Hybognathus amarus ), an imperiled, pelagic broadcast-spawning minnow. Environmental Biology of Fishes 103:1033–1044 Balon EK (1975) Reproductive guilds of fishes: a proposal and definition. Journal of the Fisheries Board of Canada 32:821–864 Balon EK (1981) Additions and amendments to the classification of reproductive styles in fishes. Environmental Biology of Fishes 6:377–389 Bergland AO, Behrman EL, O’Brien KR, et al (2014) Genomic evidence of rapid and stable adaptive oscillations over seasonal time scales in Drosophila . PLoS genetics 10:e1004775 Bilton TP, McEwan JC, Clarke SM, et al (2018) Linkage disequilibrium estimation in low coverage high-throughput sequencing data. Genetics 209:389–400 Blanchet S, Prunier JG, Paz‐Vinas I, et al (2020) A river runs through it: The causes, consequences, and management of intraspecific diversity in river networks. Evolutionary Applications 13:1195–1213 Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114–2120 Bonner TH (2000) Life history and reproductive ecology of the Arkansas River shiner and peppered chub in the Canadian River, Texas and New Mexico. Texas Tech University Bonner TH, Wilde GR (2000) Changes in the fish assemblage of the Canadian River, Texas, associated with reservoir construction. Journal of Freshwater Ecology 15:189–198 Bottrell CE, Ingersol RH, Jones RW (1964) Notes on the embryology, early development, and behavior of Hybopsis aestivalis tetranemus (Gilbert). Transactions of the American Microscopical Society 83:391–399 Breiman L (2001) Random forests. Machine Learning 45:5–32 Caeiro-Dias G, Osborne M, Turner T (2024) Time is of the essence: using archived samples in the development a GT-seq panel to preserve continuity of ongoing genetic monitoring. Authorea . Campbell NR, Harmon SA, Narum SR (2015) Genotyping‐in‐Thousands by sequencing (GT‐seq): A cost effective SNP genotyping method based on custom amplicon sequencing. Molecular Ecology Resources 15:855–867 Chavez MJ, Budy P, Pennock CA, et al (2024) Movement patterns of a small-bodied minnow suggest nomadism in a fragmented, desert river. Movement Ecology 12:52 Cheng H, Concepcion GT, Feng X, et al (2021) Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nature Methods 18:170–175 Danecek P, Auton A, Abecasis G, et al (2011) The variant call format and VCFtools. Bioinformatics 27:2156–2158 Danecek P, Bonfield JK, Liddle J, et al (2021) Twelve years of SAMtools and BCFtools. GigaScience 10:giab008 Demandt MH (2010) Temporal changes in genetic diversity of isolated populations of perch and roach. Conservation Genetics 11:249–255 DiBattista JD (2008) Patterns of genetic variation in anthropogenically impacted populations. Conservation Genetics 9:141–156 Do C, Waples RS, Peel D, et al (2014) NeEstimator v2: re‐implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources 14:209–214 Dudley RK, Platania SP (2007) Flow regulation and fragmentation imperil pelagic‐spawning riverine fishes. Ecological Applications 17:2074–2086 Eisenhour DJ (1999) Systematics of Macrhybopsis tetranema (Cypriniformes: Cyprinidae). Copeia 969–980 Franklin PA, Bašić T, Davison PI, et al (2024) Aquatic connectivity: challenges and solutions in a changing climate. Journal of Fish Biology 105:392–411 Galván‐Femenía I, Graffelman J, Barceló‐i‐Vidal C (2017) Graphics for relatedness research. Molecular Ecology Resources 17:1271–1282 Garrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. arXiv preprint arXiv:12073907 Gilbert KJ, Whitlock MC (2015) Evaluating methods for estimating local effective population size with and without migration. Evolution 69:2154–2166 Gompert Z, Springer A, Brady M, et al (2021) Genomic time‐series data show that gene flow maintains high genetic diversity despite substantial genetic drift in a butterfly species. Molecular Ecology 30:4991–5008 Hill WG (1981) Estimation of effective population size from data on linkage disequilibrium1. Genetics Research 38:209–216 Jombart T (2008) adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403–1405 Jombart T, Ahmed I (2011) adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27:3070–3071 Jones AT, Ovenden JR, Wang Y-G (2016) Improved confidence intervals for the linkage disequilibrium method for estimating effective population size. Heredity 117:217–223 Jorde PE, Ryman N (1995) Temporal allele frequency change and estimation of effective size in populations with overlapping generations. Genetics 139:1077–1090 Keenan K, McGinnity P, Cross TF, et al (2013) diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods in ecology and evolution 4:782–788 Kuczynski L, Legendre P, Grenouillet G (2018) Concomitant impacts of climate change, fragmentation and non‐native species have led to reorganization of fish communities since the 1980s. Global Ecology and Biogeography 27:213–222 Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nature Methods 9:357 Leberg P (2002) Estimating allelic richness: effects of sample size and bottlenecks. Molecular Ecology 11:2445–2449 Li H (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34:3094–3100 Li H, Handsaker B, Wysoker A, et al (2009) The sequence alignment/map format and SAMtools. Bioinformatics 25:2078–2079 Liaw A, Wiener M (2002) Classification and regression by randomForest. R news 2:18–22 Luttrell GR, Echelle AA, Fisher WL, Eisenhour DJ (1999) Declining status of two species of the Macrhybopsis aestivalis complex (Teleostei: Cyprinidae) in the Arkansas River basin and related effects of reservoirs as barriers to dispersal. Copeia 981–989 Mangiafico S (2021) package rcompanion: Functions to Support Extension Education Program Evaluation. Rutgers Cooperative Extension, New Brunswick, New Jersey Version 2:30 Marçais G, Kingsford C (2011) A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics 27:764–770 Marschall EA, Crowder LB (1996) Assessing population responses to multiple anthropogenic effects: a case study with brook trout. Ecological Applications 6:152–167 Mathieu‐Bégné E, Loot G, Chevalier M, et al (2019) Demographic and genetic collapses in spatially structured populations: Insights from a long‐term survey in wild fish metapopulations. Oikos 128:196–207 Meirmans PG (2020) Genodive version 3.0: Easy‐to‐use software for the analysis of genetic data of diploids and polyploids. Molecular Ecology Resources 20:1126–1131 Messer PW, Ellner SP, Hairston NG (2016) Can population genetics adapt to rapid evolution? Trends in Genetics 32:408–418 Morrissey MB, de Kerckhove DT (2009) The maintenance of genetic variation due to asymmetric gene flow in dendritic metapopulations. The American Naturalist 174:875–889 Nei M (1987) Molecular evolutionary genetics. Columbia University Press Nei M, Tajima F (1981) Genetic drift and estimation of effective population size. Genetics 98:625–640 Nunn AD, Moccetti P, Hänfling B, et al (2024) The genome sequence of the Eurasian minnow, Phoxinus phoxinus (Linnaeus, 1758). Wellcome Open Research 9:504 Osborne MJ, Archdeacon TP, Yackulic CB, et al (2023a) Genetic erosion in an endangered desert fish during a megadrought despite long‐term supportive breeding. Conservation Biology e14154 Osborne MJ, Caeiro‐Dias G, Turner TF (2023b) Transitioning from microsatellites to SNP‐based microhaplotypes in genetic monitoring programmes: Lessons from paired data spanning 20 years. Molecular Ecology 32:316–334 Osborne MJ, Carson EW, Turner TF (2012) Genetic monitoring and complex population dynamics: insights from a 12‐year study of the Rio Grande silvery minnow. Evolutionary Applications 5:553–574 Osborne MJ, Hatt JL, Gilbert EI, Davenport SR (2021) Still time for action: genetic conservation of imperiled South Canadian River fishes, Arkansas River Shiner ( Notropis girardi ), Peppered Chub ( Macrhybopsis tetranema ) and Plains Minnow (Hybognathus placitus). Conservation Genetics 22:927–945 Osborne MJ, Hudgell MAB, Caeiro-Dias G, Turner TF (2023c) The complete mitochondrial genomes of two imperiled species endemic to the Southwestern United States: Peppered Chub ( Macrhybopsis tetranema ) and Gila Trout ( Oncorhynchus gilae ). Mitochondrial DNA Part B 8:809–814 Paradis E (2010) pegas: an R package for population genetics with an integrated–modular approach. Bioinformatics 26:419–420 Paz-Vinas I, Loot G, Hermoso V, et al (2018) Systematic conservation planning for intraspecific genetic diversity. Proceedings of the Royal Society B: Biological Sciences 285:20172746 Pennock CA, Gido KB, Perkin JS, et al (2017) Collapsing range of an endemic Great Plains minnow, Peppered Chub Macrhybopsis tetranema . The American Midland Naturalist 177:57–68 Perkin JS, Gido KB (2011) Stream fragmentation thresholds for a reproductive guild of Great Plains fishes. Fisheries 36:371–383 Perkin JS, Gido KB, Costigan KH, et al (2015) Fragmentation and drying ratchet down Great Plains stream fish diversity. Aquatic Conservation: Marine and Freshwater Ecosystems 25:639–655 Perkin JS, Starks TA, Pennock CA, et al (2019) Extreme drought causes fish recruitment failure in a fragmented Great Plains riverscape. Ecohydrology 12:e2120 Platania SP, Altenbach CS (1998) Reproductive strategies and egg types of seven Rio Grande basin cyprinids. Copeia 559–569 Platania SP, Mortensen JG, Farrington MA, et al (2020) Dispersal of stocked Rio Grande silvery minnow ( Hybognathus amarus ) in the middle Rio Grande, New Mexico. The Southwestern Naturalist 64:31–42 Ragsdale AP, Gravel S (2020) Unbiased estimation of linkage disequilibrium from unphased data. Molecular Biology and Evolution 37:923–932 Ranallo-Benavidez TR, Jaron KS, Schatz MC (2020) GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nature Communications 11:1432 Russello MA, Waterhouse MD, Etter PD, Johnson EA (2015) From promise to practice: pairing non-invasive sampling with genomics in conservation. PeerJ 3:e1106 Ryan SF, Deines JM, Scriber JM, et al (2018) Climate-mediated hybrid zone movement revealed with genomics, museum collection, and simulation modeling. Proceedings of the National Academy of Sciences 115:E2284–E2291 Steffensmeier ZD, Mayes KB, Perkin JS (2024) Linking short-term movement rate of pelagic-broadcast spawning fishes to river fragment length and conservation status. Biological Conservation 293:110585 Stuart IG, Sharpe CP (2020) Riverine spawning, long distance larval drift, and floodplain recruitment of a pelagophilic fish: a case study of golden perch ( Macquaria ambigua ) in the arid Darling River, Australia. Aquatic Conservation: Marine and Freshwater Ecosystems 30:675–690 Torterotot J-B, Perrier C, Bergeron NE, Bernatchez L (2014) Influence of forest road culverts and waterfalls on the fine-scale distribution of brook trout genetic diversity in a boreal watershed. Transactions of the American Fisheries Society 143:1577–1591 U.S. Fish and Wildlife Service (2018) Species status assessment report for the Arkansas River shiner ( Notropis girardi ) and peppered chub ( Macrhybopsis tetranema ), version 1.0, with appendices. October 2018. Albuquerque, NM. 172 pp. U. S. Fish and Wildlife Service (2022) Endangered and threatened wildlife and plants; endangered species status for the peppered chub and designation of critical habitat. Federal Register 87:11188–11220. Waples RS (2005) Genetic estimates of contemporary effective population size: to what time periods do the estimates apply? Molecular Ecology 14:3335–3352 Waples RS, Do C (2010) Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: a largely untapped resource for applied conservation and evolution. Evolutionary Applications 3:244–262 Wilde GR, Durham BW (2008) A life history model for peppered chub, a broadcast-spawning cyprinid. Transactions of the American Fisheries Society 137:1657–1666 Willis SC, Hollenbeck CM, Puritz JB, et al (2017) Haplotyping RAD loci: an efficient method to filter paralogs and account for physical linkage. Molecular Ecology Resources 17:955–965 Yackulic CB, Archdeacon TP, Valdez RA, et al (2022) Quantifying flow and nonflow management impacts on an endangered fish by integrating data, research, and expert opinion. Ecosphere 13:e4240 Zheng X, Levine D, Shen J, et al (2012) A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28:3326–3328 Additional Declarations No competing interests reported. Supplementary Files MtetgenomictimeseriessuppmaterialConGen081425.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7863499","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":532972543,"identity":"718b4fb8-e84e-46e2-a774-eb36497982f8","order_by":0,"name":"Guilherme Caeiro-Dias","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA1klEQVRIiWNgGAWjYLCCBAYGxjYgfgDnEquF2YB4LUDA2MDAwCZBlBZ598MPHzzcYyPbx378WTXvjm0M/Ow5Bni1GJ5JMzZIeJZm3MaTY3ab98xtBsmeNwS0zGAwk0g4cDixjSGH7TZv220GgxuEbJnB/v0HWAv/82fFIC32hLTIS/CYMYC1SCSYMYNtkSCgxYAnpxjoMKBfJN4YS85tu80jceZZAX5b2o9v/PjjgI3s/P70hx/ett2W429P3oDflgNoAjx4lYNtaSCoZBSMglEwCkY8AABsLUhQJJNoGwAAAABJRU5ErkJggg==","orcid":"","institution":"University of New Mexico","correspondingAuthor":true,"prefix":"","firstName":"Guilherme","middleName":"","lastName":"Caeiro-Dias","suffix":""},{"id":532972544,"identity":"e700a9f2-a76c-46af-8468-496c9deb4fd4","order_by":1,"name":"Alexander Cameron","email":"","orcid":"","institution":"Arizona Game and Fish Department","correspondingAuthor":false,"prefix":"","firstName":"Alexander","middleName":"","lastName":"Cameron","suffix":""},{"id":532972545,"identity":"3430ad7f-7a44-4a1f-b51e-f1d23c9a2466","order_by":2,"name":"Thomas Turner","email":"","orcid":"","institution":"University of New Mexico","correspondingAuthor":false,"prefix":"","firstName":"Thomas","middleName":"","lastName":"Turner","suffix":""},{"id":532972546,"identity":"c65d5edb-d1e5-4a59-a717-1ae92e89182a","order_by":3,"name":"Megan Osborne","email":"","orcid":"","institution":"University of New Mexico","correspondingAuthor":false,"prefix":"","firstName":"Megan","middleName":"","lastName":"Osborne","suffix":""}],"badges":[],"createdAt":"2025-10-15 04:08:33","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7863499/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7863499/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":94095429,"identity":"6d395ff3-2744-460a-9adc-50e8ca6bf65f","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3789061,"visible":true,"origin":"","legend":"","description":"","filename":"MtetgenomictimeseriesConGen081425.docx","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/4bcbf0875c692ea78b243640.docx"},{"id":94095415,"identity":"b35c0dd7-0414-4c93-96b9-175cca41a684","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"json","order_by":1,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":7402,"visible":true,"origin":"","legend":"","description":"","filename":"93ebd77084e0466cb121314b8a366f70.json","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/ad1d4929b6ae31e464074bbc.json"},{"id":94096237,"identity":"74c98d80-5e22-4144-92be-5e2d0f58791e","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":110503,"visible":true,"origin":"","legend":"","description":"","filename":"MtetgenomictimeseriessuppmaterialConGen081425.docx","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/712ee49861f167fc97ac0228.docx"},{"id":94096526,"identity":"aed29a00-cbda-42b8-848b-fa551d3fdd79","added_by":"auto","created_at":"2025-10-22 10:07:05","extension":"xml","order_by":3,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":167056,"visible":true,"origin":"","legend":"","description":"","filename":"93ebd77084e0466cb121314b8a366f701enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/46dce904436cc3c842190d79.xml"},{"id":94096529,"identity":"ad82fdab-42e2-4613-86c3-07153e9b6ea6","added_by":"auto","created_at":"2025-10-22 10:07:05","extension":"emf","order_by":4,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3150864,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage1.emf","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/d10b5b84c7f3312f3a71d904.emf"},{"id":94095424,"identity":"eed09a63-59c2-4b1e-9812-f3fd790c2b8d","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":5,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":451299,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/3b883c2f1a66ed3295c027fb.png"},{"id":94097410,"identity":"c457b16d-5e24-4644-932f-e023ea4e5529","added_by":"auto","created_at":"2025-10-22 10:15:05","extension":"png","order_by":6,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":94150,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/3f05f6a79d985191c0e0c981.png"},{"id":94096239,"identity":"65c761e4-747c-499e-9616-a9bfdb139b7e","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"png","order_by":7,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":90291,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/94343650494091a4822c1b81.png"},{"id":94097411,"identity":"ff548071-b545-4d81-8321-35823cc61bc6","added_by":"auto","created_at":"2025-10-22 10:15:05","extension":"png","order_by":8,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":467491,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/bc732c001d50589d6adea353.png"},{"id":94096241,"identity":"aaa76637-54d4-46c6-895a-24775f285bae","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"png","order_by":9,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":118945,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/e95b070ea26f64a6d77cf5c8.png"},{"id":94095431,"identity":"13abd56d-1baf-4fc4-8212-eb83686074ea","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":10,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":100291,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/c06d73d6b4e27583d420c1e2.png"},{"id":94096245,"identity":"3f91c124-b2f0-4351-8712-6021d8413bf7","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"png","order_by":11,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":14945,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/52303c41904b2e5ecf39d2d2.png"},{"id":94096531,"identity":"ef294b80-b02c-4ff1-acab-60c64e2b9645","added_by":"auto","created_at":"2025-10-22 10:07:05","extension":"png","order_by":12,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":86854,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/0e5b5026681d1055002df9ad.png"},{"id":94096527,"identity":"564e62b4-3d9d-4595-bf36-d27ec3ea38a3","added_by":"auto","created_at":"2025-10-22 10:07:05","extension":"png","order_by":13,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":21583,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/4b1885ba1df3d2160c47b885.png"},{"id":94095428,"identity":"235b04c6-f5fc-42c7-a0f3-d81ba641c336","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":14,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":43184,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/0212d74d10b0b43956af2e62.png"},{"id":94095435,"identity":"1bf8f9c6-c0e3-4ee9-8b38-81cbc5d98cfa","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":15,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":248053,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/8ff85304dfae9c64d46587a1.png"},{"id":94095434,"identity":"933ee521-8b9f-4b85-b6c4-8a4adca1251e","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":16,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":31282,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/db5ae868e651cfa4caf1305e.png"},{"id":94095438,"identity":"b0275bfe-5ef2-4795-b02f-d23d08cc6b5a","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":17,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":41991,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/152f3a0959452613c9bee9b0.png"},{"id":94096247,"identity":"2698669e-63bf-447a-b8aa-df8af7a89bf7","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"xml","order_by":18,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":164910,"visible":true,"origin":"","legend":"","description":"","filename":"93ebd77084e0466cb121314b8a366f701structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/a22394d75ee783f785a01f63.xml"},{"id":94095439,"identity":"404458aa-aa78-43ef-bace-3ff9fd05b82a","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"html","order_by":19,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":173329,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/6186a75ad427b056ea86e7ee.html"},{"id":94095413,"identity":"3260ebc2-1636-49dc-b529-415696d2b532","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":58251,"visible":true,"origin":"","legend":"\u003cp\u003eMap showing remnant geographic distribution of\u003cstrong\u003e \u003c/strong\u003e\u003cem\u003eMacrhybopsis tetranema \u003c/em\u003ein South Canadian River between Ute Lake (New Mexico) and Lake Meredith (Texas) and the geographic origin of samples used for nextRAD sequencing. This region is highlighted inset map.\u003c/p\u003e","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/573d16d5b53037f4ad25d927.png"},{"id":94096236,"identity":"6c07369a-f823-4485-8cc7-b20e9661716a","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":360364,"visible":true,"origin":"","legend":"\u003cp\u003eResults of discriminant analysis of principal components (DAPC) based on both datasets used. \u003cstrong\u003ea)\u003c/strong\u003e Distribution of each individual included in the NM-TX dataset (1,908 loci) across the DAPC space; samples from New Mexico (NM) are represented with circles and those from Texas (TX) are represented with triangles. \u003cstrong\u003eb)\u003c/strong\u003e Ellipses representing the 95% confidence level (CL) for a multivariate normal distribution (MND) for the NM-TX dataset. \u003cstrong\u003ec)\u003c/strong\u003e Distribution of each individual included in the NM dataset (2,804 loci) across the DAPC space. \u003cstrong\u003ed)\u003c/strong\u003e Ellipses representing the 95% CL for a MND for the NM dataset. The percentage values in the axis of each plot refers to the variance explained by each discriminant function.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/d72439bc2b33534c79575171.png"},{"id":94096235,"identity":"30ab0da3-515d-4125-8967-dc4ae22afec6","added_by":"auto","created_at":"2025-10-22 09:59:05","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":131666,"visible":true,"origin":"","legend":"\u003cp\u003ePairwise F\u003csub\u003eST\u003c/sub\u003e comparisons \u003cstrong\u003ea)\u003c/strong\u003e obtained from the NM-TX dataset (1,908 loci) and \u003cstrong\u003eb)\u003c/strong\u003e obtained from NM dataset (2,804 loci). Significant F\u003csub\u003eST\u003c/sub\u003e values are highlighted with an asterisk (*).\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/7139bb20387512f0bf41a4cb.png"},{"id":94095416,"identity":"3644e49f-2144-4951-9700-a5eccbed3e4e","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":92526,"visible":true,"origin":"","legend":"\u003cp\u003eViolin plots and box plots showing absolute allele frequency difference (AFD) in upstream population between consecutive years using \u003cstrong\u003ea)\u003c/strong\u003e the NM-TX dataset (1,908 loci) and \u003cstrong\u003eb)\u003c/strong\u003ethe NM dataset (2,804 loci). Gray area represents the density of loci with similar ADF; dots are the boxplot ADF outlier loci.\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/8c47796d9df7566c946f4d56.png"},{"id":94095421,"identity":"5d572072-f567-4971-ab99-4ad46c5c9e16","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":452884,"visible":true,"origin":"","legend":"\u003cp\u003eMetrics of genetic diversity (allelic richness [A\u003csub\u003eR\u003c/sub\u003e], observed heterozygosity [H\u003csub\u003eO\u003c/sub\u003e], expected heterozygosity [H\u003csub\u003eS\u003c/sub\u003e], and inbreeding coefficient [F\u003csub\u003eIS\u003c/sub\u003e]) estimated across time-series with \u003cstrong\u003ea)\u003c/strong\u003e, \u003cstrong\u003eb)\u003c/strong\u003e, \u003cstrong\u003ec)\u003c/strong\u003e,\u003cstrong\u003e d)\u003c/strong\u003e the NM-TX dataset (1,908 loci) and with \u003cstrong\u003ee)\u003c/strong\u003e, \u003cstrong\u003ef)\u003c/strong\u003e, \u003cstrong\u003eg)\u003c/strong\u003e,\u003cstrong\u003e h)\u003c/strong\u003e the NM dataset (2,804 loci). Kruskal-Wallis tests and associated p-values for each metric estimated with each dataset to tested if the distribution of values were similar across time are shown on each plot; results of Wilcoxon tests for multiple comparisons for significant Kruskal-Wallis tests are provided in Tables S1 – S5, Supplementary material).\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/3ea541f4197c5f4c0f15aaa3.png"},{"id":94095419,"identity":"71610129-0c21-46f2-96bc-8021aedfe94e","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":165931,"visible":true,"origin":"","legend":"\u003cp\u003eMaximum likelihood estimations of identity-by-descent for each pair of individuals obtained from \u003cstrong\u003ea)\u003c/strong\u003eNM-TX dataset (1,647 SNPs with missingness lower than 5%) and \u003cstrong\u003eb)\u003c/strong\u003e NM dataset (2,136 with missingness lower than 5%). Each circle represents a pairwise comparison and the red line represents the maximum probability of two individuals sharing zero alleles (k0) or one allele (k1) by descent.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/e75596c7f86a1faa57304f75.png"},{"id":94095423,"identity":"a0fe3f62-166f-4dd0-b24f-93f1b46f947f","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":111219,"visible":true,"origin":"","legend":"\u003cp\u003eEffective population sizes and associated parametric 95% confidence intervals (CIs) estimated from both TX-NM (1,908 loci) and NM (2,804 loci) datasets. \u003cstrong\u003ea)\u003c/strong\u003e Linkage disequilibrium effective population size (LD N\u003csub\u003ee\u003c/sub\u003e). Note that the y-axis scale is pseudo-log transformed to distinguish smaller LD N\u003csub\u003ee\u003c/sub\u003e values and to include zero; also, the x-axis refers to the parental population that originated the population from where the samples used to estimate LD N\u003csub\u003ee\u003c/sub\u003e was collected. Point estimates and CIs values are provided in Table S4, Supplementary material. \u003cstrong\u003eb)\u003c/strong\u003e Variance effective population size (N\u003csub\u003eeV\u003c/sub\u003e) estimated using the temporal method of Nei and Tajima (1981). Point estimates and CIs values are provided in Table S5, Supplementary material. In both cases, CIs extending beyond y-axis limit and without point estimates represent infinite N\u003csub\u003ee\u003c/sub\u003e estimates.\u003c/p\u003e","description":"","filename":"floatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/35aba9f9d00ec1348326c52d.png"},{"id":96711261,"identity":"6b4a894f-b1d4-4dad-887c-90b1937a63c7","added_by":"auto","created_at":"2025-11-25 10:11:49","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2379932,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/43c2102b-8f5c-4934-ac9d-9472a636bc1b.pdf"},{"id":94095414,"identity":"412d4b8e-1490-48a4-866d-27eb193d175f","added_by":"auto","created_at":"2025-10-22 09:51:05","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":110503,"visible":true,"origin":"","legend":"","description":"","filename":"MtetgenomictimeseriessuppmaterialConGen081425.docx","url":"https://assets-eu.researchsquare.com/files/rs-7863499/v1/0e9efec3372baaccbdea5af4.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Sustaining genetic diversity in an imperiled pelagophilic fish despite genetic drift","fulltext":[{"header":"Introduction","content":"\u003cp\u003eContemporary patterns of intraspecific genetic variation are the result of demographic and evolutionary processes. Identifying specific processes generating these patterns is critically important to ecological, evolutionary and conservation studies. In particular, contemporary demographic and microevolutionary processes may impact recent patterns of genetic variation (Messer et al. \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). Time-series population-level genetic data can reveal patterns of allele frequency change over contemporary time scales and can be important for understanding recent evolutionary processes, including short term maintenance (Osborne et al. \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2012\u003c/span\u003e) and loss (Osborne et al. \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2023a\u003c/span\u003e) of genetic diversity in the face of wild population augmentation, differences in genetic declines in distinct subpopulations in a metapopulation (Mathieu-B\u0026eacute;gn\u0026eacute; et al. \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), or genetic drift after isolation (Demandt \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). Although genome-wide time-series data are slowly accumulating, they are still uncommon in the literature. Nonetheless, such studies have uncovered evidence of rapid adaptation (Bergland et al. \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2014\u003c/span\u003e), climate-mediated hybrid zone movement (Ryan et al. \u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e2018\u003c/span\u003e), maintenance of high genetic diversity due to gene flow despite genetic drift (Gompert et al. \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2021\u003c/span\u003e), and genomic signatures of genetic drift and increased inbreeding in the face of massive population declines despite population supplementation efforts (Osborne et al. \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e).\u003c/p\u003e\u003cp\u003ePatterns of genetic variation within species may also differ in space. In lotic freshwater ecosystems with dendritic connectivity (e.g., rivers, streams) several empirical studies found an increase of genetic diversity downstream (Alp et al. \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Torterotot et al. \u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e2014\u003c/span\u003e; Paz-Vinas et al. \u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). While upstream-directed colonization can replenish diversity at upstream sites, the distribution of genetic diversity across the riverscape also depends on species\u0026rsquo; intrinsic traits such as fecundity, reproductive strategy, and migratory behavior (Morrissey and de Kerckhove \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2009\u003c/span\u003e). Anthropogenic changes can further impact patterns of genetic diversity in space and time in riverine systems. Environmental alterations induced by anthropogenic activities such as habitat changes, climate change, and/or introduction of non-native species, together with stochastic events have disrupted freshwater fish population connectivity, altered population dynamics (Kuczynski et al. \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Chavez et al. \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Franklin et al. \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2024\u003c/span\u003e), and caused substantial changes in abundance (Marschall and Crowder \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e1996\u003c/span\u003e; Yackulic et al. \u003cspan citationid=\"CR81\" class=\"CitationRef\"\u003e2022\u003c/span\u003e), which in turn can affect levels of genetic diversity (DiBattista \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2008\u003c/span\u003e) and its distribution across the riverscape (Blanchet et al. \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2020\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eAcross the Great Plains of North America, fishes belonging to the pelagic broadcast spawning reproductive guild have declined in range and abundance (Perkin et al. \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Pelagophils release neutrally buoyant eggs that are fertilized once released (Platania and Altenbach \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e1998\u003c/span\u003e) and drift while developing (Bottrell et al. \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e1964\u003c/span\u003e). \u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e (peppered chub), a pelagic spawning minnow, was historically distributed throughout the Arkansas River basin in New Mexico (NM), Kansas, Texas (TX), and Oklahoma (Eisenhour \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e1999\u003c/span\u003e). This species is now extirpated from the majority of its historical range. Until about a decade ago, three isolated populations remained but a drought cycle in the late 1980\u0026rsquo;s and early 1990\u0026rsquo;s resulted in decline and eventual extirpation (last seen in 2011) from the Cimarron River (U.S. Fish and Wildlife Service \u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) and further drought between 2011 and 2013 resulted in extirpation of the Ninnescah and Arkansas river population in Kansas (Perkin et al. \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Pennock et al. \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Today a single extant population persists in a 218 km stretch of the South Canadian River between Ute Lake (NM) and Lake Meredith (TX) (Bonner and Wilde \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Pennock et al. \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). As such, \u003cem\u003eM. tetranema\u003c/em\u003e was listed as federally endangered (U.S. Fish and Wildlife Service 2022).\u003c/p\u003e\u003cp\u003e\u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e evolved life-history traits that facilitated its persistence in rivers with marked seasonality and a highly variable flow regime. Traits include small body size, short generation time (1\u0026ndash;3 years; Wilde and Durham \u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e2008\u003c/span\u003e), and non-adhesive semi-buoyant eggs characteristic of the pelagophilic reproductive guild of fishes (Balon \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e1975\u003c/span\u003e, \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e1981\u003c/span\u003e). Because eggs drift while developing, persistence of upstream populations must depend on the retention of eggs and larvae in local nursery habitats or upstream dispersal. Upstream movement was observed in post-larval life stages of several pelagophilic species (Archdeacon et al. \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Platania et al. \u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) including young-of-the-year (Stuart and Sharpe \u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). River fragmentation by human-made barriers (e.g., impoundments, dams) disrupt upstream dispersal. These barriers and the habitat changes they cause, contribute to declines of pelagophilic minnows including \u003cem\u003eM. tetranema\u003c/em\u003e (Luttrell et al. \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e1999\u003c/span\u003e; Dudley and Platania \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Perkin and Gido \u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Perkin et al. \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Archdeacon et al. \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). River fragmentation impacts habitat availability and quality, reduces stream discharge, increases frequency of stream dewatering, alters flow periodicity, and increases stream channelization (Hoagstrom and Turner 2015). Together these factors disrupt the reproductive cycle, decrease available nursery habitat, and enhance probabilities of recruitment failure in pelagic broadcast spawners (Dudley and Platania \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Archdeacon et al. \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2020\u003c/span\u003e)d \u003cem\u003etetranema\u003c/em\u003e is no exception (Wilde and Durham \u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Perkin et al. \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). For a short-lived minnow species, recruitment failure contributes to population decline and potential reduction of genetic diversity. The combination of fragmented habitats, limited geographic distribution, short lifespan and recruitment failure increase the chance of local population extirpation (Perkin et al. \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Pennock et al. \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2017\u003c/span\u003e) and, in the worst case, can lead to the extinction of this species.\u003c/p\u003e\u003cp\u003eUsing microsatellite and mtDNA data collected between 2015 and 2018, Osborne et al. (\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) found that the remnant \u003cem\u003eM. tetranema\u003c/em\u003e population maintained genetic diversity across the time-series and that effective population size (N\u003csub\u003ee\u003c/sub\u003e) estimates showed an increasing trend consistent with the demographic trajectory. We re-evaluated putatively neutral genetic variability and effective size (N\u003csub\u003ee\u003c/sub\u003e) trends using genome-wide time-series data from 2015 to 2020. Specifically, we (1) evaluated the magnitude of genetic drift from temporal sampling, (2) tested if there was a significant loss of genetic diversity across the extended time-series, and (3) assessed temporal trends of N\u003csub\u003ee\u003c/sub\u003e. Results from genome-wide data provided a more comprehensive understanding of recent changes of neutral genetic diversity in \u003cem\u003eM. tetranema\u003c/em\u003e but also raised new questions about population dynamics in the remnant range in the South Canadian River.\u003c/p\u003e"},{"header":"Material and Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003eSampling\u003c/h2\u003e\u003cp\u003eTemporal sampling for reduced-representation sequencing included archived DNA isolates from 2015 (n\u0026thinsp;=\u0026thinsp;34), 2016 (n\u0026thinsp;=\u0026thinsp;31), 2017 (n\u0026thinsp;=\u0026thinsp;21), and 2018 (n\u0026thinsp;=\u0026thinsp;9), and new samples from 2019 (n\u0026thinsp;=\u0026thinsp;41) and 2020 (n\u0026thinsp;=\u0026thinsp;43) collected from multiple localities across the South Canadian River in New Mexico (NM) from Ute Late to NM-TX border (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). This region is referend in this study as \u0026ldquo;upstream\u0026rdquo;. Because it was not possible to obtain recent collections from Texas, analysis from downstream of the distribution range was restricted to relatively small number of archived samples collected in 2017 (n\u0026thinsp;=\u0026thinsp;10) from a single site (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFor the draft genome sequencing, a \u003cem\u003eM. tetranema\u003c/em\u003e individual was collected from the South Canadian River, NM (approximate geographic coordinates: 35.389376, \u0026minus;\u0026thinsp;103.355331) in August 2021. Species identity was determined by visual examination of external morphology by experienced personnel from US Fish and Wildlife Service. The fish was euthanized with an overdose of MS222. Remaining tissue and the DNA isolate were deposited in the Museum of Southwestern Biology (MSB Cat no. 117041 for \u003cem\u003eM. tetranema\u003c/em\u003e; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://msb.unm.edu/\u003c/span\u003e\u003cspan address=\"https://msb.unm.edu/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). A complete mitochondrial genome has previously been reported for this specimen (Osborne et al. \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e).\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003eDraft genome sequencing and assembly\u003c/h3\u003e\n\u003cp\u003eHigh-molecular-weight genomic DNA was isolated from muscle and fin tissue using genomic-tips 20/G (Qiagen\u0026reg;) according to the manufacturer\u0026rsquo;s directions. Prior to sequencing, DNA quality was assessed using Qubit HS assays (Invitrogen\u0026reg;, Thermo Fisher Scientific) and using pulse-field electrophoresis on a 1% agarose gel. A genomic library was prepared using PacBio HiFi SMRTbell\u0026reg; and the whole genome was sequenced using 1 SMRT cell on a PacBio Sequel II platform at the UC Davis Genomics Core facility.\u003c/p\u003e\u003cp\u003eTo estimate \u003cem\u003eM. tetranema\u003c/em\u003e genome size, Jellyfish v. 2.2.7 (Mar\u0026ccedil;ais and Kingsford \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2011\u003c/span\u003e) was used to obtain the k-mers count which was then used as input to GenomeScope v. 2.0 online tool (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://genomescope.org/genomescope2.0/\u003c/span\u003e\u003cspan address=\"http://genomescope.org/genomescope2.0/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e; Ranallo-Benavidez et al. \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). The genome was assembled from PacBio reads with the HiFiAsm assembler v.0.18.1-r466 (Cheng et al. \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) using the default parameters. Genome assembly was improved by scaffolding the contigs assembled by HiFiAsm to scaffold-level assemblies available for other leuciscid fishes using the homology-based scaffolding software RagTag (Alonge et al. \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Briefly, assembled contigs from \u003cem\u003eM. tetranema\u003c/em\u003e were mapped to primary genomic assemblies for Spikedace and Loach Minnow (\u003cem\u003eMeda fulgida\u003c/em\u003e and \u003cem\u003eTiaroga cobitis\u003c/em\u003e; Alexandre et al., \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2023\u003c/span\u003e) and the Eurasian Minnow (\u003cem\u003ePhoxinus phoxinus\u003c/em\u003e, Nunn et al., \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). The default aligner minimap2 (Li \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) was used and mapping quality score was increased to 20. Resulting assembly gap path (AGP) files were then merged with the edge weight set to 2 (e.g., scaffold joins had to be supported by at least 2 AGPs). The online platform CoGe (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://genomevolution.org/coge/\u003c/span\u003e\u003cspan address=\"https://genomevolution.org/coge/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) was used to obtain summary statistics from both contig and scaffold-level assemblies.\u003c/p\u003e\n\u003ch3\u003eNextRAD sequencing, variant call, quality filtering and microhaplotype data\u003c/h3\u003e\n\u003cp\u003eIsolation of DNA from tissue samples (2019 and 2020 collections) was performed using E.Z.N.A. tissue DNA (Omega Bio-Tek Inc.) or Zymo Quick-DNA (Zymo Research Corp.) kits according to the procedures outlined by the manufacturer. Isolations of DNA from previous genetic monitoring (2015\u0026ndash;2018; Osborne et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) were purified using Zymo DNA clean and concentrator kits (Zymo Research Corp.) to remove phenol which can inhibit DNA sequencing. All samples were treated with RNase A following DNA extraction. Each isolate was evaluated for the presence of high molecular weight DNA and the absence of RNA contamination by electrophoresis on a 1.2% agarose gel and by quantification using Qubit HS assays. One-hundred and ninety samples with sufficient high-quality DNA were sent to SNPsaurus, LLC (University of Oregon) for library preparation and sequencing. Genomic DNA was converted into a Nextera-tagmented reductively-amplified DNA (nextRAD) library and sequenced according to Russello et al. (\u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Genomic DNA was initially fragmented with Nextera DNA Flex reagent (Illumina\u0026reg;, Inc), which also ligates short adapter sequences to the ends of the fragments. The Nextera reaction was scaled for fragmenting 30 ng of genomic DNA, although 60 ng of genomic DNA was used for input to compensate for degraded DNA and to increase fragment sizes in samples. Fragmented DNA was then amplified for 27 cycles at 74 degrees, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence GTGTAGAGCC. Only fragments starting with a sequence that can be hybridized by the selective sequence of the primer will be efficiently amplified. All samples were pooled for the nextRAD library and sequenced with 150 base pair (bp) single-end reads on one Illumina Hi-Seq 4000 lane.\u003c/p\u003e\u003cp\u003eRaw sequence reads were received from the sequencing provider, trimmed for adapter sequences, and demultiplexed by individual and by lane. Raw reads were further trimmed to remove low quality bases with Trimmomatic v. 0.39 (Bolger et al. \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). Bases on both extremes of each read were removed if quality was bellow 20 or if ambiguous (N). Each read was then scanned with a 5-base wide sliding window, cutting when the average quality per base drops below 10. After trimming, reads were discarded if smaller than 60 bases.\u003c/p\u003e\u003cp\u003eRetained reads were mapped against the scaffold level draft genome with Bowtie v. 2.4.2 (Langmead and Salzberg \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2012\u003c/span\u003e) using the \u0026lsquo;local alignment\u0026rsquo; and default \u0026lsquo;very sensitive\u0026rsquo; options. Successfully aligned reads were filtered with SAMtools v. 1.16 (Li et al. \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Danecek et al. \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2021\u003c/span\u003e) to remove reads with mapping quality lower than 20. Before variant calling, we used Picard tools v. 2.26.2 (Broad Institute 2019) to add read group (RG) flags to bam files. Genetic variants were identified using FreeBayes v. 1.3.6 (Garrison and Marth \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). FreeBayes uses base quality scores to estimate a probability for each allele. We kept a maximum of 10 raw variants from each alignment with higher probabilities and at least a base quality of five.\u003c/p\u003e\u003cp\u003eTo remove erroneous or potentially erroneous variants, we used VCFtools v. 0.1.16 (Danecek et al. \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2011\u003c/span\u003e) to filter out variants with mean depth of coverage lower than 20 and higher than 200, minor allele count less than three, minor allele frequency lower than 2%, genotype depth of coverage lower than five, and with genotype quality lower than 20. Multi-nucleotide states were decomposed into single variants with vcflib (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/ekg/vcflib\u003c/span\u003e\u003cspan address=\"https://github.com/ekg/vcflib\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) and VCFtools was used to filter out nucleotide insertions and deletions and to retain only the bi-allelic SNPs. The dataset was then filtered by missing data, keeping SNPs present in at least 80% of samples and removing individuals with more than 30% missing data. After this step, SNPs were filtered using the bash script \u003cem\u003edDocent_filters\u003c/em\u003e (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/jpuritz/dDocent/blob/master/scripts/dDocent_filters\u003c/span\u003e\u003cspan address=\"https://github.com/jpuritz/dDocent/blob/master/scripts/dDocent_filters\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) that uses vcflib and VCFtools to filter loci based on allelic balance at heterozygous genotypes, strand representation, and quality vs depth of coverage.\u003c/p\u003e\u003cp\u003eFirst, loci were removed if at heterozygous positions, the alternative allele had a coverage lower than 20% or higher than 80% compared with the reference allele as reads with alleles from heterozygous positions are expected to have similar frequencies in the same individual. Alleles with frequencies smaller than 0.01 and higher that 0.99 were not removed to account for fixed alleles. Additionally, if the quality sum of the reference or alternative allele was zero, the locus was removed. This removes positions with spurious heterozygous genotype calls. Then loci with the ratio between the mean mapping quality of the alternative and reference allele lower than 0.25 or higher than 1.75 were removed, because loci from the same genomic location should not have large discrepancy between mapping qualities of two alleles. Furthermore, loci with quality scores less than half of the total depth were excluded because excessive depth inflates FreeBayes quality scores. Of the remaining loci, the average depth and standard deviation across all individuals was calculated. Loci with depth greater than the average depth plus one standard deviation were removed if the quality score was less than two times the depth. Finally, this script removed loci with a mean depth across individuals greater than two times the mode (98) that corresponded approximately to the 95th percentile of mean depth. Subsequently potential erroneous SNPs were filtered based on Hardy-Weinberg equilibrium (HWE) expectations with the pearl script \u003cem\u003efilter_hwe_by_pop.pl\u003c/em\u003e (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/jpuritz/dDocent/blob/master/scripts/filter_hwe_by_pop.pl\u003c/span\u003e\u003cspan address=\"https://github.com/jpuritz/dDocent/blob/master/scripts/filter_hwe_by_pop.pl\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). Typically, errors would have a low p-value and would be present in many populations; SNPs present in more that 50% of the populations (here each year was considered a \u0026lsquo;population\u0026rsquo;) and with an HWE p-value lower than 0.001 were removed. We filtered out potential incorrectly assembled paralogous loci that exhibited a large variation in read depth across all individuals. Standard deviation was estimated with package \u003cem\u003estats\u003c/em\u003e implemented in R v. 4.2.1 (R Core Team 2022) and read depth with VCFtools. Additional filtering based on missing data per locus (keeping loci present in 80% of individuals) was then reapplied.\u003c/p\u003e\u003cp\u003eRemaining SNPs were used to identify haplotypes within loci (referred to as microhaplotypes). Haplotyping SNPs within a locus also eliminates possible paralogous loci while neutralizing physical linkage without losing data (Willis et al. \u003cspan citationid=\"CR80\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). This was performed with the \u003cem\u003erad_haplotyper.pl\u003c/em\u003e pearl script (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/chollenbeck/rad_haplotyper\u003c/span\u003e\u003cspan address=\"https://github.com/chollenbeck/rad_haplotyper\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) excluding microhaplotypes if considered paralogs at least in five individuals and if missing from more than 30% of individuals.\u003c/p\u003e\u003cp\u003eRetained loci were tested for deviations from HWE and for linkage disequilibrium (LD) considering individuals captured in each year as a single \u0026lsquo;population\u0026rsquo;. Departures from HWE were accessed using a chi-square test on microhaplotype data with R package \u003cem\u003epegas\u003c/em\u003e v. 1.0 (Paradis \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2010\u003c/span\u003e) and using the Bonferroni correction for multiple comparisons implemented in the R package rcompanion v. 2.4.0 (Mangiafico \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2021\u003c/span\u003e), as implemented in the R function multi_HWE_tests (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/gcaeirodias/multi\u003c/span\u003e\u003cspan address=\"https://github.com/gcaeirodias/multi\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e HWE tests; Caeiro-Dias et al. \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Estimations of LD were performed on SNP data using the SNP of each microhaplotype with higher minimum allele frequency. If a SNP was found in LD, then the entire locus was removed. Tests for LD were performed using the chi-square test implemented in the R package \u003cem\u003eGUSLD\u003c/em\u003e v. 1.0.1 (Bilton et al. \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) and the Bonferroni correction to account for multiple simultaneous tests as implemented in the R pipeline significantLD (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/gcaeirodias/significantLD\u003c/span\u003e\u003cspan address=\"https://github.com/gcaeirodias/significantLD\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e; Caeiro-Dias et al. \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). If loci were found in multiple significant LD pairs, the loci that appeared in the highest number of comparisons were discarded to keep the maximum number of loci possible. In the remainder of instances, one locus from each pair were discarded randomly. Loci were considered as deviating from HWE and to be in LD if tests were significant across the six temporal samples (p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05). The resulting dataset should represent a robust genome-wide neutral SNP dataset.\u003c/p\u003e\u003cp\u003eAfter filtering, samples from one collection (TX 2017) had increased missing data when compared to other collections (see subsection \u003cem\u003eNextRAD sequencing, variant call, quality filtering and microhaplotype data\u003c/em\u003e of Results). Missing data can bias some results. As such, two datasets were created from the filtered set of markers and individuals described before. One dataset included all the microhaplotypes, but the TX 2017 collection was removed (NM dataset). For the other dataset, we removed the loci with more than 20% missing data in TX 2017 collection while retaining all the individuals after the filtering described above (NM-TX dataset). Unless otherwise specified, the analyses described below were performed with both datasets.\u003c/p\u003e\n\u003ch3\u003eTemporal genetic variation and evaluation of genetic drift\u003c/h3\u003e\n\u003cp\u003eTemporal genetic variation was first visualized using a discriminant analysis of principal components (DAPC), which summarizes genotypes in principal components to construct linear functions that maximize among group variation while minimizing within group variation. Analysis was performed using the R package adegenet v. 1.3\u0026ndash;1 (Jombart \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Jombart and Ahmed \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). Prior to DAPC we replaced missing data within each temporal sample using the Breiman's regression random forest algorithm (Breiman \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2001\u003c/span\u003e) implemented in R package randomForest v. 4.6\u0026ndash;14 (Liaw and Wiener \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2002\u003c/span\u003e). Values of missing data were predicted from 500 independently constructed regression trees and 50 bootstrap iterations with default bootstrap sample size. This was preferred over the default \u0026ldquo;mean method\u0026rdquo; (i.e., missing genotypes are replaced by the average estimated across the data set) implemented in adegenet to ensure that we did not artificially increase similarity of allele frequencies across temporal samples. An initial DAPC was performed using temporal samples as groups, allele frequencies centered but not scaled, retaining all PCA and DA axes, and keeping other options as default. The \u003cem\u003ea-score\u003c/em\u003e method was used to select the optimal number of principal components to retain for the final DAPC, using the maximum number of PCs, all DAs and the other options were set as default. The final DAPC was performed using the optimal number of PCs, two DAs and keeping the other default options.\u003c/p\u003e\u003cp\u003eWe assessed genetic differentiation over time by estimating pairwise F\u003csub\u003eST\u003c/sub\u003e between all temporal collections and p-values using 1,000 bootstrap iterations over loci as implemented in GenoDive v. 3.06 (Meirmans \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Significant values of F\u003csub\u003eST\u003c/sub\u003e between consecutive years or those separated by a few years, reflect differences in allele frequencies and can be a sign of genetic drift. Next, we directly quantified allele frequency changes in each locus across time by calculating absolute allele frequency differences (AFD) between consecutive years and identified outliers of frequency changes in each pair of years. This analysis was conducted to evaluate the proportion of the loci with inflated changes in allele frequency. When using the NM-TX dataset we did not include TX 2017 because there were no consecutive years to estimate allele frequency differences for this location. Allele frequencies were estimated with GenoDive and for each locus, allele frequency differences were averaged across alleles.\u003c/p\u003e\n\u003ch3\u003eTemporal genetic diversity, inbreeding, and relatedness\u003c/h3\u003e\n\u003cp\u003eFor each temporal collection, we estimated standard genetic diversity and inbreeding metrics. We used R package diverRsity v. 1.9.90 (Keenan et al. \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2013\u003c/span\u003e) to estimate expected heterozygosity (HS; Nei \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e1987\u003c/span\u003e), observed heterozygosity (H\u003csub\u003eO\u003c/sub\u003e), allelic richness (A\u003csub\u003eR\u003c/sub\u003e), and inbreeding coefficient (F\u003csub\u003eIS\u003c/sub\u003e). All metrics were estimated for each population and by locus within each population. If changes to genomic diversity occurred across time, we expect to see significant changes in the distribution of values of those metrics across temporal collections. To evaluate this, we tested if the distribution of values for each metric estimated with both datasets were similar across populations using non-parametric Kruskal-Wallis tests, because diversity measures did not follow a normal distribution and/or variance across populations was not homogeneous (data not shown). When tests were significant, we performed pairwise Wilcoxon tests between populations using a Bonferroni correction for multiple comparisons. All tests were performed with the built-in R package stats.\u003c/p\u003e\u003cp\u003eMaximum likelihood estimation (MLE) of identity-by-descent (IBD) for each pair of individuals was estimated with the R package SNPrelate v. 1.32.2 (Zheng et al. \u003cspan citationid=\"CR82\" class=\"CitationRef\"\u003e2012\u003c/span\u003e) using the \u003cem\u003eExpectation-Maximization\u003c/em\u003e algorithm. Because SNPrelate does not take loci with more than two alleles as input, we selected one SNP at random for each locus. SNPs with more than 5% missing data were excluded as it can have an impact on relatedness analyses. Resulting MLE-IBD coefficients (k0 and k1) were used to construct k0,k1 plots (e.g., Galv\u0026aacute;n-Femen\u0026iacute;a et al. \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2017\u003c/span\u003e), where k0 and k1 are the probabilities that two individuals share zero or one IBD alleles, respectively. This analysis was performed to evaluate if the overall relatedness changed overtime.\u003c/p\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003eTemporal trends in effective population size\u003c/h2\u003e\u003cp\u003eLinkage disequilibrium effective size (LD N\u003csub\u003ee\u003c/sub\u003e) was estimated from both datasets for every year using the method from Hill (\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e1981\u003c/span\u003e) implemented in NeEstimator v. 2.1 (Do et al. \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) excluding microhaplotypes with allele frequencies lower than 2%. This metric reflects the effective number of parents that produced the progeny from which the sample was drawn (see Waples \u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e2005\u003c/span\u003e). To estimate the 95% confidence intervals (CIs) the jackknife approach may be more suitable for data sets with large numbers of loci (Jones et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2016\u003c/span\u003e) as loci may not be entirely independent (Gilbert and Whitlock \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). However, datasets in this study were filtered to avoid linkage and indeed linkage among loci was limited as r\u003csup\u003e2\u003c/sup\u003e per population with more than 20 samples was always lower than 0.06. Hence, a parametric approach was used to estimate the 95% CIs.\u003c/p\u003e\u003cp\u003eFor consecutive temporal samples, variance effective population size (N\u003csub\u003eeV\u003c/sub\u003e) and 95% CIs were estimated using the temporal method of Nei and Tajima (\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e1981\u003c/span\u003e) as implemented in NeEstimator v. 2.1, using the Plan I option (sampling with replacement), excluding microhaplotyes with allele frequencies lower than 2%, as highly polymorphic loci with many rare alleles can result in biased estimates of N\u003csub\u003eeV\u003c/sub\u003e (Hedrick 1999; Turner et al. 2001). The parametric approach was used to estimate the 95% CIs. To estimate N\u003csub\u003eeV\u003c/sub\u003e we also used both datasets, but TX 2017 was excluded from NM-TX dataset as we only had a single temporal sample from TX.\u003c/p\u003e\u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e\u003ch2\u003eDraft genome sequencing\u003c/h2\u003e\u003cp\u003eA total of 2,244,201 raw PacBio reads were generated, containing 29.72 Giga bases (Gb) of sequence data. Genome size was estimated from the k-mers count to be 1.04 Gb. Contig assembly resulted in 960 sequences, ranging from 6,390 to 38,396,035 bp (mean\u0026thinsp;=\u0026thinsp;1,331,581, mode\u0026thinsp;=\u0026thinsp;30,561) with a total length of 1,259,117,741 bp (1.26 Gb) and an N50 of 18,070,155 bp. After scaffolding the number of sequences decreased to 923, ranging from 6,390 to 48,550,500 bp (mean\u0026thinsp;=\u0026thinsp;1,364,162, mode\u0026thinsp;=\u0026thinsp;29,197) with a total length of 1,259,121,441 bp and an N50 of 24,795,794 bp. This corresponds to a 23.6x genome coverage if we assume that the total assembly size is similar to the true genome size, which is expected to be about 1.2 Gb based on knowledge on genomes of other related species (e.g., Alexandre et al. \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2023\u003c/span\u003e).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003eNextRAD sequencing, variant call, quality filtering and microhaplotype data\u003c/h2\u003e\u003cp\u003eAfter trimming, read alignment to the draft genome yielded an average of 3.5\u0026nbsp;million (M) reads per individual (minimum\u0026thinsp;=\u0026thinsp;0.1 M; maximum\u0026thinsp;=\u0026thinsp;5 M). FreeBayes identified 1.6 M raw variants (including SNPs, multi-nucleotide polymorphisms, indels and other complex variants) across the 190 individuals. After filtering, a total of 2,804 loci containing 6,725 SNPs across 187 individuals with a maximum of 30% missing data. Average depth per locus and per individual was 46.1 (ranging from 20.4 to 98.5 and 13.4 to 79.8, respectively). Missing data estimates across all loci was on average 0.14 for TX 2017 while for other collections ranged from 0.02 to 0.06 (Figure \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e, Supplementary material). From this dataset two others were generated. The NM dataset included all 2,804 loci and 179 individuals after TX 2017 was excluded. The NM-TX dataset was restricted to 1,908 loci with less than 20% missing data and all 187 individuals. Dataset details are presented in Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e.\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003e\u003cb\u003e\u0026ndash;\u003c/b\u003e Number of loci and samples from each temporal collection retained in each analyzed dataset, after bioinformatics filtering.\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"9\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c8\" colnum=\"8\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c9\" colnum=\"9\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u0026nbsp;\u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/th\u003e\u003cth align=\"left\" colspan=\"6\" nameend=\"c9\" namest=\"c4\"\u003e\u003cp\u003eYear\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e2015\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e2016\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e\u003cp\u003e2017\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e2018\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e2019\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e2020\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\" morerows=\"2\" rowspan=\"3\"\u003e\u003cp\u003e\u003cb\u003eDataset\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eNM (2,804 loci)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eUpstream\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e34\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e31\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e\u003cp\u003e21\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e41\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e43\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\" morerows=\"1\" rowspan=\"2\"\u003e\u003cp\u003eNM-TX (1,908 loci)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eUpstream\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e34\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e31\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e\u003cp\u003e21\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e41\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e43\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eDownstream\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c6\"\u003e\u003cp\u003e8\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003ePatterns of temporal genomic variation\u003c/h2\u003e\u003cp\u003eRegardless of the dataset, a considerable number of samples overlapped across the DAPC space (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), however the DAPC performed with the NM-TX dataset (1,908 loci) shows less variability among NM samples (Figs.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea and \u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb) than the DAPC performed with the NM dataset (2,804 loci; Figs.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ec and \u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ed). Using all 2,804 loci, NM samples largely overlapped from 2015 until 2018, but most of the samples showed a shift in genetic variation in 2019 and a further shift in 2020.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eThe estimated pairwise F\u003csub\u003eST\u003c/sub\u003e values between all possible comparisons and using fewer loci (NM-TX dataset) were very low (\u0026le;0.003) and mostly not statistically different from zero (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea). The exceptions were three comparisons involving mostly NM 2016 and NM 2020. When using the NM dataset, the magnitude of the F\u003csub\u003eST\u003c/sub\u003e values was similar (\u0026le;0.002) but eight comparisons were significant, including the consecutive years 2015\u0026ndash;2016, 2016\u0026ndash;2017, and 2019\u0026ndash;2020. Five comparisons involved NM 2020, and four included NM 2016 (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb). Regardless of the dataset used, NM 2016 and NM 2020 had more significant F\u003csub\u003eST\u003c/sub\u003e differences when compared with other temporal/spatial samples. However, it is worth mentioning that the significance of F\u003csub\u003eST\u003c/sub\u003e is impacted by low sample size, such as in TX 2017 and NM 2018.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eMean AFD between consecutive years using both datasets showed overall small changes (ranging from 0.04 to 0.07 in both cases), but some loci exhibited substantial fluctuations (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). Overall, more loci and bigger changes occurred between 2017\u0026ndash;2018 and 2018\u0026ndash;2019 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e and Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), but those values might be biased due to small sample size in 2018. However, excluding those comparisons, frequency changes of 0.1 or higher were noted in 79 to 173 loci with NM-TX dataset and 117 to 267 loci with NM dataset (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). Moreover, 76 to 86 loci were outliers in terms of AFD between consecutive years when using the NM-TX dataset and 105 to 130 when using the NM dataset (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). Most of the outlier loci were shared between datasets (Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). Moreover, regardless of the dataset, only two loci located in different contigs showed AFD of 0.1 or higher across all consecutive years. However, allele frequency changes do not show any particular trend (Figure S2, Supplementary material).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab2\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003e\u003cb\u003e\u0026ndash;\u003c/b\u003e Number of loci with absolute allele frequency differences (AFD) between consecutive years distributed in five classes of 0.1. The highest absolute allele frequency difference between two consecutive years was 0.41.\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"7\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u003cp\u003eDataset\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e\u003cp\u003eAFD Class\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e\u003cp\u003e2015\u0026ndash;2016\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c4\"\u003e\u003cp\u003e2016\u0026ndash;2017\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c5\"\u003e\u003cp\u003e2017\u0026ndash;2018\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c6\"\u003e\u003cp\u003e2018\u0026ndash;2019\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c7\"\u003e\u003cp\u003e2019\u0026ndash;2020\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\" morerows=\"5\" rowspan=\"6\"\u003e\u003cp\u003eNM (2,804 loci)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0,0.1[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e2605 (92.9%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e2537 (90.48%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e2236 (79.74%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e2296 (81.9%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e2687 (95.83%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.1,0.2[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e192 (6.85%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e251 (8.95%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e489 (17.44%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e452 (16.1%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e116 (4.13%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.2,0.3[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e6 (0.21%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e16 (0.57%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e66 (2.35%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e48 (1.7%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e1 (0.04%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.3,0.4[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e1 (0.04%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e12 (0.43%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e8 (0.3%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.4,0.5[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e1 (0.04%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eOutliers\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e105 (3.74%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e105 (3.74%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e130 (4.64%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e126 (4.49%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e124 (4.42%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\" morerows=\"5\" rowspan=\"6\"\u003e\u003cp\u003eNM-TX (1,908 loci)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0,0.1[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e1763 (92.4%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e1735 (91%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e1536 (80.5%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e1578 (82.7%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e1829 (95.9%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.1,0.2[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e141 (7.4%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e165 (8.6%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e324 (16.98%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e285 (14.94%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e79 (4.1%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.2,0.3[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e4 (0.2%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e8 (0.4%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e38 (2%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e38 (1.99%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.3,0.4[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e10 (0.52%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e7 (0.37%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e[0.4,0.5[\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e0\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eOutliers\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e76 (3.98%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e64 (3.35%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e92 (4.82%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e88 (4.61%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e79 (4.14%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colspan=\"2\" nameend=\"c2\" namest=\"c1\"\u003e\u003cp\u003eOutliers common to both datasets\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e76\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e63\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e86\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e83\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e79\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003eTemporal genetic diversity, inbreeding, and relatedness\u003c/h2\u003e\u003cp\u003eGenetic diversity metrics (A\u003csub\u003eR\u003c/sub\u003e, H\u003csub\u003eO\u003c/sub\u003e, and H\u003csub\u003eS\u003c/sub\u003e) estimated with both datasets were very similar across time and generally not significantly different between temporal collections (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea\u0026ndash;c and \u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ee\u0026ndash;g). However, Kruskal-Wallis tests were significant for A\u003csub\u003eR\u003c/sub\u003e estimated with both dataset (NM-TX dataset p-value\u0026thinsp;=\u0026thinsp;7.87x10\u003csup\u003e\u0026minus;\u0026thinsp;18\u003c/sup\u003e, NM dataset p-value\u0026thinsp;=\u0026thinsp;1.2x10\u003csup\u003e\u0026minus;\u0026thinsp;11\u003c/sup\u003e; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea and \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed) and for H\u003csub\u003eO\u003c/sub\u003e estimated with the NM dataset (p-value\u0026thinsp;=\u0026thinsp;0.02; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ee). Pairwise Wilcoxon tests revealed that A\u003csub\u003eR\u003c/sub\u003e values were significantly different when comparing TX 2017 or NM 2018 to other populations (Tables S1 and S2, Supplementary material) and revealed H\u003csub\u003eO\u003c/sub\u003e as significantly different between NM 2018 and 2020 (Tables S3, Supplementary material). In both cases, genetic diversity metrics were significant when one of the populations had smaller sample size suggesting that sample size was a possible source of bias.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eAlthough the average inbreeding coefficients estimated with both datasets were close to zero (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed and \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eh), most pairwise comparisons were significantly different (Tables S4 and S5, Supplementary material). In most cases differences between collections were subtle (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed and \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eh), but an obvious lower genome-wide F\u003csub\u003eIS\u003c/sub\u003e was found in TX 2017 (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed), which had the lowest Wilcoxon test p-values with all other temporal collections (except when compared to NM 2018; Tables S4, Supplementary material). These results can be indicative of an excess of heterozygotes in many loci.\u003c/p\u003e\u003cp\u003eRelatedness analyses were performed with a single SNP per locus from the NM-TX and NM datasets. Removing SNPs with missing data higher than 5% resulted in 1,647 and 2,136 SNPs for the NM-TX and NM datasets, respectively. Results with both datasets were very similar (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). In general, sampled individuals were unrelated to each other (theoretical expectations for totally unrelated individuals are k0 \u0026asymp; 1 and k1 \u0026asymp; 0). A few pairs exhibit a probability of shared alleles similar to the expectations for first cousins (k0\u0026thinsp;=\u0026thinsp;0.75, k1\u0026thinsp;=\u0026thinsp;0.25) but values of k0\u0026thinsp;\u0026gt;\u0026thinsp;0.75 and k1\u0026thinsp;\u0026lt;\u0026thinsp;0.25 may also include more distant relatives. Two individuals exhibited probabilities of shared alleles between what is expected for first cousins and second-degree relationships (k0\u0026thinsp;=\u0026thinsp;0.5, k1\u0026thinsp;=\u0026thinsp;0.5; half-siblings, avuncular, grandchild-grandparent). Among those few pairs of potential relatives (from second degree to distant relatives), highest relatedness was between an individual from NM 2017 and an individual from NM 2020 (k0\u0026thinsp;=\u0026thinsp;0.57, k1\u0026thinsp;=\u0026thinsp;0.43 with NM-TX dataset; k0\u0026thinsp;=\u0026thinsp;0.62, k1\u0026thinsp;=\u0026thinsp;0.34 with NM dataset;) but all the other pairs with k0 \u0026le; 0.8 were from NM 2015 and or NM 2016. In addition, two other pairs of individuals from NM 2019 exhibited k0 \u0026asymp; 0 and k1 \u0026asymp; 0.12, meaning that the probability of sharing zero alleles by descent was close to 0% and the probability of sharing one allele was about 12%. Because this analysis used biallelic data, the probability of sharing two alleles by descent was around 88%, which means that those two pairs were either identical twins or, most likely, resulted from matings between related individuals.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\u003ch2\u003eTemporal trends in effective population size\u003c/h2\u003e\u003cp\u003eEstimates of LD N\u003csub\u003ee\u003c/sub\u003e and 95% CIs in NM were consistent between both datasets (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003ea and Table S6, Supplementary material). This estimator suggests that the population in New Mexico increased from 2014 (LD N\u003csub\u003ee\u003c/sub\u003e = 2,255, estimated from 2015 sample) to 2015 (LD N\u003csub\u003ee\u003c/sub\u003e = 3,666, estimated from 2016 sample). Whether the trend was maintained until 2017 is not clear because point estimates and CIs of LD N\u003csub\u003ee\u003c/sub\u003e for 2016 and 2017 parental generations were infinite. Infinite values mean that LD N\u003csub\u003ee\u003c/sub\u003e values were very high or alternatively the power to estimate LD N\u003csub\u003ee\u003c/sub\u003e was low due to small sample sizes. In 2018, LD N\u003csub\u003ee\u003c/sub\u003e declined drastically (estimated from the 2019 sample) followed by a notable increase in 2019 (estimated from the 2020 sample).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eEstimates of N\u003csub\u003eeV\u003c/sub\u003e in NM obtained with both datasets were very small (\u0026lt;\u0026thinsp;351; Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eb and Table S7, Supplementary material). Only the point estimate obtained with NM-TX dataset between 2017 and 2018 was infinite (this is also likely to because of the small 2018 sample size). Estimates obtained with NM-TX dataset were similar across time-series while the NM dataset suggested a small increase in 2017\u0026ndash;2018. Nevertheless, the confidence intervals obtained from both datasets largely overlapped, suggesting that N\u003csub\u003eeV\u003c/sub\u003e never exceeded 300, except between 2017 and 2018 (NM dataset N\u003csub\u003eeV\u003c/sub\u003e = 351; CI = [196, 1,374]). However, the sample size for those two years were the smallest and these values should be interpreted with caution.\u003c/p\u003e\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eIn this study, temporal trends of genomic variation, genetic diversity metrics and effective population size in \u003cem\u003eM. tetranema\u003c/em\u003e were evaluated from samples collected from 2015 to 2020. Results from genomic variation and N\u003csub\u003ee\u003c/sub\u003e across time suggests that \u003cem\u003eM. tetranema\u003c/em\u003e upstream (NM) is affected by genetic drift, but most genetic diversity metrics did not change significantly across the time-series. Moreover, relatedness metrics suggest that the large majority of the breeding populations upstream are unrelated.\u003c/p\u003e\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\u003ch2\u003eEvidence of upstream genetic drift from temporal genomic variation and effective population size\u003c/h2\u003e\u003cp\u003eThe DAPC results using different datasets exhibited some differences. The NM dataset (2,804 loci) samples from 2019 and 2020 showed a higher dispersion and were shifted compared to all other samples while the NM-TX dataset (1908 loci) detected more subtle shifts. Notably, the dataset with more loci (NM dataset) revealed some changes in genomic variation not detected with the reduced dataset (NM-TX dataset). Consistent with these results, most pairwise F\u003csub\u003eST\u003c/sub\u003e estimates were significantly different from zero when analyzing the NM dataset but not when using the NM-TX dataset. However, outlier AFD loci across time (i.e., loci with higher differences) were essentially the same in both datasets. This result suggests that the differences observed in the DAPC and F\u003csub\u003eST\u003c/sub\u003e between datasets is likely a result of many loci with moderate AFD values present in the NM dataset, particularly after 2018, rather than few outlier loci (i.e., higher frequency changes). Moreover, regardless of the magnitude of AFD no particular trend was identified, i.e., changes were random. The moderate to large values of AFD detected in our datasets are expected when the population experiences genetic drift (Gompert et al. \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2021\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eAcross time-series, N\u003csub\u003eeV\u003c/sub\u003e estimates were similar and consistently low, regardless of the dataset used, suggesting substantial changes in allele frequencies between consecutive years. These results were very similar to those obtained with microsatellite data (Osborne et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Low N\u003csub\u003eeV\u003c/sub\u003e estimates are a sign that genetic drift is currently a strong evolutionary force (Nei and Tajima \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e1981\u003c/span\u003e; Jorde and Ryman \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e1995\u003c/span\u003e; Gompert et al. \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Osborne et al. \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e). Unlike N\u003csub\u003eeV\u003c/sub\u003e, estimates of LD N\u003csub\u003ee\u003c/sub\u003e in New Mexico (upstream) showed a small increase in N\u003csub\u003ee\u003c/sub\u003e from 2014 to 2015. These results are also consistent with microsatellite data that showed a progressive increase in LD N\u003csub\u003ee\u003c/sub\u003e between 2014 and 2017 in NM (Osborne et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). However, the values reported here are higher and with finite confidence intervals (e.g., LD N\u003csub\u003ee\u003c/sub\u003e [NM dataset 2014]\u0026thinsp;=\u0026thinsp;2,255.3, CI = [1,917.5; 2,736.5]; LD N\u003csub\u003ee\u003c/sub\u003e [microsatellites 2014]\u0026thinsp;=\u0026thinsp;1,386, CI = [265; \u0026infin;]) due to the larger number of loci. Improved precision is an advantage offered by genomic time-series data when compared to microsatellite temporal data (Osborne et al. \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e). However, it is worth noting that when analyzing a high number of loci and relatively small number of individuals, estimates of N\u003csub\u003ee\u003c/sub\u003e based on linkage disequilibrium may be inflated (Waples and Do \u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Ragsdale and Gravel \u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Following the gradual increase, LD N\u003csub\u003ee\u003c/sub\u003e suffered an accentuated decline in 2018 (estimated from 2019 sample) and in the following year it was the highest recorded across the entire time-series. While the general increasing trend of LD N\u003csub\u003ee\u003c/sub\u003e described here is consistent with relative abundance data (catch-per-unit-effort, CPUE) obtained between 2012 and 2019 for the same population (Osborne et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e), the LD N\u003csub\u003ee\u003c/sub\u003e decrease in 2018 was contrary to the increased CPUE of \u003cem\u003eM. tetranema\u003c/em\u003e over that period. This suggests that despite increasing abundance, in some years, reproductive success may be limited to a small fraction of the population. For example, this can happen if the population is not entirely panmictic, contrary to what has generally been assumed (e.g., Osborne et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Also, the LD N\u003csub\u003ee\u003c/sub\u003e decrease in 2018 is in line with the biggest changes detected with DAPC after 2018. Moreover, LD N\u003csub\u003ee\u003c/sub\u003e fluctuations are expected in pelagophilic fish species inhabiting desert rivers due to the highly variable environmental conditions (Osborne et al. \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2012\u003c/span\u003e, \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e, \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eOverall, the results presented here demonstrate that genetic drift is an overarching evolutionary force driving trends in genetic variation in \u003cem\u003eM. tetranema\u003c/em\u003e. Its predominance is explained by ecological and demographic factors. For example, a regional drought between 2011 and 2013 resulted in extirpation of the population from the Ninnescah and Arkansas rivers in Kansas (Perkin et al. \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Pennock et al. \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2017\u003c/span\u003e) and led to decreases in \u003cem\u003eM. tetranema\u003c/em\u003e abundance and distribution in the South Canadian River with lowest abundance occurring in 2012 (Pennock et al. \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Osborne et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Low abundance likely increased the effect of genetic drift in the population. Yet, population recovery after 2012 did not reduce the influence of genetic drift as suggested by the results presented here. In fact, some analysis suggest that genetic drift might have increased in more recent years despite increased abundance. Genetic drift and differential reproductive success (i.e., non-panmictic population) in the upstream portion of the species\u0026rsquo; range maybe driven by downstream biased gene flow due to egg and larval drift, but this prediction lacks explicit testing in the current study.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec17\" class=\"Section2\"\u003e\u003ch2\u003eMaintenance of genetic diversity in face of genetic drift\u003c/h2\u003e\u003cp\u003eDespite signs of genetic drift, significant changes were not observed in most genetic diversity metrics between temporal collections. Exceptions were A\u003csub\u003eR\u003c/sub\u003e and H\u003csub\u003eO\u003c/sub\u003e in a single instance when sample size of one of the compared collections was small. Estimates of H\u003csub\u003eS\u003c/sub\u003e and generally H\u003csub\u003eO\u003c/sub\u003e were more robust to bias introduced by small sample sizes while A\u003csub\u003eR\u003c/sub\u003e was more sensitive as expected (Leberg \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2002\u003c/span\u003e). Regarding F\u003csub\u003eIS\u003c/sub\u003e, many distributions from both datasets are statistically different, but in general follow similar patterns with F\u003csub\u003eIS\u003c/sub\u003e close to zero for most of the loci suggesting that most \u003cem\u003eM. tetranema\u003c/em\u003e genome does not experience excesses of heterozygosity or homozygosity across time. Despite some significant tests likely due to small sample sizes, the results presented here suggest that genetic diversity in NM across time remained similar.\u003c/p\u003e\u003cp\u003eGenetic drift and maintenance of genetic diversity in upstream might be related with movement ecology in this species. Previous field observations reported high numbers of \u003cem\u003eM. tetranema\u003c/em\u003e larvae and young-of-the-year downstream while upstream population is dominated by adults (S. Davenport and J. Hatt, personal communication). While there\u0026rsquo;s a lack of robust data, such observation might be related with downstream drift of fertilized eggs while developing (Bottrell et al. \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e1964\u003c/span\u003e) if eggs and larvae retention in upstream nursery habitats is limited. Upstream spawning migration by larger, sexually mature \u003cem\u003eM. tetranema\u003c/em\u003e prior and during reproductive season was previously suggested (Bonner \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2000\u003c/span\u003e), and a recent mark-recapture study found that adults exhibit upstream-biased movement in the South Canadian River (Steffensmeier et al. \u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Together, genetic and ecological data suggest that \u003cem\u003eM. tetranema\u003c/em\u003e upstream experiences genetic drift likely due of downstream-biased gene flow, while upstream movement of adults may explain low relatedness and maintenance of genetic diversity upstream in the face of genetic drift.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec18\" class=\"Section2\"\u003e\u003ch2\u003eConservation implications\u003c/h2\u003e\u003cp\u003eAlthough the results presented in this study do not suggest immediate detrimental effects on genomic diversity, abrupt demographic changes, particularly in human-mediated environmental changes is a real threat in fish with similar life history traits to \u003cem\u003eM. tetranema\u003c/em\u003e (e.g., Osborne et al. \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2023a\u003c/span\u003e). Moreover, the maintenance of genetic diversity and N\u003csub\u003ee\u003c/sub\u003e upstream might be linked to upstream movement of adults. In fragmented habitats, another pelagophilic species experienced genetic diversity and N\u003csub\u003ee\u003c/sub\u003e declines upstream with persistence sustained by population augmentation (Osborne et al. In Review). Maintaining connectivity in the South Canadian River where \u003cem\u003eM. tetranema\u003c/em\u003e persists should be a priority.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec19\" class=\"Section2\"\u003e\u003ch2\u003eFuture research\u003c/h2\u003e\u003cp\u003eTemporal sampling upstream and downstream with increased sample size is needed to test the hypotheses related with asymmetric gene flow and distribution of genomic diversity across the riverscape. This would allow more robust comparisons of genetic variability and the effect of genetic drift across time and space and would also allow estimation of the magnitude and direction of movement between up and downstream sites. Moreover, the transition from microsatellite-based genetic monitoring to SNP-based assessments will be an important component of future conservation and management efforts. This requires the development of a panel of markers that can be consistently monitored across time. A SNP-based panel can be used for both long-term genetic monitoring and to obtain genetic data to test various hypotheses. Genotyping-in-Thousands by sequencing (GT-seq, Campbell et al. \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2015\u003c/span\u003e) is an efficient and cost-effective method of targeted SNP genotyping that uses multiplexed PCR amplicon sequencing and facilitates simultaneous amplification of hundreds of targeted genetic loci in hundreds to thousands of samples. The SNP and microhaplotype data developed in this study offer an excellent opportunity to identify and select a set of loci that can be used to achieve those research goals.\u003c/p\u003e\u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe sincerely thank Joanna Hatt, Andrew Monie, John Caldwell, Eliza Gilbert (New Mexico Department of Game and Fish), David Kitcheyan, Stephen Davenport, and Daniel Fenner (U.S. Fish and Wildlife Service) for samples collection. We thank Emily DeArmon and Alexandra Snyder (UNM Museum of Southwestern Biology Division of Fishes) for expert curatorial services. High performance computing was conducted at the UNM Center for Advanced Research Computing that is supported, in part, by the National Science Foundation. This study was funded by the New Mexico Department of Game and Fish. Founded in 1889, the University of New Mexico sits on the traditional homelands of the Pueblo of Sandia. The original peoples of New Mexico \u0026ndash; Pueblo, Navajo, and Apache \u0026ndash; since time immemorial, have deep connections to the land and have made significant contributions to the broader community statewide. We honor the land itself and those who remain stewards of this land throughout the generations and also acknowledge our committed relationship to Indigenous peoples. We gratefully recognize our history.\u003c/p\u003e\u003cp\u003eFunding\u003c/p\u003e\n\u003cp\u003eThis study was funded by the New Mexico Department of Game and Fish.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;Competing Interests\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eThe authors declare no conflicts of interest.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;Author Contributions\u003c/p\u003e\n\u003cp\u003eG.C.D. was responsible for project design, bioinformatics, data analysis, and writing the manuscript. A.C. was responsible for bioinformatics and critical revision of the manuscript, T.T. was responsible for critical revision of the manuscript. M.O. was responsible for project conception, design and coordination, laboratory work, critical revision of the manuscript, and obtaining funding.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eData Availability\u003c/p\u003e\n\u003cp\u003eIndividual genotype data from microhaplotypes are available from authors upon request. Raw sequence reads from nextRAD-seq (190 individuals), are deposited in the NCBI Sequence Read Archive (SRA), with the BioProject accession number XXXXXX (to be added upon acceptance).\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAlexandre NM, Cameron AC, Tian D, et al (2023) Chromosome-level reference genomes of two imperiled desert fishes: spikedace (\u003cem\u003eMeda fulgida\u003c/em\u003e) and loach minnow (\u003cem\u003eTiaroga cobitis\u003c/em\u003e). G3: Genes, Genomes, Genetics 13:jkad157\u003c/li\u003e\n\u003cli\u003eAlonge M, Lebeigle L, Kirsche M, et al (2022) Automated assembly scaffolding using RagTag elevates a new tomato system for high-throughput genome editing. Genome biology 23:258\u003c/li\u003e\n\u003cli\u003eAlp M, Keller I, Westram AM, Robinson CT (2012) How river structure and biological traits influence gene flow: a population genetic study of two stream invertebrates with differing dispersal abilities. Freshwater Biology 57:969\u0026ndash;981\u003c/li\u003e\n\u003cli\u003eArchdeacon TP, Davenport SR, Grant JD, Henry EB (2018) Mass upstream dispersal of pelagic-broadcast spawning cyprinids in the Rio Grande and Pecos River, New Mexico. Western North American Naturalist 78:100\u0026ndash;105\u003c/li\u003e\n\u003cli\u003eArchdeacon TP, Diver-Franssen TA, Bertrand NG, Grant JD (2020) Drought results in recruitment failure of Rio Grande silvery minnow (\u003cem\u003eHybognathus amarus\u003c/em\u003e), an imperiled, pelagic broadcast-spawning minnow. Environmental Biology of Fishes 103:1033\u0026ndash;1044\u003c/li\u003e\n\u003cli\u003eBalon EK (1975) Reproductive guilds of fishes: a proposal and definition. Journal of the Fisheries Board of Canada 32:821\u0026ndash;864\u003c/li\u003e\n\u003cli\u003eBalon EK (1981) Additions and amendments to the classification of reproductive styles in fishes. Environmental Biology of Fishes 6:377\u0026ndash;389\u003c/li\u003e\n\u003cli\u003eBergland AO, Behrman EL, O\u0026rsquo;Brien KR, et al (2014) Genomic evidence of rapid and stable adaptive oscillations over seasonal time scales in \u003cem\u003eDrosophila\u003c/em\u003e. PLoS genetics 10:e1004775\u003c/li\u003e\n\u003cli\u003eBilton TP, McEwan JC, Clarke SM, et al (2018) Linkage disequilibrium estimation in low coverage high-throughput sequencing data. Genetics 209:389\u0026ndash;400\u003c/li\u003e\n\u003cli\u003eBlanchet S, Prunier JG, Paz‐Vinas I, et al (2020) A river runs through it: The causes, consequences, and management of intraspecific diversity in river networks. Evolutionary Applications 13:1195\u0026ndash;1213\u003c/li\u003e\n\u003cli\u003eBolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114\u0026ndash;2120\u003c/li\u003e\n\u003cli\u003eBonner TH (2000) Life history and reproductive ecology of the Arkansas River shiner and peppered chub in the Canadian River, Texas and New Mexico. Texas Tech University\u003c/li\u003e\n\u003cli\u003eBonner TH, Wilde GR (2000) Changes in the fish assemblage of the Canadian River, Texas, associated with reservoir construction. Journal of Freshwater Ecology 15:189\u0026ndash;198\u003c/li\u003e\n\u003cli\u003eBottrell CE, Ingersol RH, Jones RW (1964) Notes on the embryology, early development, and behavior of \u003cem\u003eHybopsis aestivalis tetranemus\u003c/em\u003e (Gilbert). Transactions of the American Microscopical Society 83:391\u0026ndash;399\u003c/li\u003e\n\u003cli\u003eBreiman L (2001) Random forests. Machine Learning 45:5\u0026ndash;32\u003c/li\u003e\n\u003cli\u003eCaeiro-Dias G, Osborne M, Turner T (2024) Time is of the essence: using archived samples in the development a GT-seq panel to preserve continuity of ongoing genetic monitoring. \u003cem\u003eAuthorea\u003c/em\u003e.\u003c/li\u003e\n\u003cli\u003eCampbell NR, Harmon SA, Narum SR (2015) Genotyping‐in‐Thousands by sequencing (GT‐seq): A cost effective SNP genotyping method based on custom amplicon sequencing. Molecular Ecology Resources 15:855\u0026ndash;867\u003c/li\u003e\n\u003cli\u003eChavez MJ, Budy P, Pennock CA, et al (2024) Movement patterns of a small-bodied minnow suggest nomadism in a fragmented, desert river. Movement Ecology 12:52\u003c/li\u003e\n\u003cli\u003eCheng H, Concepcion GT, Feng X, et al (2021) Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nature Methods 18:170\u0026ndash;175\u003c/li\u003e\n\u003cli\u003eDanecek P, Auton A, Abecasis G, et al (2011) The variant call format and VCFtools. Bioinformatics 27:2156\u0026ndash;2158\u003c/li\u003e\n\u003cli\u003eDanecek P, Bonfield JK, Liddle J, et al (2021) Twelve years of SAMtools and BCFtools. GigaScience 10:giab008\u003c/li\u003e\n\u003cli\u003eDemandt MH (2010) Temporal changes in genetic diversity of isolated populations of perch and roach. Conservation Genetics 11:249\u0026ndash;255\u003c/li\u003e\n\u003cli\u003eDiBattista JD (2008) Patterns of genetic variation in anthropogenically impacted populations. Conservation Genetics 9:141\u0026ndash;156\u003c/li\u003e\n\u003cli\u003eDo C, Waples RS, Peel D, et al (2014) NeEstimator v2: re‐implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources 14:209\u0026ndash;214\u003c/li\u003e\n\u003cli\u003eDudley RK, Platania SP (2007) Flow regulation and fragmentation imperil pelagic‐spawning riverine fishes. Ecological Applications 17:2074\u0026ndash;2086\u003c/li\u003e\n\u003cli\u003eEisenhour DJ (1999) Systematics of \u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e (Cypriniformes: Cyprinidae). Copeia 969\u0026ndash;980\u003c/li\u003e\n\u003cli\u003eFranklin PA, Ba\u0026scaron;ić T, Davison PI, et al (2024) Aquatic connectivity: challenges and solutions in a changing climate. Journal of Fish Biology 105:392\u0026ndash;411\u003c/li\u003e\n\u003cli\u003eGalv\u0026aacute;n‐Femen\u0026iacute;a I, Graffelman J, Barcel\u0026oacute;‐i‐Vidal C (2017) Graphics for relatedness research. Molecular Ecology Resources 17:1271\u0026ndash;1282\u003c/li\u003e\n\u003cli\u003eGarrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. arXiv preprint arXiv:12073907\u003c/li\u003e\n\u003cli\u003eGilbert KJ, Whitlock MC (2015) Evaluating methods for estimating local effective population size with and without migration. Evolution 69:2154\u0026ndash;2166\u003c/li\u003e\n\u003cli\u003eGompert Z, Springer A, Brady M, et al (2021) Genomic time‐series data show that gene flow maintains high genetic diversity despite substantial genetic drift in a butterfly species. Molecular Ecology 30:4991\u0026ndash;5008\u003c/li\u003e\n\u003cli\u003eHill WG (1981) Estimation of effective population size from data on linkage disequilibrium1. Genetics Research 38:209\u0026ndash;216\u003c/li\u003e\n\u003cli\u003eJombart T (2008) adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403\u0026ndash;1405\u003c/li\u003e\n\u003cli\u003eJombart T, Ahmed I (2011) adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27:3070\u0026ndash;3071\u003c/li\u003e\n\u003cli\u003eJones AT, Ovenden JR, Wang Y-G (2016) Improved confidence intervals for the linkage disequilibrium method for estimating effective population size. Heredity 117:217\u0026ndash;223\u003c/li\u003e\n\u003cli\u003eJorde PE, Ryman N (1995) Temporal allele frequency change and estimation of effective size in populations with overlapping generations. Genetics 139:1077\u0026ndash;1090\u003c/li\u003e\n\u003cli\u003eKeenan K, McGinnity P, Cross TF, et al (2013) diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods in ecology and evolution 4:782\u0026ndash;788\u003c/li\u003e\n\u003cli\u003eKuczynski L, Legendre P, Grenouillet G (2018) Concomitant impacts of climate change, fragmentation and non‐native species have led to reorganization of fish communities since the 1980s. Global Ecology and Biogeography 27:213\u0026ndash;222\u003c/li\u003e\n\u003cli\u003eLangmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nature Methods 9:357\u003c/li\u003e\n\u003cli\u003eLeberg P (2002) Estimating allelic richness: effects of sample size and bottlenecks. Molecular Ecology 11:2445\u0026ndash;2449\u003c/li\u003e\n\u003cli\u003eLi H (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34:3094\u0026ndash;3100\u003c/li\u003e\n\u003cli\u003eLi H, Handsaker B, Wysoker A, et al (2009) The sequence alignment/map format and SAMtools. Bioinformatics 25:2078\u0026ndash;2079\u003c/li\u003e\n\u003cli\u003eLiaw A, Wiener M (2002) Classification and regression by randomForest. R news 2:18\u0026ndash;22\u003c/li\u003e\n\u003cli\u003eLuttrell GR, Echelle AA, Fisher WL, Eisenhour DJ (1999) Declining status of two species of the Macrhybopsis aestivalis complex (Teleostei: Cyprinidae) in the Arkansas River basin and related effects of reservoirs as barriers to dispersal. Copeia 981\u0026ndash;989\u003c/li\u003e\n\u003cli\u003eMangiafico S (2021) package rcompanion: Functions to Support Extension Education Program Evaluation. Rutgers Cooperative Extension, New Brunswick, New Jersey Version 2:30\u003c/li\u003e\n\u003cli\u003eMar\u0026ccedil;ais G, Kingsford C (2011) A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics 27:764\u0026ndash;770\u003c/li\u003e\n\u003cli\u003eMarschall EA, Crowder LB (1996) Assessing population responses to multiple anthropogenic effects: a case study with brook trout. Ecological Applications 6:152\u0026ndash;167\u003c/li\u003e\n\u003cli\u003eMathieu‐B\u0026eacute;gn\u0026eacute; E, Loot G, Chevalier M, et al (2019) Demographic and genetic collapses in spatially structured populations: Insights from a long‐term survey in wild fish metapopulations. Oikos 128:196\u0026ndash;207\u003c/li\u003e\n\u003cli\u003eMeirmans PG (2020) Genodive version 3.0: Easy‐to‐use software for the analysis of genetic data of diploids and polyploids. Molecular Ecology Resources 20:1126\u0026ndash;1131\u003c/li\u003e\n\u003cli\u003eMesser PW, Ellner SP, Hairston NG (2016) Can population genetics adapt to rapid evolution? Trends in Genetics 32:408\u0026ndash;418\u003c/li\u003e\n\u003cli\u003eMorrissey MB, de Kerckhove DT (2009) The maintenance of genetic variation due to asymmetric gene flow in dendritic metapopulations. The American Naturalist 174:875\u0026ndash;889\u003c/li\u003e\n\u003cli\u003eNei M (1987) Molecular evolutionary genetics. Columbia University Press\u003c/li\u003e\n\u003cli\u003eNei M, Tajima F (1981) Genetic drift and estimation of effective population size. Genetics 98:625\u0026ndash;640\u003c/li\u003e\n\u003cli\u003eNunn AD, Moccetti P, H\u0026auml;nfling B, et al (2024) The genome sequence of the Eurasian minnow, \u003cem\u003ePhoxinus phoxinus\u003c/em\u003e (Linnaeus, 1758). Wellcome Open Research 9:504\u003c/li\u003e\n\u003cli\u003eOsborne MJ, Archdeacon TP, Yackulic CB, et al (2023a) Genetic erosion in an endangered desert fish during a megadrought despite long‐term supportive breeding. Conservation Biology e14154\u003c/li\u003e\n\u003cli\u003eOsborne MJ, Caeiro‐Dias G, Turner TF (2023b) Transitioning from microsatellites to SNP‐based microhaplotypes in genetic monitoring programmes: Lessons from paired data spanning 20 years. Molecular Ecology 32:316\u0026ndash;334\u003c/li\u003e\n\u003cli\u003eOsborne MJ, Carson EW, Turner TF (2012) Genetic monitoring and complex population dynamics: insights from a 12‐year study of the Rio Grande silvery minnow. Evolutionary Applications 5:553\u0026ndash;574\u003c/li\u003e\n\u003cli\u003eOsborne MJ, Hatt JL, Gilbert EI, Davenport SR (2021) Still time for action: genetic conservation of imperiled South Canadian River fishes, Arkansas River Shiner (\u003cem\u003eNotropis girardi\u003c/em\u003e), Peppered Chub (\u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e) and Plains Minnow (Hybognathus placitus). Conservation Genetics 22:927\u0026ndash;945\u003c/li\u003e\n\u003cli\u003eOsborne MJ, Hudgell MAB, Caeiro-Dias G, Turner TF (2023c) The complete mitochondrial genomes of two imperiled species endemic to the Southwestern United States: Peppered Chub (\u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e) and Gila Trout (\u003cem\u003eOncorhynchus gilae\u003c/em\u003e). Mitochondrial DNA Part B 8:809\u0026ndash;814\u003c/li\u003e\n\u003cli\u003eParadis E (2010) pegas: an R package for population genetics with an integrated\u0026ndash;modular approach. Bioinformatics 26:419\u0026ndash;420\u003c/li\u003e\n\u003cli\u003ePaz-Vinas I, Loot G, Hermoso V, et al (2018) Systematic conservation planning for intraspecific genetic diversity. Proceedings of the Royal Society B: Biological Sciences 285:20172746\u003c/li\u003e\n\u003cli\u003ePennock CA, Gido KB, Perkin JS, et al (2017) Collapsing range of an endemic Great Plains minnow, Peppered Chub \u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e. The American Midland Naturalist 177:57\u0026ndash;68\u003c/li\u003e\n\u003cli\u003ePerkin JS, Gido KB (2011) Stream fragmentation thresholds for a reproductive guild of Great Plains fishes. Fisheries 36:371\u0026ndash;383\u003c/li\u003e\n\u003cli\u003ePerkin JS, Gido KB, Costigan KH, et al (2015) Fragmentation and drying ratchet down Great Plains stream fish diversity. Aquatic Conservation: Marine and Freshwater Ecosystems 25:639\u0026ndash;655\u003c/li\u003e\n\u003cli\u003ePerkin JS, Starks TA, Pennock CA, et al (2019) Extreme drought causes fish recruitment failure in a fragmented Great Plains riverscape. Ecohydrology 12:e2120\u003c/li\u003e\n\u003cli\u003ePlatania SP, Altenbach CS (1998) Reproductive strategies and egg types of seven Rio Grande basin cyprinids. Copeia 559\u0026ndash;569\u003c/li\u003e\n\u003cli\u003ePlatania SP, Mortensen JG, Farrington MA, et al (2020) Dispersal of stocked Rio Grande silvery minnow (\u003cem\u003eHybognathus amarus\u003c/em\u003e) in the middle Rio Grande, New Mexico. The Southwestern Naturalist 64:31\u0026ndash;42\u003c/li\u003e\n\u003cli\u003eRagsdale AP, Gravel S (2020) Unbiased estimation of linkage disequilibrium from unphased data. Molecular Biology and Evolution 37:923\u0026ndash;932\u003c/li\u003e\n\u003cli\u003eRanallo-Benavidez TR, Jaron KS, Schatz MC (2020) GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nature Communications 11:1432\u003c/li\u003e\n\u003cli\u003eRussello MA, Waterhouse MD, Etter PD, Johnson EA (2015) From promise to practice: pairing non-invasive sampling with genomics in conservation. PeerJ 3:e1106\u003c/li\u003e\n\u003cli\u003eRyan SF, Deines JM, Scriber JM, et al (2018) Climate-mediated hybrid zone movement revealed with genomics, museum collection, and simulation modeling. Proceedings of the National Academy of Sciences 115:E2284\u0026ndash;E2291\u003c/li\u003e\n\u003cli\u003eSteffensmeier ZD, Mayes KB, Perkin JS (2024) Linking short-term movement rate of pelagic-broadcast spawning fishes to river fragment length and conservation status. Biological Conservation 293:110585\u003c/li\u003e\n\u003cli\u003eStuart IG, Sharpe CP (2020) Riverine spawning, long distance larval drift, and floodplain recruitment of a pelagophilic fish: a case study of golden perch (\u003cem\u003eMacquaria ambigua\u003c/em\u003e) in the arid Darling River, Australia. Aquatic Conservation: Marine and Freshwater Ecosystems 30:675\u0026ndash;690\u003c/li\u003e\n\u003cli\u003eTorterotot J-B, Perrier C, Bergeron NE, Bernatchez L (2014) Influence of forest road culverts and waterfalls on the fine-scale distribution of brook trout genetic diversity in a boreal watershed. Transactions of the American Fisheries Society 143:1577\u0026ndash;1591\u003c/li\u003e\n\u003cli\u003eU.S. Fish and Wildlife Service (2018) Species status assessment report for the Arkansas River shiner (\u003cem\u003eNotropis girardi\u003c/em\u003e) and peppered chub (\u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e), version 1.0, with appendices. October 2018. Albuquerque, NM. 172 pp.\u003c/li\u003e\n\u003cli\u003eU. S. Fish and Wildlife Service (2022) Endangered and threatened wildlife and plants; endangered species status for the peppered chub and designation of critical habitat. Federal Register 87:11188\u0026ndash;11220.\u003c/li\u003e\n\u003cli\u003eWaples RS (2005) Genetic estimates of contemporary effective population size: to what time periods do the estimates apply? Molecular Ecology 14:3335\u0026ndash;3352\u003c/li\u003e\n\u003cli\u003eWaples RS, Do C (2010) Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: a largely untapped resource for applied conservation and evolution. Evolutionary Applications 3:244\u0026ndash;262\u003c/li\u003e\n\u003cli\u003eWilde GR, Durham BW (2008) A life history model for peppered chub, a broadcast-spawning cyprinid. Transactions of the American Fisheries Society 137:1657\u0026ndash;1666\u003c/li\u003e\n\u003cli\u003eWillis SC, Hollenbeck CM, Puritz JB, et al (2017) Haplotyping RAD loci: an efficient method to filter paralogs and account for physical linkage. Molecular Ecology Resources 17:955\u0026ndash;965\u003c/li\u003e\n\u003cli\u003eYackulic CB, Archdeacon TP, Valdez RA, et al (2022) Quantifying flow and nonflow management impacts on an endangered fish by integrating data, research, and expert opinion. Ecosphere 13:e4240\u003c/li\u003e\n\u003cli\u003eZheng X, Levine D, Shen J, et al (2012) A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28:3326\u0026ndash;3328\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"genomic time-series, neutral genetic diversity, peppered chub, Macrhybopsis tetranema, South Canadian River","lastPublishedDoi":"10.21203/rs.3.rs-7863499/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7863499/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003ePatterns of genetic variation within species may differ in space and time, particularly in the face of rapid anthropogenic changes. Genomic time-series population-level data can reveal patterns of allele frequency change and can help understanding on-going evolutionary processes. This study used a time-series spanning 6 years to evaluate genomic variability and trends in effective population size, and to assess the influence of genetic drift on \u003cem\u003eMacrhybopsis tetranema\u003c/em\u003e genomic diversity. Temporal trends were evaluated from samples collected from 2015 to 2020. Results from genomic variation and variance effective population size across time suggests \u003cem\u003eM. tetranema\u003c/em\u003e upstream (New Mexico) are affected by genetic drift, but genetic diversity did not change significantly across the time-series. Absolute allele frequency differences across time-series revealed that drift is a result of many loci with moderate allele frequency changes, particularly after 2018, rather than few loci with high frequency changes. Also, linkage disequilibrium effective population size showed fluctuations and a general increasing trend. Overall, the results presented here demonstrated that genetic drift is an overarching evolutionary force driving trends in genetic variation in \u003cem\u003eM. tetranema\u003c/em\u003e. Together, genetic data reported here and ecological data from recent studies, suggest that \u003cem\u003eM. tetranema\u003c/em\u003e in the upstream portion of the remnant distribution experiences genetic drift likely due of downstream-biased gene flow, while upstream movement of adults can explain low relatedness and maintenance of genetic diversity upstream in the face of genetic drift.\u003c/p\u003e","manuscriptTitle":"Sustaining genetic diversity in an imperiled pelagophilic fish despite genetic drift","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-10-22 09:51:00","doi":"10.21203/rs.3.rs-7863499/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"1d9ca7e3-3d7a-408a-b8e9-5037ea2ac6bc","owner":[],"postedDate":"October 22nd, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2025-11-25T07:24:04+00:00","versionOfRecord":[],"versionCreatedAt":"2025-10-22 09:51:00","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7863499","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7863499","identity":"rs-7863499","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}