Accepted Article

Received Date : 05-Sep-2014 Revised Date : 01-Mar-2015 Accepted Date : 04-Mar-2015 Article type : Original Article

Past climate change drives current genetic structure of an endangered freshwater mussel species Kentaro Inoue (KI)1, Brian K. Lang (BKL)2, David J. Berg (DJB)3 1

Department of Biology, Miami University, 700 E. High Street, Oxford, Ohio 45056 23 Llanito Road, Bernalillo, New Mexico 87004 3 Department of Biology, Miami University, 1601 University Boulevard, Hamilton, OH 45011 2

Keywords: approximate Bayesian computation (ABC), aquatic connectivity, ecological niche modeling, Chihuahuan Desert, desert streams Corresponding author: K. Inoue, Department of Biology, Miami University, 700 E. High Street, Oxford, Ohio 45056, USA fax: 513-529-6900, email: [email protected] Running title: Climate change alters aquatic connectivity

Abstract Historical-to-recent climate change and anthropogenic disturbance affect species distributions and genetic structure. The Rio Grande watershed of the United States and Mexico encompasses ecosystems that are intensively exploited, resulting in substantial degradation of aquatic habitats. While significant anthropogenic disturbances in the Rio Grande are recent, inhospitable conditions for freshwater organisms likely existed prior to such disturbances. A combination of anthropogenic and past climate factors may contribute to current distributions of aquatic fauna in the Rio Grande basin. We used mitochondrial DNA and 18 microsatellite loci to infer evolutionary history and genetic structure of an endangered freshwater mussel, Popenaias

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/mec.13156 This article is protected by copyright. All rights reserved.

Accepted Article

popeii, throughout the Rio Grande drainage. We estimated spatial connectivity and gene flow across extant populations of P. popeii and used ecological niche models (ENMs) and approximate Bayesian computation (ABC) to infer its evolutionary history during the Pleistocene.

STRUCTURE

results recovered regional and local population clusters in the Rio

Grande. ENMs predicted drastic reductions in suitable habitat during the last glacial maximum. ABC analyses suggested that regional population structure likely arose in this species during the mid-to-late-Pleistocene, and was followed by a late-Pleistocene population bottleneck in New Mexico populations. The local population structure arose relatively recently, perhaps due to anthropogenic factors. Popenaias popeii, one of the few freshwater mussel species native to the Rio Grande basin, is a case study for understanding how both geological and anthropogenic factors shape current population genetic structure. Conservation strategies for this species should account for the fragmented nature of contemporary populations.

Introduction A major interest in evolution and population genetics is the impact of past climate

conditions and vicariant events on current species’ distributions and the genetic structure of natural populations (Avise 2000; Hewitt 2000). Climate change and geological processes often cause strong directional selection in natural populations (Hewitt 2000). For example, during the Pleistocene interglacial periods, times of drastic climate change altered freshwater systems by rerouting rivers and inundating large areas of land with meltwater (Hocutt et al. 1978; Mayden

1988; Strange & Burr 1997). Consequently in the Northern Hemisphere, organisms were extirpated over northern portions of their ranges, survived in isolated populations, and/or dispersed to new locations (Hewitt 2000). In fragmented populations, neutral forces such as

This article is protected by copyright. All rights reserved.

Accepted Article

genetic drift may have shaped genetic structure. Thus, current species distributions and genetic structure at both intra- and interspecific levels were formed through isolation of populations via vicariant events followed by recolonization with secondary contact during Pleistocene glacialinterglacial cycles (Mäkinen & Merilä 2008; Bossu et al. 2013; Inoue et al. 2014b). Studies of

evolutionary responses to vicariant events are especially timely, given rapid changes in ecological conditions due to climate change and anthropogenic disturbance.

Recently, the development of advanced statistical methods has led to improved ability to

reconstruct the ways in which past climate change and geological events have shaped existing species’ distributions and phylogeographic patterns (Knowles 2009; Fordham et al. 2014).

Coalescent-based population genetic approaches allow researchers to infer the evolutionary histories of species. For example, with the development of coalescent methods via approximate Bayesian computation (ABC), we can test hypothesized demographic and evolutionary scenarios by creating many simulated datasets, measuring similarity between those datasets and empirical data, and estimating the posterior probability that a given scenario is a true model (Beaumont 2010). By using approximations of likelihood inference, ABC can manage much more complex

demographic and evolutionary models when full-likelihood methods are not efficient (Nielsen & Beaumont 2009). Ecological niche models (ENMs), which predict suitable habitats for a given species by estimating the relationship between species occurrences and environmental characteristics of locations, also provide powerful tools for evolutionary biologists, biogeographers, paleoecologists, and conservation biologists (Elith et al. 2006; Fordham et al. 2014). The combination of coalescent-based approaches and ENMs promises to shed light on a variety of issues, including the evolutionary histories of species, effects of population

This article is protected by copyright. All rights reserved.

Accepted Article

the relationships into pairs of localities within and among rivers, only DEST showed slight positive isolation-by-distance patterns within rivers (Fig. S1E).

Population genetic structure Popenaias popeii showed evidence of significant range-wide population genetic

structure. The STRUCTURE analysis recovered explicit boundaries between the Black River and Devils River/Rio Grande at k = 2, indicated by lnP(k) = −12576.0 and ∆k = 1305.5 (Figs 4 and

S2, Supporting Information). At k = 3, STRUCTURE further split the BR sites into two clusters, and this assignment scheme received the highest log-likelihood estimate, lnP(k) = −12473.7, and the second best ∆k estimate, ∆k = 63.4 (Figs 4 and S2). The BR sites were split between three upstream sites (hereafter BR-u) and five downstream sites (hereafter BR-d) at k = 3, where BR-d

had admixture of BR-u genotypes (Fig. 4). We found no evidence of admixture between BR and RG populations when k = 2 or k = 3. Generally, ∆k is thought to capture the strongest patterns of

population structure (Evanno et al. 2005; Coulon et al. 2008), and the highest lnP(k) can be useful to investigate finer population structure (e.g., Coulon et al. 2008; Fisher-Reid et al. 2013;

but see Kalinowski 2011). The STRUCTURE analysis on the BR sites alone showed the split between BR-u and BR-d at k = 2 received the highest log-likelihood and ∆k estimates (lnP(k) = −6965.0, ∆k = 5.0); this assignment scheme was consistent with the overall dataset. Thus, our

STRUCTURE

results showed that both the best likelihood and the second best ∆k reflect fine-scale

structure in the Black River. Because the BR-u and BR-d groups of sites are separated by a lowwater crossing which may inhibit fish movement (see Discussion), we concluded that k = 3 represents the most biologically relevant clustering scheme for P. popeii.

This article is protected by copyright. All rights reserved.

Accepted Article

Berg Lab, and the Center for Bioinformatics and Functional Genomics at Miami University. We thank A. Walters for help with ecological niche modeling and L. Burlakova for providing us with collection records for P. popeii. We also thank private landowners for access to the river.

Comments from C. Franz Berg and four anonymous reviewers greatly improved the manuscript. Funding was provided by the New Mexico Department of Game and Fish, and the US Fish and Wildlife Service.

References Alvarado-Serrano DF, Knowles LL (2014) Ecological niche models in phylogeographic studies: applications, advances and precautions. Molecular Ecology Resources, 14, 233-248. Antao T, Lopes A, Lopes RJ, Beja-Pereira A, Luikart G (2008) LOSITAN: a workbench to detect molecular adaptation based on a Fst-outlier method. BMC Bioinformatics, 9, 323. Avise JC (2000) Phylogeography: The History and Formation of Species Harvard University Press, Cambridge, Massachusetts. Balloux F, Lehmann L (2003) Random mating with a finite number of matings. Genetics, 165, 2313-2315. Beaumont MA (2010) Approximate Bayesian computation in evolution and ecology. Annual Review of Ecology, Evolution, and Systematics, 41, 379-406. Beerli P (2006) Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics, 22, 341-345. Bequaert JC, Miller WB (1973) The Mollusks of the Arid Southwest: With an Arizona Check List University of Arizona Press, Tucson, Arizona. Bogan AE (1996) Popenaias popeii. In: IUCN 2013. IUCN Red List of Threatened Species. Bond N, Thomson J, Reich P, Stein J (2011) Using species distribution models to infer potential climate change-induced range shifts of freshwater fish in south-eastern Australia. Marine and Freshwater Research, 62, 1043-1061. Bossu CM, Beaulieu JM, Ceas PA, Near TJ (2013) Explicit tests of palaeodrainage connections of southeastern North America and the historical biogeography of Orangethroat Darters (Percidae: Etheostoma: Ceasia). Molecular Ecology, 22, 5397-5417. Brown JL (2014) SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods in Ecology and Evolution, 5, 694-700. Carman SM (2007) Texas hornshell Popenaias popeii recovery plan. New Mexico Department of Game and Fish, Conservation Services Division, Santa Fe, New Mexico. Carstens BC, Richards CL (2007) Integrating coalescent and ecological niche modeling in comparative phylogeography. Evolution, 61, 1439-1454. Chagnon G, de la Cheneliere V (1998) Status assessment and conservation action plan for Lampsilis radiata. McGill University, Montreal, Quebec, Canada. Christian AD, Monroe EM, Asher AM, M. LJ, Berg DJ (2007) Methods of DNA extraction and PCR amplification for individual freshwater mussel (Bivalvia: Unionidae) glochidia, with

This article is protected by copyright. All rights reserved.

Accepted Article

the first report of multiple paternity in these organisms. Molecular Ecology Resources, 7, 570–573. Clement M, Posada D, Crandall KA (2000) TCS: a computer program to estimate gene genealogies. Molecular Ecology, 9, 1657-1659. Contreras-Balderas S, Lozano-Vilano ML (1994) Water, endangered fishes, and development perspectives in arid lands of Mexico. Conservation Biology, 8, 379-387. Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R, Marin JM, Estoup A (2014) DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics, 30, 1187-1189. Cornuet JM, Ravigne V, Estoup A (2010) Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0). BMC Bioinformatics, 11, 401. Coulon A, Fitzpatrick JM, Bowman R, Stith BM, Makarewich CA, Stenzler LM, Lovette IJ (2008) Congruent population structure inferred from dispersal behaviour and intensive genetic surveys of the threathened Florida scrub-jay (Aphelocoma coerulenscens). Molecular Ecology, 17, 1685-1701. Davis JR (1980) Species composition and diversity of benthic macroinvertebrate population of the Pecos River, Texas. Southwestern Naturalist, 25, 241-256. Earl DA, vonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources, 4, 359-361. Edwards RJ, Garrett GP, Marsh-Matthews E (2002) Conservation and status of the fish communities inhabiting the Rio Conchos basin and middle Rio Grande, Mexico and U.S.A. Reviews in Fish Biology and Fisheries, 12, 119-132. Elith J, Graham CH, Anderson RP, Dudik M, Ferrier S, Guisan A, Hijmans RJ, Huettmann F, Leathwick JR, Lehmann A, Li J, Lohmann LG, Loiselle BA, Manion G, Moritz C, Nakamura M, Nakazawa Y, Overton JM, Peterson AT, Phillips SJ, Richardson K, Scachetti-Pereira R, Schapire RE, Soberon J, Williams S, Wisz MS, Zimmermann NE (2006) Novel methods improve prediction of species' distributions from occurrence data. Ecography, 29, 129-151. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611-2620. Fetzner JW, Jr., Crandall KA (2003) Linear habitats and the nested clade analysis: an empirical evaluation of geographic versus river distances using an Ozark crayfish (Decapoda: Cambaridae). Evolution, 57, 2101-2118. Fisher-Reid MC, Engstrom TN, Kuczynski CA, Stephens PR, Wiens JJ (2013) Parapatric divergence of sympatric morphs in a salamander: incipient speciation on Long Island? Molecular Ecology, 22, 4681-4694. Fordham DA, Brook BW, Moritz C, Nogues-Bravo D (2014) Better forecasts of range dynamics using genetic data. Trends in Ecology & Evolution, 29, 436-443. Galloway WE, Whiteaker TL, Ganey-Curry P (2011) History of Cenozoic North American drainage basin evolution, sediment yield, and accumulation in the Gulf of Mexico basin. Geosphere, 7, 938-973. Goudet J (1995) FSTAT (version 1.2): a computer program to calculate F-statistics. Journal of Heredity, 86, 485-486.

This article is protected by copyright. All rights reserved.

Accepted Article

Materials and Methods Sampling, genetic data collection, and analyses of population genetic diversity Using snorkeling, SCUBA, and tactile methods, we sampled 254 P. popeii from two

tributaries and one main-stem region of the Rio Grande basin, including 193 individuals from

eight locations in the Black River, New Mexico; 58 individuals from five locations in the mainstem Rio Grande, Texas; and three individuals from three locations in the Devils River, Texas (Table 1; Fig. 1). Given the small sample size, we pooled samples from the Devils River into a single “population” prior to analyses. We sampled nondestructively using tissue swabs, collected from between foot and mantle tissues by rubbing mucus and epidermal cells using sample collection swabs (Epicentre Biotechnologies, Madison, WI), and then returned mussels to the riverbed. Samples were preserved in 95% ethanol and stored at −20˚C. Total genomic DNA was extracted using ArchivePure DNA Cell/Tissue Kits (5 Prime, Gaithersburg, MD), diluted to 10 ng/µL, and used as a template in polymerase chain reactions (PCR) for mtDNA and microsatellite analyses.

We used PRIMER3 (Untergasser et al. 2012) to design primers (forward: 5′-

TGTGGGGTGAATCATTCCTT-3′ and reverse: 5′-TAAACCTCAGGATGCCCAAA-3′) from a

complete mtDNA genome of the mussel Lampsilis ornata (GenBank accession number: NC_005335). These primers amplified about 810 basepairs of part of the cytochrome oxidase II gene and the cytochrome oxidase I gene (hereafter collectively abbreviated as COX). We followed procedures for conditions for PCR, sequencing, and post-sequencing analyses described in Inoue et al. (2014b).

This article is protected by copyright. All rights reserved.

Accepted Article

Using DNASP v5.10 (Librado & Rozas 2009), we estimated population genetic indices from mtDNA sequences, including number of haplotypes (H), mean number of basepair differences (K), and mean nucleotide diversity (π) over the pooled dataset and within each locality and river because we did not have a priori population genetic information. To correct for sample-size bias, we estimated rarefied number of haplotypes (HR) with a standardized

sample size of six individuals using the VEGAN package (Oksanen et al. 2015) in R v3.0.3 (R Development Core Team 2011). We then used TCS v1.21 (Clement et al. 2000) to build a 95%

confidence parsimony network from COX haplotypes. Assigning the shortest path from the most frequent haplotype simplified multiple connections between haplotypes (Fetzner & Crandall 2003).

We genotyped populations of P. popeii at 20 tetra-nucleotide microsatellite loci (Inoue et

al. 2013). Forward primers for each PCR were labeled with a 5′ fluorescent tag (6-FAM, NED, PET, or VIC) for visualization. We performed five sets of multiplex PCR (Plex1: Tetra17-1941; Plex2: Tetra01-09-14-24-30-36; Plex3: Tetra02-03-22-23; Plex4: Tetra05-31-40; and Plex5:

Tetra08-15-33-37) designed by MULTIPLEXMANAGER (Holleley & Geerts 2009). Thermal cycling began with initial denaturing at 95˚C for 2 min; 35 cycles of 94˚C for 30 s, 60˚C for 1.5 min, and 72˚C for 1 min; and final extension at 72˚C for 30 min. We used previously published procedures for fragment analyses, allele scoring, and assignment of integer numbers to DNA fragment sizes (Inoue et al. 2014b).

This article is protected by copyright. All rights reserved.

Accepted Article

We tested for the presence of null alleles and large allele dropout using MICRO-CHECKER (van Oosterhout et al. 2004). We checked for microsatellite loci under directional or balancing selection using LOSITAN (Antao et al. 2008). Using GENEPOP v4.0.10 (Rousset 2008), we

conducted exact tests for pairwise linkage disequilibrium and deviation from Hardy-Weinberg expectation (HWE) for each locality, along with several population genetic indices (mean number of alleles per locus, NA; observed and expected heterozygosities, HO and HE; and number

of private alleles, NP) using GENALEX v6.3 (Peakall & Smouse 2006). We used rarefaction to a standardized sample size of six individuals to correct mean allelic richness (rarefied number of alleles per locus; AR) in FSTAT v2.9.3 (Goudet 1995).

To assess the effect of geographic distance on genetic structure, we examined the

correlation of matrices of pairwise genetic differences (φST for mtDNA sequences and DEST for

microsatellites) and geographic distances between pairs of localities using Mantel tests. We estimated the pairwise genetic differences using GENALEX. We measured total distance between pairs of localities along the rivers (river distance) using ArcGIS v10.2 (ESRI, Inc.). We excluded localities with small population size (0.75; Table S1), indicating overall adequate model

performance. In the ENMs for P. popeii, maximum temperature of warmest month (41.0%) and mean diurnal temperature range (24.2%) were the best predictors of suitable environments,

This article is protected by copyright. All rights reserved.

Accepted Article

followed by precipitation seasonality (15.4%) and precipitation of warmest quarter (13.9%; Table S1). The present-day ENM model predicted suitable environment that closely matched the known distribution of P. popeii in the USA (Fig. 5). Additionally, the ENMs predicted potential suitable environment in the lower Pecos River and the Rio Salado. ENMs of the LGM model predicted drastic reduction of suitable environment with low suitability scores. However, the LIG predicted suitable environment similar to the present-day model, but with an even greater total area of predicted suitable environment (Fig. 5). These results suggest that P. popeii populations were potentially continuously distributed and much more widespread throughout the Rio Grande drainage during the LIG, but the range was reduced and populations fragmented during the LGM.

The ENMs for the host species predicted wider distributional ranges of suitable

environments than those of P. popeii (Fig. S3, Supporting Information). Like ENMs for P. popeii, maximum temperature of warmest month best predicted suitable environments for the host species; however, additional environmental variables varied among host species (Table S1). Consistent with the ENM results for P. popeii, the ENM models of the paleodistributions of the host species predicted a wide range of suitable environments in the Rio Grande drainage during the LIG, followed by drastic reduction of suitable environments during the LGM.

Approximate Bayesian computation analysis of demographic history Our ABC analysis using DIYABC identified Scenario 3, which specified a population

bottleneck in the ancestral Black River, as the most highly supported scenario, with strong posterior probability (0.740; Table 5). All other scenarios received much lower support. Based

This article is protected by copyright. All rights reserved.

Accepted Article

on our ABC analyses of pseudo-observed datasets, the Type I error rate for this most-probable scenario was relatively high (Table 5), but the Type II error rate was low. Although a high Type I error indicated that this scenario was less frequently chosen when it was the true scenario, the high posterior probability and low Type II error suggest greater confidence that this is the best of the proposed scenarios.

The RMAEs were relatively high for NBR and te2, but low for the rest of the parameters,

indicating most parameters estimated by ABC were reliable values (Table 2). The high posterior probability of Scenario 3 indicated that the divergence from MRCAP occurred approximately 9020 generations ago (95% CI: 2370 – 18,900; RMAE = 0.185). This was followed by a severe population bottleneck event in the ancestral BR population 3670 generations ago (95% CI: 614 – 11,000; RMAE = 0.316), and divergence between BR-u and BR-d about 125 generations ago (95% CI: 12 – 410; RMAE = 0.298). Although no species-specific generation time for P. popeii is known, when using an average generation time of 8.9 years (based on 8.1 - 9.6 years for Lampsilis radiata; Chagnon & de la Cheneliere 1998), we inferred that the time of the divergence from MRCAP, the population bottleneck, and the divergence between BR-u and BRd are 80.3 ka (21.1 – 168.2 ka), 32.7 ka (5.5 – 97.9 ka), and 1.1 ka (0.1 – 3.6 ka), respectively. During the bottleneck event, the ancestral BR population fell by 96% (from Ne = 39,400 to 1630 individuals; Table 2). Posterior distributions of current Ne varied among populations; the largest

Ne was RG with 70,300 individuals (41,600 – 91,900; RMAE = 0.200), whereas the two BR populations were small (1190 individuals for BR-u and 4680 individuals for BR-d). Mean mutation rate over 18 microsatellite loci was estimated to be 1.53 × 10-4 per generation (7.48 ×

10-5 – 3.39 × 10-4, RMAE = 0.214).

This article is protected by copyright. All rights reserved.

Accepted Article

Discussion Spatial and demographic history of Popenaias popeii during the Pleistocene We found that P. popeii populations were highly structured among Rio Grande drainages

at a regional scale, and within the Black River at a more local scale. Based on our population genetic modeling and ENM analyses, we conclude that this pattern probably derives from historical demographic events associated with climate change. For example, the ENM results indicated considerable suitable environment for P. popeii during the LIG (ca. 120 – 140 ka; Fig. 5), suggesting spatial connectivity among populations, followed by drastic reductions in suitable environment in the LGM (ca. 12 ka; Fig. 5). These distributional shifts are nearly congruent with time estimates from coalescent simulations; ABC analyses suggested that population divergence from the MRCAP occurred during the mid-to-late Pleistocene (ca. 80 ka; Table 2) followed by a severe population bottleneck in the ancestral BR population (ca. 33 ka; Table 2), but not in the RG population. Even though P. popeii might not have had continuous suitable environments throughout the Pecos River during the LIG, the continuous distribution of suitable environments for all host species indicates that P. popeii was capable of dispersing throughout the drainage. In fact, estimates of historic gene flow were high but rather unidirectional from RG populations to BR populations, suggesting that RG populations were source populations for the Black River. Although MIGRATE-N and ABC analyses have different assumptions, where the former assumes populations in migration-drift equilibrium and the latter assumes divergence without gene flow, both analyses indicate that RG and BR populations diverged during the midto-late-Pleistocene and RG populations retained most of the ancestral polymorphisms. Furthermore, BR individuals share a single COX haplotype—the most common haplotype in the RG population—also consistent with historical migration between the rivers and recent

This article is protected by copyright. All rights reserved.

Accepted Article

divergence. Apparently, subsequent glaciation reduced the extent of suitable environment during the LGM (ca. 21 ka); such environments for P. popeii and its hosts were small and fragmented, presumably leading to greatly decreased levels of dispersal and gene flow among populations and a substantial bottleneck in the ancestral BR population.

Although glaciers advanced only to southern New Mexico’s Sierra Blanca (Smith &

Miller 1986), continental glaciation caused cooler summers and a significant increase in winter precipitation in much of the Chihuahuan Desert during the Pleistocene glacial-interglacial periods (Metcalfe et al. 2000). Paleoecological studies have shown that the Pleistocene’s dynamic climate induced distributional shifts in unglaciated areas. For example, glacialinterglacial cycles prompted severe range contractions and population reductions in currently widespread species of fishes (Miller 1977) and more geographically restricted populations of land snails (Bequaert & Miller 1973; Metcalf 1997) in the Chihuahuan Desert. More recent studies using genetic and ENM approaches identified population isolation followed by species diversification in small mammals in North American deserts (Jezkova et al. 2009; Mantooth et al. 2013). Similarly, our results indicate that the LGM’s cooler climate reduced suitable environment for P. popeii, likely isolating the ancestral BR and RG populations. The subsequent

bottleneck substantially reduced Ne of the ancestral BR population; a lack of connectivity with

the Rio Grande likely caused current relatively low levels of genetic variation and led to complete lineage sorting of COX haplotypes in the BR populations. In contrast, we did not detect bottlenecks in RG populations from ABC analyses, indicating that RG populations maintained higher Ne, and thus, higher genetic diversity through to the present. We estimated the

timing of historical demographic events based on generalized mutation rates and average

This article is protected by copyright. All rights reserved.

Accepted Article

generation time for L. radiata, which might differ from species-specific values for P. popeii. Nevertheless, our results indicate that the BR and RG populations diverged during the mid-tolate-Pleistocene and the BR population experienced a severe population bottleneck during the late Pleistocene.

Vicariant events associated with tectonic activity and climate change can also reroute

rivers, shifting distributions of aquatic organisms. Such events frequently occurred in Rio Grande drainages during the Pleistocene (Thomas 1972; Galloway et al. 2011). Because we

used current Rio Grande drainages to project historic ENMs, the ENMs might overestimate suitable environment for P. popeii and its hosts. For example, pluvial lakes formed along the mid-Rio Grande during the Pleistocene (Metcalfe et al. 2000; Galloway et al. 2011).

Furthermore, the upper Rio Grande (above El Paso, Texas) was not connected with the lower Rio Grande. Instead, its headwaters, along with the upper Pecos River, were tributaries of the Canadian River (Mississippi River drainage) and headwaters of the Brazos River (Gulf of Mexico drainage) during the mid- and late-Pleistocene (Thomas 1972; Smith & Miller 1986). Thus, although the upper Pecos River and Rio Grande provided suitable environments for P. popeii and its host species, they likely did not occur there due to lack of stream connections. Additionally, our ENMs were based only on bioclimatic data. Although bioclimatic ENMs have

been used for freshwater fishes to predict their contemporary ranges and range shifts under future climate change scenarios (e.g., Bond et al. 2011), given the scale of resolution of such data, our bioclimatic ENMs might only give a rough indication of which parts of drainages are suitable for these species. Previous studies have found that P. popeii inhabits distinctive microhabitats, including undercut riverbanks and the bases of large boulders (Lang 2001; Karatayev et al. 2012;

This article is protected by copyright. All rights reserved.

Accepted Article

Inoue et al. 2014a), that are not distributed uniformly throughout these drainages (Inoue et al. 2014a). Additional geological and landscape information may increase accuracy of the ENMs. Currently, such information during the Pleistocene is not available for this region.

Contemporary gene flow and genetic structure Although the BR and RG populations diverged during the Pleistocene, we detected finer

spatial-scale genetic structure within segments of the Black River, perhaps driven by anthropogenic factors (Fig. 4). ABC analyses indicated recent divergence between BR-u and BR-d (0.1 – 3.6 ka; Table 2); however, because DIYABC does not allow migration among populations in demographic scenarios, the divergence time between the populations is likely overestimated. In fact, STRUCTURE and BAYESASS analyses indicated that these populations are not completely isolated.

STRUCTURE

analyses showed admixture of two clusters in BR-d (Fig.

4), whereas estimates of contemporary gene flow indicated active but unidirectional (upstream) gene flow between BR populations (Table 4). A possible explanation for the inconsistent results between STRUCTURE and BAYESASS analyses is that BR-u was recently colonized by a few individuals from BR-d, with all individuals in BR-u being comprised of partial BR-d genetic characters. In fact, population genetic indices showed that the BR-u individuals had consistently lower genetic diversity and Ne than the BR-d individuals (Tables 1 and 2) and genetic divergence

between the populations was statistically significant, but low (Table 3). The climatic history of the late Holocene (last 4000 years) in southeast New Mexico shows sequences of drier and wetter climates (Polyak & Asmerom 2001). While climate fluctuation might have caused BR populations to differentiate, it is unclear how this occurs with such small-spatial-scale isolation. Alternatively, because estimates from BAYESASS reflect migration rates over the last few

This article is protected by copyright. All rights reserved.

Accepted Article

generations and there is genetic divergence between BR-u and BR-d, a few individual of P. popeii might have recently colonized BR-u and are currently isolated by a low-water culvert crossing (built 1932-1936; Fig. 1), which can prevent fish dispersal (Warren & Pardew 1998). Future studies of host fish movements may identify the means by which fine-scale population structure is maintained in BR populations of P. popeii.

Because the lower Pecos River connects BR and RG populations, gene flow could occur

in contemporary time. However, genetic analyses and ENMs under current bioclimatic conditions support discontinuity between these populations, most likely because impoundments and occasional intermittent conditions prevent dispersal of host fishes. Red Bluff and Amistad

reservoirs (impounded in 1936 and 1969, respectively; Fig. 1) are major impoundments between BR and RG populations, and a segment of river with intermittent flow exists downstream of Red Bluff Reservoir. Historic surveys in the lower Pecos River periodically reported the presence of a single or a few individuals of P. popeii; however, the most recent survey did not (Karatayev et al. 2012). Intermittent conditions are common in many of these drainages due to seasonal drought and use of river water for irrigation; these alter both water quality and quantity (Hubbs et al. 1977; Davis 1980). Studies of benthic macroinvertebrates in the lower Pecos River showed that species richness was significantly lower than in upper reaches due to high salinity and extreme physiochemical fluctuations (Davis 1980). These conditions continue until below Sheffield, Texas, where numerous spring-fed creeks contribute groundwater that dilutes releases from naturally saline Red Bluff Reservoir (Davis 1980).

This article is protected by copyright. All rights reserved.

Accepted Article

Conservation implications Population genetic analyses and ENMs indicate that the BR and RG populations are

currently disconnected. Because the BR populations have been isolated from RG populations since the mid-Pleistocene and contain low genetic diversity, the BR and RG populations likely comprise separate evolutionarily significant units that must be considered when developing conservation strategies and species recovery plans. Furthermore, although current ENMs for P. popeii indicate stable suitable environments in the lower Rio Grande, landscape features (such as topographic complexity) and anthropogenic land- and water-uses can influence potential dispersal and availability of high-quality habitat. Over the past several decades, at least 30 springs have gone dry in the Rio Grande drainages (Contreras-Balderas & Lozano-Vilano 1994). Rising human demand for water will likely exacerbate intermittent conditions. A long-term demographic study revealed that reduced river discharge is associated with significantly decreased survival of P. popeii (Inoue et al. 2014a). Significant changes in hydrological regimes will likely cause significant declines in P. popeii populations. If the region’s human population continues to grow, intensifying anthropogenic threats, P. popeii has a low probability of persistence even in the suitable environments predicted from ENMs.

Given a previous study that estimated 48,006 individuals inhabiting the Black River

(Inoue et al. 2014a), we can estimate the ratio of genetically effective population size (Ne) and census population size (Nc) in the BR populations. With an estimated Ne of 5870 individuals in the Black River from ABC analyses, the Ne/Nc ratio is 0.12. Such estimates of the Ne/Nc ratio of

a target population are critical to understanding whether a population can persist and maintain adaptive genetic variance. In general, the Ne/Nc ratio can be as low as 10-4 to 10-3 and fluctuate

This article is protected by copyright. All rights reserved.

Accepted Article

enormously throughout time in fish and invertebrates that have very high fecundity and juvenile mortality (i.e., Type III survivorship curve; Luikart et al. 2010; Hare et al. 2011). Although freshwater mussels fit these criteria, promiscuous reproduction through factors such as multiple paternity (Christian et al. 2007) may maintain higher Ne/Nc ratios (Balloux & Lehmann 2003;

Snook et al. 2009). The Ne/Nc ratio can help unravel the relative importance of environmental,

ecological, and genetic factors that drive population persistence (Palstra & Fraser 2012).

Conclusion The Chihuahuan Desert is one of the world’s most biologically rich and diverse deserts

(Olson & Dinerstein 1998); its ecosystems harbor high numbers of endemics that are often evolutionarily unique (Metcalf & Smartt 1997; Kodric-Brown & Brown 2007). Unfortunately, desert ecosystems are fragile and recover slowly from perturbations, so they are susceptible to human activities. Our study integrated evolutionary history, estimates of distribution of suitable environment over time, and genetic connectivity among extant populations, to elucidate factors affecting current population genetic structure in a desert river. We found that climate change in the Pleistocene strongly influenced current genetic structure of P. popeii in the Rio Grande

drainage. Future changes in climate and habitat due to anthropogenic activities mean that P. popeii populations will continue to be isolated from one another. Resource managers must consider environmental and genetic features when developing species-recovery strategies. Acknowledgements We thank L. Burlakova, A. Karatayev, T. Miller, Y. Zhang, and T. Nobles for help with sample collections. Field and laboratory assistance was provided by the New Mexico Department of Game and Fish, the US Fish and Wildlife Service, numerous undergraduate technicians in the

This article is protected by copyright. All rights reserved.

Accepted Article

Berg Lab, and the Center for Bioinformatics and Functional Genomics at Miami University. We thank A. Walters for help with ecological niche modeling and L. Burlakova for providing us with collection records for P. popeii. We also thank private landowners for access to the river.

Comments from C. Franz Berg and four anonymous reviewers greatly improved the manuscript. Funding was provided by the New Mexico Department of Game and Fish, and the US Fish and Wildlife Service.

References Alvarado-Serrano DF, Knowles LL (2014) Ecological niche models in phylogeographic studies: applications, advances and precautions. Molecular Ecology Resources, 14, 233-248. Antao T, Lopes A, Lopes RJ, Beja-Pereira A, Luikart G (2008) LOSITAN: a workbench to detect molecular adaptation based on a Fst-outlier method. BMC Bioinformatics, 9, 323. Avise JC (2000) Phylogeography: The History and Formation of Species Harvard University Press, Cambridge, Massachusetts. Balloux F, Lehmann L (2003) Random mating with a finite number of matings. Genetics, 165, 2313-2315. Beaumont MA (2010) Approximate Bayesian computation in evolution and ecology. Annual Review of Ecology, Evolution, and Systematics, 41, 379-406. Beerli P (2006) Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics, 22, 341-345. Bequaert JC, Miller WB (1973) The Mollusks of the Arid Southwest: With an Arizona Check List University of Arizona Press, Tucson, Arizona. Bogan AE (1996) Popenaias popeii. In: IUCN 2013. IUCN Red List of Threatened Species. Bond N, Thomson J, Reich P, Stein J (2011) Using species distribution models to infer potential climate change-induced range shifts of freshwater fish in south-eastern Australia. Marine and Freshwater Research, 62, 1043-1061. Bossu CM, Beaulieu JM, Ceas PA, Near TJ (2013) Explicit tests of palaeodrainage connections of southeastern North America and the historical biogeography of Orangethroat Darters (Percidae: Etheostoma: Ceasia). Molecular Ecology, 22, 5397-5417. Brown JL (2014) SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods in Ecology and Evolution, 5, 694-700. Carman SM (2007) Texas hornshell Popenaias popeii recovery plan. New Mexico Department of Game and Fish, Conservation Services Division, Santa Fe, New Mexico. Carstens BC, Richards CL (2007) Integrating coalescent and ecological niche modeling in comparative phylogeography. Evolution, 61, 1439-1454. Chagnon G, de la Cheneliere V (1998) Status assessment and conservation action plan for Lampsilis radiata. McGill University, Montreal, Quebec, Canada. Christian AD, Monroe EM, Asher AM, M. LJ, Berg DJ (2007) Methods of DNA extraction and PCR amplification for individual freshwater mussel (Bivalvia: Unionidae) glochidia, with

This article is protected by copyright. All rights reserved.

Accepted Article

the first report of multiple paternity in these organisms. Molecular Ecology Resources, 7, 570–573. Clement M, Posada D, Crandall KA (2000) TCS: a computer program to estimate gene genealogies. Molecular Ecology, 9, 1657-1659. Contreras-Balderas S, Lozano-Vilano ML (1994) Water, endangered fishes, and development perspectives in arid lands of Mexico. Conservation Biology, 8, 379-387. Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R, Marin JM, Estoup A (2014) DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics, 30, 1187-1189. Cornuet JM, Ravigne V, Estoup A (2010) Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0). BMC Bioinformatics, 11, 401. Coulon A, Fitzpatrick JM, Bowman R, Stith BM, Makarewich CA, Stenzler LM, Lovette IJ (2008) Congruent population structure inferred from dispersal behaviour and intensive genetic surveys of the threathened Florida scrub-jay (Aphelocoma coerulenscens). Molecular Ecology, 17, 1685-1701. Davis JR (1980) Species composition and diversity of benthic macroinvertebrate population of the Pecos River, Texas. Southwestern Naturalist, 25, 241-256. Earl DA, vonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources, 4, 359-361. Edwards RJ, Garrett GP, Marsh-Matthews E (2002) Conservation and status of the fish communities inhabiting the Rio Conchos basin and middle Rio Grande, Mexico and U.S.A. Reviews in Fish Biology and Fisheries, 12, 119-132. Elith J, Graham CH, Anderson RP, Dudik M, Ferrier S, Guisan A, Hijmans RJ, Huettmann F, Leathwick JR, Lehmann A, Li J, Lohmann LG, Loiselle BA, Manion G, Moritz C, Nakamura M, Nakazawa Y, Overton JM, Peterson AT, Phillips SJ, Richardson K, Scachetti-Pereira R, Schapire RE, Soberon J, Williams S, Wisz MS, Zimmermann NE (2006) Novel methods improve prediction of species' distributions from occurrence data. Ecography, 29, 129-151. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611-2620. Fetzner JW, Jr., Crandall KA (2003) Linear habitats and the nested clade analysis: an empirical evaluation of geographic versus river distances using an Ozark crayfish (Decapoda: Cambaridae). Evolution, 57, 2101-2118. Fisher-Reid MC, Engstrom TN, Kuczynski CA, Stephens PR, Wiens JJ (2013) Parapatric divergence of sympatric morphs in a salamander: incipient speciation on Long Island? Molecular Ecology, 22, 4681-4694. Fordham DA, Brook BW, Moritz C, Nogues-Bravo D (2014) Better forecasts of range dynamics using genetic data. Trends in Ecology & Evolution, 29, 436-443. Galloway WE, Whiteaker TL, Ganey-Curry P (2011) History of Cenozoic North American drainage basin evolution, sediment yield, and accumulation in the Gulf of Mexico basin. Geosphere, 7, 938-973. Goudet J (1995) FSTAT (version 1.2): a computer program to calculate F-statistics. Journal of Heredity, 86, 485-486.

This article is protected by copyright. All rights reserved.

Accepted Article

Hare MP, Nunney L, Schwartz MK, Ruzzante DE, Burford M, Waples RS, Ruegg K, Palstra F (2011) Understanding and estimating effective population size for practical application in marine species management. Conservation Biology, 25, 438–449. Hewitt G (2000) The genetic legacy of the Quanternary ice ages. Nature, 405, 907-913. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25, 19651978. Hocutt CH, Denoncourt RF, Stauffer JR JR (1978) Fishes of the Greenbrier River, West Virginia, with drainage history of the central Appalachians. Journal of Biogeography, 5, 59-80. Holleley CE, Geerts PG (2009) Multiplex Manger 1.0: a cross-platform computer program that plans and optimizes multiplex PCR. Biotechniques, 46, 511-517. Hubbs C, Miller RR, Edwards RJ, Thompson KW, Marsh E, Garrett GP, Powell GL, Norris DJ, Zerr RW (1977) Fishes inhabiting the Rio Grande, Texas and Mexico, between El Paso and the Pecos confluence. In: Importance, Preservation and Management of Riparian Habitat: A Symposium eds. Johnson RR, Jones DA). USDA Forest Service, Fort Collins, Colorado. Inoue K, Lang BK, Berg DJ (2013) Development and characterization of 20 polymorphic microsatellite markers for the Texas hornshell, Popenaias popeii (Bivalvia: Unionidae), through next-generation sequencing. Conservation Genetics Resources, 5, 195-198. Inoue K, Levine TD, Lang BK, Berg DJ (2014a) Long-term mark-and-recapture study of a freshwater mussel reveals patterns of habitat use and an association between survival and river discharge. Freshwater Biology, 59, 1872-1883. Inoue K, Monroe EM, Elderkin CL, Berg DJ (2014b) Phylogeographic and population genetic analyses reveal Pleistocene isolation followed by high gene flow in a wide-ranging, but endangered, freshwater mussel. Heredity, 112, 282-290. Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23, 1801-1806. Jezkova T, Jaeger JR, Marshall ZL, Riddle BR (2009) Pleistocene impacts on the phylogeography of the desesrt pocket mouse (Chaetodipus penicillatus). Journal of Mammalogy, 90, 306-320. Kalinowski ST (2011) The computer program STRUCTURE does not reliably identify the main genetic clusters within species: simulations and implications for human population structure. Heredity, 106, 625-632. Karatayev AY, Miller TD, Burlakova LE (2012) Long-term changes in unionid assemblages in the Rio Grande, one of the world's top 10 rivers at risk. Aquatic Conservation: Marine and Freshwater Ecosystems, 22, 206-219. Kim D, Conway KW (2014) Phylogeography of Rhinichthys cataractae (Teleostei: Cyprinidae): pre-glacial colonization across the Continental Divide and Pleistocene diversification within the Rio Grande drainage. Biological Journal of the Linnean Society, 111, 317-333. Knowles LL (2009) Statistical phylogeography. Annual Review of Ecology, Evolution and Systematics, 40, 593-612. Kodric-Brown A, Brown JH (2007) Native fishes, exotic mammals, and the conservation of desert springs. Frontiers in Ecology and the Environment, 5, 549-553.

This article is protected by copyright. All rights reserved.

Accepted Article

Kramer-Schadt S, Niedballa J, Pilgrim JD, Schröder B, Lindenborn J, Reinfelder V, Stillfried M, Heckmann I, Scharf AK, Augeri DM, Cheyne SM, Hearn AJ, Ross J, Macdonald DW, Mathai J, Eaton J, Marshall AJ, Semiadi G, Rustam R, Bernard H, Alfred R, Samejima H, Duckworth JW, Breitenmoser-Wuersten C, Belant JL, Hofer H, Wilting A, Robertson M (2013) The importance of correcting for sampling bias in MaxEnt species distribution models. Diversity and Distributions, 19, 1366-1379. Kuhner MK (2006) LAMARC 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics, 22, 768-770. Lang BK (2001) Status of the Texas hornshell and native freshwater mussels (Unionoidea) in the Rio Grande and Pecos River of New Mexico and Texas. New Mexico Department of Game and Fish, Santa Fe, New Mexico. Levine TD, Lang BK, Berg DJ (2012) Physiological and ecological hosts of Popenaias popeii (Bivalvia: Unionidae): laboratory studies identify more hosts than field studies. Freshwater Biology, 57, 1854-1864. Librado P, Rozas J (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics, 25, 1451-1452. Luikart G, Ryman N, Tallmon DA, Schwartz MK, Allendorf FW (2010) Estimation of census and effective population sizes: the increasing usefulness of DNA-based approaches. Conservation Genetics, 11, 355-373. Lydeard C, Cowie RH, Ponder WF, Bogan AE, Bouchet P, Clark SA, Cummings KS, Frest TJ, Gargominy O, Herbert DG, Hershler R, Perez KE, Roth B, Seddon M, Strog EE, Thompson FG (2004) The global decline of nonmarine mollusks. Bioscience, 54, 321330. Mäkinen HS, Merilä J (2008) Mitochondrial DNA phylogeography of the three-spined stickleback (Gasterosteus aculeatus) in Europe - evidence for multiple glacial refugia. Molecular Phylogenetics and Evolution, 46, 167-182. Mantooth SJ, Hafner DJ, Bryson RW, Jr., Riddle BR (2013) Phylogeographic diversification of antelope squirrels (Ammospermophilus) across North American deserts. Biological Journal of the Linnean Society, 109, 949-967. Mayden RL (1988) Vicariance biogeography, parsimony, and evolution in North America freshwater fishes. Systematic Zoology, 37, 329-355. Metcalf AL (1997) Land snail of New Mexico from an historical zoogeographic point of view. New mexico Museum of Natural History and Science, Bulletin, 10, 71-108. Metcalf AL, Smartt RA (1997) Land snails of New Mexico. New mexico Museum of Natural History and Science, Bulletin, 10, 1-145. Metcalfe SE, O'Hara SL, Caballero M, Davies SJ (2000) Records of Late Pleistocene–Holocene climatic change in Mexico — a review. Quanternary Science Reviews, 19, 699-721. Miller RR (1977) Composition and derivation of the native fish fauna of the Chihuahuan Desert region. In: Transactions of the Symposium on the Biological Resources of the Chihuahuan Desert Region, United States and Mexico (eds. Wauer RH, Riskind DH), pp. 365-381. Sul Ross State University, Alpine, Texas. Nielsen R, Beaumont MA (2009) Statistical inferences in phylogeography. Molecular Ecology, 18, 1034-1047. Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, O'Hara RB, Simpson GL, Solymos P, Stevens MHH, Wagner H (2015) Vegan: community ecology package, R package version 2.2-1, available at http://cran.r-project.org/web/packages/vegan/.

This article is protected by copyright. All rights reserved.

Accepted Article

Olson DM, Dinerstein E (1998) The Global 200: a representation approach to conserving the Earth's most biologically valuable ecoregions. Conservation Biology, 12, 502-515. Otto-Bliesner BL, Marshall SJ, Overpeck JT, Miller GH, Hu A (2006) Simulating Arctic climate warmth and icefield retreat in the last interglaciation. Science, 311, 1751-1753. Palstra FP, Fraser DJ (2012) Effective/census population size ratio estimation: a compendium and appraisal. Ecology and Evolution, 2, 2357-2365. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6, 288-295. Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190, 231-259. Phillips SJ, Dudik M, Elith J, Graham CH, Lehmann A, Leathwick JR, Ferrier S (2009) Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecological Applications, 19, 181-197. Polyak VJ, Asmerom Y (2001) Late Holocene climate and cultural changes in the southwestern United States. Science, 294, 148-151. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genetype data. Genetics, 155, 945-959. R Development Core Team (2011) R: A language and environment for statistical computing. R Fundation for Statistical Computing, Vienna, Austria. Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Molecular Ecology Notes, 4, 137-138. Rousset F (2008) GENEPOP'007: a complete re-implementation of the GENEPOP software for Windows and Linux. Molecular Ecology Resources, 8, 103-106. Schmandt J (2002) Bi-national water issues in the Rio Grande/Río Bravo basin. Water Policy, 4, 137-155. Schwalb AN, Morris TJ, Mandrak NE, Cottenie K (2013) Distribution of unionid freshwater mussels depends on the distribution of host fishes on a regional scale. Diversity and Distributions, 19, 446-454. Seidel RA, Lang BK, Berg DJ (2009) Phylogeographic analysis reveals multiple cryptic species of amphipods (Crustacea: Amphipoda) in Chihuahuan Desert springs. Biological Conservation, 142, 2303-2313. Smith ML, Miller RR (1986) The evolution of the Rio Grande basin as inferred from its fish fauna. In: The Zoogeography of North American Freshwater Fishes (eds. Hocutt CH, Wiley EO). John Wiley & Sons, New York. Snook RR, Brustle L, Slate J (2009) A test and review of the role of effective population size on experimental sexual selection patterns. Evolution, 63, 1923-1933. Strange RM, Burr BM (1997) Intraspecific phylogeography of North American highland fishes: A test of the Pleistocene vicariance hypothesis. Evolution, 51, 885-897. Strenth NE (2004) New records of the Texas hornshell Popenaias popeii (Bivalvia: Unionidae) from Texas and northern Mexico. The Texas Journal of Science, 56, 223. Thomas RG (1972) The Geomorphic Evolution of the Pecos River System The Baylor University Press, Waco, TX. Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3—new capabilities and interfaces. Nucleic Acids Research, 40, e115. US Fish and Wildlife Service (2013) Endangered and threatened wildlife and plants; review of native species that are cadidates for listing as endangered or threatened; annual notice of

This article is protected by copyright. All rights reserved.

Accepted Article

findings on resubmitted petitions; annual description of progress on listing actions. Federal Register, 78, 70104-70162. van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: Software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4, 535-538. Warren ML, Pardew MG (1998) Road crossings as barriers to small-stream fish movement. Transactions of the American Fisheries Society, 127, 637-644. Williams JD, Warren ML, Jr., Cummings KS, Harris JL, Neves RJ (1993) Conservation status of freshwater mussels of the United States and Canada. Fisheries, 18, 6-22. Wilson GA, Rannala B (2003) Bayesian inference of recent migration rates using multilocus genotypes. Genetics, 163, 1177-1191. Wong CM, Williams CE, Pittock J, Collier U, Schelle P (2007) World's top 10 rivers at risk. WWF International, Gland, Swizerland.

Data Accessibility COX sequences are deposited in GenBank: accession numbers KP779655 – KP779900. Specimen information; aligned mtDNA sequences, in FASTA format; microsatellite genotype data, in GENEPOP format; input and parameter files for STRUCTURE; and occurrence records (both raw and MAXENT-formatted files) for P. popeii and its hosts for the ENMs are available on Dryad, doi:10.5061/dryad.pm870. Author Contributions KI designed the study, conducted field sampling, performed laboratory and data analyses, and drafted the manuscript. BKL designed the study and conducted field analyses. DJB designed the study and contributed to field sampling, data analysis, and writing of the manuscript.

Figure legends Fig. 1. Map of the Rio Grande drainages in the southwest USA indicating major tributaries, reservoirs, and Popenaias popeii sampling sites (circles). Colors correspond to the parsimony network of COX sequences in Fig. 3. The magnified inset map shows Black River sampling locations. Fig. 2. Four demographic scenarios tested using DIYABC: 1) all populations have constant size over time; 2) a small number of individuals founded the ancestral Black River population followed by population expansion at te2; 3) a large number of founders colonized the ancestral Black River population followed by a population bottleneck at te2; and 4) Rio Grande populations experienced population bottleneck at te2 followed by population expansion at te1, in addition to Scenario 3. All scenarios assume divergence between the ancestral BR and RG populations at t2, divergence between BR-u and BR-d at t1, and STRUCTURE-defined populations at the present time (t0). Timing of events is shown at right. See Table 2 for detailed labels. Fig. 3. Parsimony network of COX sequences for Popenaias popeii. Each circle represents a unique haplotype; lines between haplotypes represent single-base-pair changes; black dots are inferred extinct or unsampled haplotypes. Haplotype frequency is relative to the size and number in the circle. The most common haplotype is shared between BR (n = 185) and RG (n =

This article is protected by copyright. All rights reserved.

Accepted Article

10) individuals. Circle colors represent localities (Black River = dark gray; Devils River = light gray; Rio Grande = white).

Fig. 4. Bar plots obtained from STRUCTURE, assigning individuals into k = 2 and k = 3 clusters. For k = 2, clusters consisted of BR (light gray) and RG (black) individuals. For k = 3, an additional division was observed in the BR populations (light gray, upstream sites; medium gray, downstream sites). Site name abbreviations at the bottom of the plot are a priori population assignments.

Fig. 5. Potential distribution of Popenaias popeii identified using ecological niche modeling under current bioclimatic conditions (1950 to 2000) and projections back to two paleoclimatic periods (last interglacial, 120 – 140 ka; last glacial maximum, 21 ka). Models included major rivers of the Rio Grande watershed in the USA and Mexico. Black dots on the Present map represent occurrence points included in the ENMs. Scale bars in the bottom-left corners represent 200 km.

Table 1 Descriptive statistics for COX sequences and 18 microsatellite loci for Popenaias popeii populations in the Black River, New Mexico and the Devils River and Rio Grande, Texas. COX Microsatellites HR K π N A AR N P HO HE 3 1 -0 0 2.0 -0 0.380 0.314 24 1 1.0 0 0 3.3 2.8 0 0.469 0.483 19 1 1.0 0 0 3.6 2.9 1 0.469 0.477 55 1 1.0 0 0 4.5 3.1 2 0.503 0.507 20 1 1.0 0 0 3.9 3.1 0 0.527 0.518 41 1 1.0 0 0 4.6 3.1 4 0.507 0.532 10 1 1.0 0 0 3.4 3.0 0 0.450 0.480 21 1 1.0 0 0 3.9 3.1 0 0.513 0.518 193 1 1.0 0 0 6.0 3.0 8 0.495 0.528 3 Devils 1 -0 0 3.4 -2 0.556 0.608 29 21 5.9 3.56 0.0050 14.2 7.7 Rio Grande 40 0.813 0.881 5 4 -- 2.00 0.0028 6.2 -5 0.800 0.793 2 2 -- 5.00 0.0070 3.4 -0 0.800 0.693 16 13 5.6 3.39 0.0047 11.4 7.6 11 0.805 0.872 6 5 5.0 4.60 0.0064 6.7 6.7 3 0.731 0.798 58 33 5.5 3.49 0.0049 16.7 7.3 167 0.806 0.902 254 34 2.2 1.41 0.0020 17.3 4.5 Overall 68 0.566 0.694 AR, rarefied allelic richness; COX, ~810 basepairs of part of the cytochrome oxidase II gene and the cytochrome oxidase I gene; H, number of haplotypes; HE, mean expected heterozygosity; HO, mean observed heterozygosity; HR, rarefied number of haplotypes; K, mean number of base pair differences between all pairs of individuals; n, number of mussels sampled; NA, mean number of observed alleles; NP, number of private alleles; π, nucleotide diversity River Black

Site ID CC Fall Magby DV RF5 RF1 BS RF3 All BR DR RG1 RG2 RG3 RG4 RG5 All RG

n

H

This article is protected by copyright. All rights reserved.

Accepted Article

Table 2 Prior distributions of parameters for each scenario (Fig. 2) and posterior parameter values for Scenario 3. The unit of time is generations, except that values in parentheses are 1000 years before present (ka) calculated with an average generation time of 8.9 years (based on 8.1 9.6 years for Lampsilis radiata; Chagnon & de la Cheneliere 1998).

Population Black River upstream population (BR-u) Black River downstream population (BR-d) Rio Grade population (RG) Ancestral Black River population Ancestral Rio Grande population Founder population into the Black River Bottlenecked population in the Black River Bottlenecked population in the Rio Grande Most ancestral population Time of divergence within the Black River

Parameter NBR-u NBR-d NRG NBR NA-RG NF-BR NB-BR NB-RG NA t1

Time of expansion/bottleneck event Time of expansion/bottleneck event

te1 te2

Time of the ancestral divergence

t2

Mean mutation rate per generation

µ

Prior distribution Distribution ~U(10,5000) ~U(10,50000 ~U(10,100000) ~U(10,50000) ~U(10,100000) ~U(10,50000) ~U(10,5000) ~U(10,100000) ~LU(10,100000) ~U(1,500)

Posterior parameter estimates Median 95% credible interval RMAE 1190 194 4230 0.212 4680 3460 4980 0.222 70300 41600 91900 0.200 39400 5800 87200 0.422 --------1630 344 4490 0.254 ----81000 45900 99100 0.170 125 12 410 0.298 (1.113 ka) (0.107 ka) (3.649 ka) ~U(1,20000) ----~U(1,20000) 3670 614 11000 0.316 (32.663 ka) (5.465 ka) (97.900 ka) ~U(1,20000) 9020 2370 18900 0.185 (80.270 ka) (21.093 ka) (168.210 ka) ~U(1×10-5,1×10-3) 1.53×10-4 7.48×10-5 3.39×10-4 0.214

LU = log-uniform distribution; RMAE = relative median of absolute error; U = uniform distribution

Table 3 Pairwise FST (below diagonal) and DEST (above diagonal) for 18 microsatellite loci from STRUCTURE-defined populations of Popenaias popeii. All values are statistically significantly different from 0 at P