Fundulus stellifer Coosa River

Fundulus stellifer Coosa River View in CoL

1 East Fork Amite River 31 ° 05 0 55 ″ 90 ° 43 0 11 ″

2 Homochitto River 31 ° 39 0 47 ″ 90 ° 43 0 14 ″

3 Caddo River 34 ° 23 0 10 ″ 93 ° 36 0 30 ″

4 South Fork Ouachita River 34 ° 33 0 29 ″ 93 ° 41 0 46 ″

5 Ten Mile Creek 34 ° 32 0 43 ″ 92 ° 45 0 15 ″

6 Elk River 36 ° 37 0 25 ″ 94 ° 35 0 25 ″

7 Elk River 36 ° 33 0 37 ″ 94 ° 25 0 26 ″

8 Richland Creek 36 ° 02 0 49 ″ 93 ° 58 0 19 ″

9 Kings River 36 ° 23 0 39 ″ 93 ° 38 0 10 ″ 10 Strawberry River 36 ° 05 0 56 ″ 91 ° 36 0 32 ″ 11 Black River 37 ° 25 0 00 ″ 90 ° 49 0 31 ″ 12 Current River 36 ° 37 0 04 ″ 90 ° 50 0 22 ″ 13 Big River 37 ° 48 0 46 ″ 90 ° 46 0 20 ″ 14 Huzzah Creek 37 ° 56 0 52 ″ 91 ° 10 0 39 ″ 15 Big Piney River 37 ° 13 0 11 ″ 92 ° 00 0 17 ″ 16 Leatherwood Creek 39 ° 34 0 56 ″ 85 ° 59 0 05 ″ 17 Sugar Creek 39 ° 37 0 26 ″ 85 ° 56 0 46 ″ 18 Trammel Fork 36 ° 45 0 08 ″ 86 ° 17 0 15 ″ 19 Falling Timber Creek 36 ° 55 0 27 ″ 85 ° 48 0 19 ″ 20 Hurricane Creek 36 ° 35 0 58 ″ 85 ° 40 0 08 ″ 21 Turnbull Creek 36 ° 06 0 03 ″ 87 ° 07 0 35 ″ 22 East Fork Stones River 35 ° 56 0 30 ″ 86 ° 22 0 36 ″ 23 Otter Creek 36 ° 42 0 47 ″ 84 ° 57 0 42 ″ 24 Bear Creek 34 ° 38 0 01 ″ 88 ° 09 0 21 ″ 25 Buffalo River 35 ° 27 0 48 ″ 87 ° 32 0 07 ″ 26 Duck River 35 ° 28 0 59 ″ 86 ° 27 0 47 ″ 27 Little River 35 ° 47 0 07 ″ 83 ° 53 0 01 ″ 28 Little Pigeon River 35 ° 52 0 12 ″ 83 ° 34 0 03 ″ 29 Clinch River 36 ° 31 0 25 ″ 83 ° 09 0 20 ″ 30 Hachemedega Creek 32 ° 50 0 42 ″ 86 ° 13 0 36 ″

5 5 5 5 6 5 5 5 5 4 3 3 5 5

10 5 5 3 3 3 5 3 2 5 5 4 4 4 2 2

1

1

1

1

1 1

1 1

1

1

2 1 1 1 1

2

Individuals nDNA) (2 1 1 1

Individuals) (cytb 2 1 2 1 2 2 1 1 2

″ Longitude (W) ° ″ 85 36 51 ° ″ 35 92 13 ° ″ 59 93 39 ° ″ 79 17 57 ° ″ 80 49 49 ″ ° 34 81 40 0 0 0 0 0 0 0 ° 101 48 57

″ ″ ″ ″ ″ ″ ″ Latitude (N) 0 ° 59 08 32 0 ° 34 10 19 ° 0 47 14 46 ° 0 36 11 21 ° 0 11 27 51 0 ° 29 05 12 0 ° 24 38 56

,) 445402 Fisheries Creek) Creek River Josie Leg Creek River Saline Lacs Mille Lake FJ GenBank (Conservation. Charles (Inc Country Line Lake Okeechobee George Lake Willow Creek

Locality 31 Outgroup Outgroup Outgroup Outgroup Outgroup Outgroup Outgroup Outgroup

River Drainage Tallapoosa

Continued bifax chrysotus diaphanus heteroclitus julisia rathbuni seminolis zebrinus Table. 1 Species Fundulus Fundulus Fundulus Fundulus Fundulus Fundulus Fundulus Fundulus template DNA, 2.75 lL water, 6.25 lL GoTaq Green Master Mix (Promega, Madison, WI), 1.0 lL forward primer, and 1.0 lL reverse primer. Cytb was amplified using the GLU (5 0 -GACTTGAAGAACCA CCGTTG-3 0) and THR (5 0 -TCCGACATTCGGTTTAG AAG-3 0) primers described in Near, Porterfield & Page (2000) and using the following thermal conditions: denaturation at 94 ° C for 30 s, annealing at 50 ° C for 30 s, and extension at 72 ° C for 90 s, repeated for 25 cycles.

The four nDNA intron loci were amplified using the following primers: stx5a – 10F (5 0 -GGAGGAG ACKGACTGGAAGT-3 0) and 10R (5 0 -GCAGAACATY GARAGCACMA-3 0), ncl1 – 428 bp F8 (5 0 -SGCCAGG TTGATYTTCTTRT-3 0) and R8 (5 0 -CCAGTCTGCTSC AGGACAAY-3 0), rpsa – 3F (5 0 -ATTGTTGCCATYG ARAAYCC-3 0) and 3R (5 0 -GCWGCCTGRATCTGA TTGGT-3 0), and rps3 – F4 (5 0 -CTACAAGCTGCTSGG AGGMC-3 0) and R4 (5 0 -TAGTTSACKGGGTCTCC RCT-3 0) ( Halas & Simons, 2014, D. Halas pers. comm.). All introns were amplified with an initial denaturation at 95 ° C for 4 min, followed by denaturation at 95 ° C for 40 s, varying annealing temperatures for 40 s, and 7 ° C for 90 s for 25 cycles. Annealing temperature for ncl1, stx5a, rpsa, and rps3 was 50, 51, 52, and 56 ° C, respectively. Owing to initial difficulty phasing rpsa and rps3 sequences of some individuals, internal primers were developed using PRIMER3 (http://frodo.wi.mit.edu/primer3/). Additional rpsa sequences were amplified using the intron thermal profile with the primer pairs F101 (5 0 -GTAAACGGATCGGGGTTTCT-3 0) and R800 (5 0 - AAGGCCCTTTTCACTTTTCA-3 0), and F146 (5 0 -TG ACTGGGGTATGAGAAGCTC-3 0) and R728 (5 0 -CA CGCTTTCTAACCTCCCTTT-3 0) and an annealing temperature of 50 and 54 ° C, respectively.

PCR products were purified using Exonuclease 1 and shrimp alkaline phosphatase (USB Corporation, Cleveland, OH) at the manufacturer’s suggested thermal profiles. Automated Sanger sequencing of purified PCR products was performed using ABI Prism BigDye Terminator v. 3.1 chemistry (Applied Biosystems, Foster City, CA) at the Biomedical Genomics Center DNA Sequencing and Analysis Facility at the University of Minnesota.

EDITING AND ALIGNING SEQUENCES

Complementary heavy and light strands were aligned into contiguous sequences (contigs) and edited in GENEIOUS v. 6.1.6 (www.geneious.com; Biomatters Ltd., Auckland, New Zealand). Length heterozygotes found in nuclear introns were phased by eye with help from CHAMPURU v. 1.0 ( Flot, 2007, available online at http://www.mnhn.fr/jfflot/ champuru/). Consensus sequences of contigs were aligned using the MUSCLE ( Edgar, 2004) clustering algorithm as implemented in GENEIOUS v. 6.1.6. Each alignment was trimmed to make sequences near uniform in length. In intron alignments, large indels found only in outgroups were removed. This reduced the length of rps3 and stx5a. Nuclear introns were tested for recombination using the phi test ( Bruen, Philippe & Bryant, 2006) as implemented in SplitsTree v. 4.13.1 ( Huson & Bryant, 2006). For some analyses, identical cytb haplotypes were removed using ElimDupes (http://hcv.lanl. gov/content/sequence/ELIMDUPES/elimdupes.html).

PHYLOGENETIC ANALYSIS

Two data sets were created for cytb: all individuals and unique haplotypes. The best-fitting partitioning scheme and nucleotide substitution models were determined using PartitionFinder v. 1.01 ( Lanfear et al., 2012) based on Bayesian information criterion scores. Bayesian analysis of each cytb data set was conducted using MrBayes v. 3.2.1 ( Ronquist et al., 2012) on the CIPRES Science Gateway portal ( Miller, Pfeiffer & Schwartz, 2010). The Metropolis coupled Markov chain Monte Carlo (MCMCMC) command was used for two simultaneous runs with four chains (three heated chains, one cold) per 15 000 000 generations, sampling every 1000. Log files were checked in the program TRACER v. 1.5 (http://beast. bio.ed.ac.uk/Tracer) to assess convergence of runs and burn-in was set to remove the first 20% of sampled trees. We performed maximum likelihood (ML) analyses in the program GARLI v. 2.0 ( Zwickl, 2006) on the CIPRES Science Gateway portal ( Miller et al., 2010). The gene tree with the best likelihood score was selected from five search replicates. The nodes of the best ML tree found by the aforementioned five search replicates were annotated with the proportion of nodes found by 1000 bootstrap replicates using SumTrees v. 3.3. 1 in the DendroPy v. 3.11 package ( Sukumaran & Holder, 2010).

To account for genetic diversity observed in the cytb gene tree, we subsampled six F. catenatus, ten Fundulus sp. cf. catenatus , two F. bifax and F. stellifer , and one F. julisia, F. diaphanus , and F. rathbuni for the nuclear intron data set. Phylogenetic analysis of each of the four nuclear intron data sets followed the same methods of partitioning scheme and nucleotide substitution model selection, Bayesian, and ML analyses.

SPECIES TREE ANALYSIS

Species tree analysis of the subsampled individuals (six Fundulus catenatus , ten F. sp. cf. catenatus , two F. bifax and F. stellifer , and one F. julisia, F. diaphanus , and F. rathbuni ) was conducted using *BEAST v. 1.7.5 ( Heled & Drummond, 2010). Species were designated based on the hypothesis of relationship suggested by the cytb gene tree. The following methods were used for an all loci (cytb + nuclear introns) and just nuclear intron data set. An appropriate clock model was determined by performing molecular clock likelihood ratio tests for each locus in PAUP* v. 4.0b10 ( Swofford, 2003). A Yule process speciation prior was used for branching rates. We applied similar partitioning schemes and nucleotide substitution models to analyses of individual genes. Ten independent runs of 50 000 000 generations each were conducted, sampling every 1000 generations. The MCMCMC log files were analysed in TRACER v. 1.5 (http://beast.bio.ed.ac.uk/Tracer) to assess convergence of the runs, ensure proper mixing, and determine an appropriate burn-in (first 10% of sampled trees). LogCombiner v. 1.7.4 (http://beast.bio.ed.ac.uk/ LogCombiner) was used to remove burn-in and combine files; the resulting 10 000 trees were used to produce a maximum clade credibility tree using TreeAnnotator v. 1.7.4 (http://beast.bio.ed.ac.uk/ TreeAnnotator).

SPECIES DELIMITATION USING NDNA

To assess the deep split recovered in both the cytb gene tree and species tree between Fundulus catenatus and F. sp. cf. catenatus , we used the program BAYESIAN PHYLOGENETICS AND PHYLOGEOGRAPHY (BP&P v. 2.2: Yang & Rannala, 2010) to compare a one-species model ( F. catenatus combined with F. sp. cf. catenatus ) and a two-species model ( F. catenatus and F. sp. cf. catenatus ). In BP&P, we used the reversible-jump Markov chain Monte Carlo method ( Rannala & Yang, 2013) to delimit species and assess cryptic diversity. BP&P requires a guide tree, species group membership definitions, and sequence alignments. The guide tree and species group memberships were assigned according to our hypothesis of species limits based on the cytb gene tree ( Fig. 2 View Figure 2 ). The phased nuclear introns of six F. catenatus and ten F. sp. cf. catenatus individuals were provided as the sequence alignments. To evaluate the influence of some priors and settings we ran multiple runs for 500 000 generations, sampling every five generations, and used a burn-in of 50 000 generations. We considered three different combinations of prior distributions of two parameters known to influence the posterior probability for models, ancestral population size (Θ) and root age (s) ( Yang & Rannala, 2010): large ancestral population size with deep divergence amongst species, small ancestral population size with shallow divergence amongst species, and large ancestral population size with small divergence amongst species ( Leache & Fujita, 2010). As suggested by the user’s manual we set the mutation rate parameter, ‘locusrate’, at 15 to account for similar rates amongst noncoding loci and ran analysis utilizing both provided reversible-jump Markov chain Monte Carlo algorithms to ensure similar results. A posterior probability ≥ 0.95 at the splitting event of F. catenatus and F. sp. cf. catenatus was considered strong evidence of two species (following Leache & Fujita, 2010).

To further explore evidence for multiple species, we utilized a species delimitation method using Bayes factors described by Grummer, Bryson & Reeder (2014). This method compares marginal likelihood scores, estimated with stepping-stone sampling ( Xie et al., 2011), and path sampling ( Lartillot & Philippe, 2006), using Bayes factors. We estimated species trees for both species delimitation models using the methods presented in the Species Tree Analysis section. The sole difference was that the initial species tree inference was run for 75 000 000 generations, sampling every 1000 generations. Following *BEAST analysis, the initial 20% of sampling was removed as burn-in, and both path sampling and stepping-stone sampling were executed for a chain length of 7 500 000 generations for 300 paths (totalling 225 000 000 generations). This was carried out five times to provide evidence of consistent results. Resulting marginal likelihood scores were compared as 2Ln (Bayes factors) (where Bayes factor = marginal likelihood score of one-species model – marginal likelihood score of a two-species model). Kass & Raftery (1995) suggested considering 2Ln (Bayes factors) scores of 0 – 2, 2 – 6, 6 – 10, and> 10 as ‘not worth more than a bare mention’, ‘positive’ support, ‘strong’ support, and ‘very strong’ support, respectively.

POPULATION STATISTICS AND MEDIAN- JOINING NETWORK

The following summary statistics were assessed using DnaSP v. 5.10.1 ( Librado & Rozas, 2009): singletons, polymorphic sites, parsimony informative sites, number of haplotypes, haplotype diversity, and nucleotide diversity.

A cytb haplotype network was constructed by the median-joining method (Bandelt, Forster & Rohl, 1999) in the program NETWORK v. 4.6.1.2 (fluxusengineering.com). The data set contained all individuals sampled and was trimmed to the shortest sequences (934 bp) because large amounts of missing data in a sequence led to a large number of unverifiable ‘unique’ haplotypes. All default settings were used with the one exception of the switching parameters frequency> 1 criterion being set to active.

NEUTRALITY TESTS

We explored possible population expansion using a coalescent-based approach with Fu’s F S ( Fu, 1997) and Ramos-Onsins and Rozas’ R 2 ( Ramos-Onsins & Rozas, 2002). We used both methods as we have varying sample sizes, and tests run by Ramos-Onsins & Rozas (2002) found R 2 to function better for small sample sizes, whereas F S behaved best for large sample sizes. The following populations were selected based on median-joining network: F. sp. cf. catenatus (excluding Mississippi and Ouachita), F. sp. cf. catenatus (Ouachita) , and F. catenatus . The southern Mississippi population of F. sp. cf. catenatus was excluded from this analysis owing to limited sampling (ten individuals). F S and R 2 were calculated and coalescent simulations were run given segregating sites for 10 000 replicates in the program DnaSP v. 5.10.1. Significantly negative values of F S (P -value <0.02) and significant R 2 (<0.05) suggest an excess of rare haplotypes, indicative of non-neutral processes such as recent demographic expansion or genetic hitchhiking ( Fu, 1997; Ramos-Onsins & Rozas, 2002).