Citation

Material Information

Title:
Adding loci improves phylogeographic resolution in red mangroves despite increased missing data: comparing microsatellites and RAD-Seq and investigating loci filtering
Creator:
Hodel, Richard
Publisher:
Nature
Publication Date:
Language:
English
Physical Description:
Journal Article

Notes

Abstract:
The widespread adoption of RAD-Seq data in phylogeography means genealogical relationships previously evaluated using relatively few genetic markers can now be addressed with thousands of loci. One challenge, however, is that RAD-Seq generates complete genotypes for only a small subset of loci or individuals. Simulations indicate that loci with missing data can produce biased estimates of key population genetic parameters, although the influence of such biases in empirical studies is not well understood. Here we compare microsatellite data (8 loci) and RAD-Seq data (six datasets ranging from 239 to 25,198 loci) from red mangroves (Rhizophora mangle) in Florida to evaluate how different levels of data filtering influence phylogeographic inferences. For all datasets, we calculated population genetic statistics and evaluated population structure, and for RAD-Seq datasets, we additionally examined population structure using coalescence. We found higher F ST using microsatellites, but that RAD-Seq-based estimates approached those based on microsatellites as more loci with more missing data were included. Analyses of RAD-Seq datasets resolved the classic Gulf-Atlantic coastal phylogeographic break, which was not significant in the microsatellite analyses. Applying multiple levels of filtering to RAD-Seq datasets can provide a more complete picture of potential biases in the data and elucidate subtle phylogeographic patterns.
Acquisition:
Collected for University of Florida's Institutional Repository by the UFIR Self-Submittal tool. Submitted by Richard Hodel.

Record Information

Source Institution:
University of Florida Institutional Repository
Holding Location:
University of Florida
Rights Management:
Copyright Creator/Rights holder. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.

Downloads

This item is only available as the following downloads:


Full Text

PAGE 1

1 Adding loci improves phylogeographic resolution in red mangroves despite increased missing data: comparing microsatellites and RAD-Seq and ‹‡•–‹‰ƒ–‹‰Ž‘…‹¤Ž–‡”‹‰Richard G. J. Hodel n n wxŠ‹…Šƒ‘ Chenxy†ƒ Pa ytonw–—ƒ”–t McDanielwxz Pamela Soltis n n xz & Douglas E. SoltiswxzThe widespread adoption of RAD-Seq data in phylogeography means genealogical relationships previously evaluated using relatively few genetic markers can now be addressed with thousands of Ž‘…‹‡…ŠƒŽŽ‡‰‡Š‘™‡‡”‹•–Šƒ–‡“‰‡‡”ƒ–‡•…‘’Ž‡–‡‰‡‘–›’‡•ˆ‘”‘Ž›ƒ•ƒŽŽ•—„•‡– of loci or individuals. Simulations indicate that loci with missing data can produce biased estimates of ‡›’‘’—Žƒ–‹‘‰‡‡–‹…’ƒ”ƒ‡–‡”•ƒŽ–Š‘—‰Š–Š‡‹—‡…‡‘ˆ•—…Š„‹ƒ•‡•‹‡’‹”‹…ƒŽ•–—†‹‡•‹•‘– ™‡ŽŽ—†‡”•–‘‘†‡”‡™‡…‘’ƒ”‡‹…”‘•ƒ–‡ŽŽ‹–‡†ƒ–ƒ~Ž‘…‹fƒ†‡“†ƒ–ƒ•‹š†ƒ–ƒ•‡–•”ƒ‰‹‰ ˆ”‘xy–‘x{w~Ž‘…‹fˆ”‘”‡†ƒ‰”‘‡•Rhizophora mangle f‹tŽ‘”‹†ƒ–‘‡ƒŽ—ƒ–‡Š‘™†‹¡‡”‡– Ž‡‡Ž•‘ˆ†ƒ–ƒ¤Ž–‡”‹‰‹—‡…‡’Š›Ž‘‰‡‘‰”ƒ’Š‹…‹ˆ‡”‡…‡•t‘”ƒŽŽ†ƒ–ƒ•‡–•™‡…ƒŽ…—Žƒ–‡†’‘’—Žƒ–‹‘ ‰‡‡–‹…•–ƒ–‹•–‹…•ƒ†‡ƒŽ—ƒ–‡†’‘’—Žƒ–‹‘•–”—…–—”‡ƒ†ˆ‘”‡“†ƒ–ƒ•‡–•™‡ƒ††‹–‹‘ƒŽŽ› ‡šƒ‹‡†’‘’—Žƒ–‹‘•–”—…–—”‡—•‹‰…‘ƒŽ‡•…‡…‡‡ˆ‘—†Š‹‰Š‡” FST—•‹‰‹…”‘•ƒ–‡ŽŽ‹–‡•„—– that RAD-Seq-based estimates approached those based on microsatellites as more loci with more missing data were included. Analyses of RAD-Seq datasets resolved the classic Gulf-Atlantic coastal ’Š›Ž‘‰‡‘‰”ƒ’Š‹…„”‡ƒ™Š‹…Š™ƒ•‘–•‹‰‹¤…ƒ–‹–Š‡‹…”‘•ƒ–‡ŽŽ‹–‡ƒƒŽ›•‡•’’Ž›‹‰—Ž–‹’Ž‡ Ž‡‡Ž•‘ˆ¤Ž–‡”‹‰–‘‡“†ƒ–ƒ•‡–•…ƒ’”‘‹†‡ƒ‘”‡…‘’Ž‡–‡’‹…–—”‡‘ˆ’‘–‡–‹ƒŽ„‹ƒ•‡•‹–Š‡ data and elucidate subtle phylogeographic patterns. Choice of molecular markers remains a critically important consideration when designing a phylogeographic, phylogenetic, or population genetic study, as researchers must optimize the amount of informative genetic data they can obtain for a xed and typically modest cost. In phylogeographic studies, theoretical considerations impact decisions regarding whether to include more individuals or more loci. Microsatellites (or simple sequence repeats, SSRs) have been one of the workhorses of phylogeographic studies for over two decades—their high var iability made them popular for distinguishing between closely related conspecic or congeneric individuals1 – 3. Microsatellite markers are now being gradually replaced by RAD-Seq data for phylogeographic inference4. ere are advantages and disadvantages to using microsatellites in phylogeographic studies2 3 5. Microsatellites are a known quantity; hundreds of thousands of studies that use SSRs are in the literature—primers are already available for many groups. In addition, many user-friendly soware packages are available for all aspects of microsatellite analysis, from loci development to population genetic inference3. If primers are already developed for the taxa of interest, microsatellites can be inexpensive to implement. Additionally, if initial results necessitate adding a few additional individuals and/or loci, project costs will increase linearly with microsatellites. However, w‡’ƒ”–‡–‘ˆ‹‘Ž‘‰›‹‡”•‹–›‘ˆtŽ‘”‹†ƒnƒ‹‡•‹ŽŽ‡tyx|wwxFlorida Museum of Natural History ‹‡”•‹–›‘ˆtŽ‘”‹†ƒnƒ‹‡•‹ŽŽ‡tyx|wwy‘ŽŽ‡‰‡‘ˆ‹ˆ‡…‹‡…‡•ƒ†‡…Š‘Ž‘‰›‘‰Œ‹‹‡”•‹–› Šƒ‰Šƒ‹xvvvxŠ‹ƒzŠ‡n‡‡–‹…•f•–‹–—–‡‹‡”•‹–›‘ˆtŽ‘”‹†ƒnƒ‹‡•‹ŽŽ‡tyx|wv‘””‡•’‘†‡…‡ ƒ†”‡“—‡•–•ˆ‘”ƒ–‡”‹ƒŽ••Š‘—Ž†„‡ƒ††”‡••‡†–‘nr‡ƒ‹Ž ”‹…Š‹‡Š‘†‡Ž;‰ƒ‹Ž…‘ ) Received: 14 July 2017 Accepted: 15 November 2017 Published: xx xx xxxxOPEN

PAGE 2

2 there are caveats to using SSRs. Perhaps most importantly, a limited number of loci (usually t t 25) can feasibly be employed in a typical microsatellite study. Also, the mutational properties of SSRs are unusually high and almost certainly do not reect those of the genome as a whole. us, the property that makes microsatellites excellent for distinguishing dierent individuals may inate statistics such as FST and heterozygosity relative to the rest of the genome. Furthermore, microsatellites can be just as expensive to implement as newer high-throughput sequencing (HTS) techniques if there are no existing genetic resources (e.g., no primers already developed, or no available transcriptomic or genomic resources)6. e use of RAD-Seq data has increased greatly over the past decade, largely because thousands of loci can be generated simultaneously for hundreds of individuals for a xed, known cost7. RAD-Seq uses restriction enzymes (REs) to create a reduced representation library of the genome; single-nucleotide polymorphisms (SNPs) in regions of DNA between restriction sites are used to distinguish between individuals8. Barcoding to allow e cient multiplexing during sequencing keeps costs down, which can be as little as $40 per individual for thousands of loci, assuming judicious sharing of reagents, and a well-designed plan for multiplexing individuals9 – 11. Microsatellite genotyping has a similar cost per individual, assuming primers are not developed, but many fewer loci are obtained3. SNPs have several advantages over microsatellites, as they are less likely to exhibit homoplasy than SSRs12. Despite advantages, there are also several caveats to using RAD-Seq. Unless there is a reference genome, loci obtained using RAD-Seq are anonymous, and some loci may be non-neutral7. Additionally, biases may be introduced at several stages in a RAD-Seq protocol: (1) digestion with REs samples a non-random portion of the genome due to biases in base composition; this is potentially worse if methylation sensitive enzymes are used; (2) polymorphisms in restriction sites that can lead to segregating presence/absence polymorphisms that are very difcult to detect without very deep sequencing and negating the cost-savings of using RAD-Seq in the rst place7 13; (3) preferential PCR amplication of some loci during library construction necessarily reduces coverage of other loci13; (4) sequencing errors and/or low sequencing depth leads to incorrect genotype calling7; and (5) false loci are constructed due to the misassembly of paralogous reads14, 15. Many potential problems are resolved by multiple PCR steps to even out loci coverage and by improvements in soware when processing loci, but concerns remain that RE-based reduced representation methods do not capture a representative snapshot of the genome16. One other concern with RAD-Seq loci is that manual data curation is impossible, and errors may go undetected even by the most careful researchers14,17,18. Finally, the biggest potential problem when using RAD-Seq is that low coverage and high proportions of missing data can make it dicult to infer heterozygotes accurately. Previous studies have compared results from SNPs and SSRs, revealing that microsatellites provide much more information—up to an order of magnitude more—on a per-marker basis than SNPs19, 20. However, SNP studies typically use several orders of magnitude more markers than an average SSR study. Evidence has shown that the large number of loci in SNP studies can eectively allow for more powerful inferences, even though the information at each locus is less than that in microsatellite markers21. Because of the low number of loci used in SSR studies, the standard practice is to aim to minimize missing data. However, the nature of current library preparation and sequencing means that higher percentages of missing data are an unavoidable part of RAD-Seq studies. Simulation studies have shown that the large amounts of missing data in RAD-Seq studies can inate FST estimates due to allelic dropout13, 18. As more loci were included in these simulations, FST appeared to increase because many loci had genotype data for only one or a few individuals. In many such loci FST t t 1 because by chance the few individuals sampled were homogeneous within populations but dierent between populations, leading to high average FST. Heterozygosity can be similarly inated if the more frequent allele is likely to be absent (e.g., because mutations in the restriction site, which lead to allelic dropout, are oen in ancestral alleles that occur at a high frequency)18. Arnold, et al .13 conrmed results from Gautier, et al .18 and also concluded that other summary statistics, including and could be inaccurately estimated from loci with missing data. In spite of these problems, more recent simulation studies have indicated that missing data in RAD-Seq studies may not lead to incorrect inference, and in fact including loci with missing data can be advantageous for identifying shallow divergences22. Convention in phylogeographic studies oen is to require 75 or 80% of individuals to have data for a given locus—otherwise that locus is discarded from the analyses (e.g., refs23 – 28). Presumably, requiring a locus to be present in a certain number of individuals will eliminate loci with high missing data that may be the cause of misestimated parameters13,18. However, the choice of a cuto is arbitrary and is typically not justied in phylogeographic studies— the number of SNPs is virtually always reported as a single xed value (e.g., “we identied a total of 4,234 SNPs,” Jackson, et al .24). In reality, the various parameter values that determine how many loci are constructed and retained in SNP alignment methods means that there is a range of loci that could conceivably be included in a study27,29. To date, no phylogeographic study has investigated the effect of varying amounts of missing data in an empirical RAD-Seq dataset, even those explicitly comparing RAD-Seq-generated SNPs and microsatellites30, 31. To remedy this knowledge gap, we investigate the phylogeography of red mangroves ( Rhizophora mangle L., Rhizophoraceae) in Florida, using both an existing microsatellite dataset32, and new RAD-Seq SNP datasets that vary in number of loci and the percentage of missing data. We ltered RAD-Seq loci to generate a dataset that would approximate the number of loci and amount of missing data typically used in RAD-Seq phylogeography studies, and we also generated datasets with more or less stringent ltering to test the eects of increasing or decreasing the number of loci and percentage of missing data. Specically, we address the following questions: (1) In RAD-Seq datasets, how are phylogeographic inferences aected by the number of loci used? (2) In RAD-Seq datasets, how are phylogeographic inferences aected by the percentage of missing data? (3) What are the impor tant dierences in performance between microsatellites and RAD-Seq data in population genetic and phylogeographic inference? (4) Do RAD-Seq data reveal any novel phylogeographic inferences not already recovered by microsatellites for red mangroves in Florida?

PAGE 3

3 To address these questions, we used 96 red mangrove (Rhizophora mangle) individuals collected from 12 sampling locations on the coasts of Florida (Table 1 Fig. 1 ). Red mangroves are salt-tolerant trees that occur in coastal estuarine environments throughout the neotropics, experiencing high temperatures, frequent inundation, saline conditions, and periodic wave action associated with the coastal environment33. Red mangroves provide a variety of ecosystem services, including ltering water, providing habitat to animals, stabilizing shorelines, and protecting coastal environments from frequent wave action and occasional storm surges. us, red mangroves are important conservation targets—for which phylogeographic data can improve conservation strategies—making red mangroves a valuable study system. Sampling Location Code Latitude (N) Longitude (W) % Loci Missing Bahia Honda Key BHKFl 24.55286 81.76776 73.5 Convoy Point CvPFl 25.46347 80.33133 81.2 Cape Canaveral CpCFl 28.82173 80.75594 83.0 Hollywood HwdFl 26.03841 80.11780 79.4 Islamorada IsmFl 24.90031 80.65690 81.0 Key Largo KyLFl 25.09569 80.42957 88.9 Melbourne MlbFl 28.07435 80.60526 79.8 New Port Richey NPRFl 28.25432 82.75723 69.5 Seahorse Key ShKFl 29.10040 83.06185 65.8 Terra Ceia Bay TCBFl 27.59172 82.57524 81.7 Vaca Key VKyFl 24.71154 81.06992 85.1 West Palm Beach WPBFl 26.67505 80.04259 83.9Table 1. e twelve sampling locations (each containing eight individuals), their codes, GPS coordinates, and the percentage of loci that have missing data for each sampling location before any ltering. Figure 1. e 12 sampling locations (each with eight individuals) are indicated by orange circles. Sampling location codes are provided in Table 1. e map was generated using R (citation: R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/), and the R package ‘maps’ (citation: Original S code by Richard A. Becker, Allan R. Wilks. R version by Ray Brownrigg. Enhancements by omas P Minka and Alex Deckmyn. (2017). maps: Draw Geographical Maps. R package version 3.2.0. https://CRAN.R-project.org/packagemaps ).

PAGE 4

4 Further analysis of phylogeographic patterns in red mangroves and other species occurring in the Florida pen insula is also warranted. Although previous studies of many coastal and marine taxa revealed a phylogeographic discontinuity at or near the southern tip of Florida34,35, recent work on red mangroves using microsatellites failed to identify such a pattern32, 36. Dierent types of molecular markers could reveal new phylogeographic insights, due to broader sampling of the genome, and provide a predictive framework for understanding how genetic variation in this iconic species will respond to climate change. Finally, red mangroves are an ideal system for compar ing the performance of alternative genetic markers, given previous analyses of microsatellite loci32, 36, 37 and the size of the genome (approximately 1 t Gb; Hodel unpublished data, based on ow cytometry observations), enabling a rigorous test of the RAD-Seq method. Genome size is a necessary consideration with RAD-Seq; as genome size increases, the number of loci shared among many individuals for a given sequencing depth decreases.ResultsDatasets. Seven datasets, ranging from 8 loci (SSR_8) to 25,198 loci (RAD_25198; Table 2 ), were used to investigate in depth how variation in number of loci and percent missing data impacted phylogeographic infer ences. ese were selected from all possible datasets, for which basic statistics were calculated (Supplementary Figure 1). e name of each of the seven datasets contains information about locus type (RAD or SSR) and number of loci in the dataset. e smallest RAD dataset contained 239 loci (RAD_239), and the percentage of missing data for RAD datasets ranged from 11.7% to 78.1% (Table 3 ). e dataset RAD_1180, which required a locus to be present in 75 of 96 individuals (78.1%), most closely mimicked the amount of loci ltering typically used in a phylogeographic study. erefore, in our analyses, we used this as a baseline dataset against which to compare other RAD datasets. Across sampling locations, the proportion of missing data was relatively uniform (Table 1 ); percentage of missing loci in the data matrix for a given sampling location ranged from 65.8% to 88.9%.Population genetic analyses. Measures of heterozygosity were not signicantly dierent between the microsatellite dataset and the RAD datasets; average HO was 0.431 for the microsatellite dataset and 0.392 in RAD_1180, with a range from 0.354 to 0.477 across all RAD datasets (Table 3 ). Average HE was 0.388 for the microsatellite dataset and 0.307 for RAD_1180 and ranged from 0.300 to 0.340 for all RAD-Seq datasets (Table 3 ). Average FST for microsatellites was 0.124, which was signicantly greater than average FST for only one of the RAD datasets—the smallest (RAD_239; Table 3 ). Within the RAD datasets, average HO was signicantly greater in RAD_25198 than all others, and it was signicantly lower in RAD_6255 than in all others; HO did not predictably increase or decrease as the number of loci increased (Table 3 ). Additionally, within RAD datasets, average HE was signicantly greater in RAD_25198 than all others. FST ranged from 0.046 to 0.108 among the RAD datasets (Table 3). ere was no signicant dierence in FST in the three smallest RAD datasets, but the three largest RAD datasets all had increased FST relative to the smaller datasets (Table 3 ). e dataset with the largest value of FST was RAD_6255 (Table 3 ). Average FIS using microsatellites was not signicantly dierent than FIS calculated using RAD datasets; within RAD datasets, FIS generally increased as more loci were added, although RAD_25198 had the lowest value of FIS (Table 3). Many of the population genetic statistics were disproportionately aected by loci with very low or very high values of FST, FIS, or heterozygosity (Fig. 2). e eect of extreme loci was particularly evident in the larger datasets (RAD_6255 and RAD_25198), in which there were large numbers of loci with extreme values (e.g., FST of 1.0; Fig. 2).Pairwise FST. e values of pairwise FST for each sampling location relative to other sampling locations were remarkably consistent across datasets (Table 4). For most sampling locations pairwise FST estimated by SSRs was approximately twice as large as RAD dataset estimates. For every dataset, pairwise FST between Seahorse Key and all other sampling locations was the highest. For every RAD dataset, Islamorada had the lowest value for pairwise FST, but the SSR dataset identied West Palm Beach as the sampling location with the lowest estimate of pairwise FST. Even as the amount of missing data increased, the pairwise FST estimates remained consistent; RAD_25198 produced similar values to smaller RAD datasets (Table 4).FIS by sampling location. Cape Canaveral was the location with the highest FIS using the microsatellite data (SSR_8), followed by Key Largo and Seahorse Key (Table 5). Meanwhile, for all RAD datasets, Seahorse Key had one of the lowest (i.e., most negative) FIS values among all populations. Within the RAD datasets, the number of loci and/or amount of missing data aected FIS. For example, in Key Largo, the largest dataset yielded a value Dataset Individuals required to retain a locus Number of loci % individuals required to retain a locus RAD_239 83 239 86.5 RAD_1180 75 1180 78.1 RAD_2317 65 2317 67.7 RAD_3831 50 3831 52.1 RAD_6255 30 6255 31.3 RAD_25198 1 25198 1.0Table 2. e seven data sets used in this study; RAD-Seq data sets were generated by ltering loci from largest data set (RAD_25198). For all data sets (six RAD and one microsatellite), the total number of loci used is indicated.

PAGE 5

5 of 0.015, while the smallest dataset had a value of 0.194. is was not a large absolute change in FIS, but the interpretation of this statistic changed based on whether it is positive or negative (higher values indicate a greater level of inbreeding). In general, within RAD datasets, FIS increased as loci were added, although this trend was not universal, especially in the largest RAD dataset. For instance, in Bahia Honda Key, FIS was lowest in the largest dataset RAD_25198 (25,198 loci, 78.1% missing data). Conversely, in Islamorada, FIS was lowest in the smallest dataset (RAD_239, 11.7% missing data).Heterozygosity by sampling location. Observed heterozygosity for each sampling location ranged from 0.320 (Seahorse Key) to 0.451 (Hollywood) when averaged across all datasets (Table 6 ). For most datasets, Seahorse Key was the sampling location with the lowest HO, although notably RAD_25198 identied six other sampling locations with lower HO than Seahorse Key (Table 6 ). Similarly, most datasets reported Hollywood as the sampling location with the highest HO, but SSR_8 found Convoy Point and Islamorada had higher HO than Hollywood, and RAD_25198 identied ve other sampling locations with greater HO, with Key Largo having the highest HO (Table 6 ). For most sampling locations, measures of HO, when ranked relative to other sampling locations, remained similar across all RAD datasets except RAD_25198. Interestingly, the values of HO ranked relative to other sampling locations were more similar between SSR_8 and the ve smallest RAD datasets (RAD_239-RAD_6255) than any of the ve smallest RAD datasets were to RAD_25198 (Table 6).PCA and SVDQuartets. e PCA analysis revealed that microsatellite data did not identify clear groupings of individuals based on sampling location or other geographical divisions (Fig. 3 ). Similarly, RAD_239 did not dierentiate the samples into discrete clusters. However, RAD_1180, RAD_2317, RAD_3831, RAD_6255, and RAD_25198 all divided the samples into two groups with minimal overlap in the PCA visualization: one group was Gulf Coast samples, and the other group was Atlantic Coast samples (Fig. 3 ). Closer inspection of the PCAs revealed that most of the Cape Canaveral individuals formed a discrete cluster intermediate between the two other clusters (Gulf and Atlantic). Most RAD datasets had sucient resolution to place Cape Canaveral between the Gulf and Atlantic clusters, but the use of a small number of loci (i.e., RAD_239) was unable to show this relationship. Furthermore, the two largest datasets, RAD_6255 and RAD_25198, showed Cape Canaveral individuals clustering more closely to the Atlantic than the Gulf cluster. e 50% majority-rule consensus bootstrap trees generated with SVDQuartets showed substantial variation between datasets when inferring genealogical relationships between individuals and/or sampling locations (Fig. 4 ). In many cases, dataset RAD_239 did not identify genealogical relationships that were recovered with other datasets with more loci. However, certain key relationships among individuals were consistently shown in multiple datasets with thousands of loci. In every dataset except RAD_239 (i.e., every dataset with at least 1180 loci), all Seahorse Key (ShKFl) samples formed a clade (Fig. 4 ). In four datasets (RAD_2317, RAD_3831, RAD_6255, RAD_25198), all Gulf Coast (NPRFl, ShKFl, TCBFl) samples (except one individual: NPRFl_R8), together with all Cape Canaveral (CpCFl) samples, formed a clade that is sister to all remaining Atlantic Coast samples plus NPRFl_R8. Interestingly, this relationship was not recovered in RAD_1180, the dataset with ‘ideal’ ltering of loci—but all datasets with more loci (and therefore more missing data) did recover the relationship. Each RAD dataset had nodes with varying levels of bootstrap support (Fig. 4 ). Datasets with fewer loci showed few nodes with bootstrap support 70%; RAD_239 had three such nodes. More loci resulted in more nodes with bootstrap support 70%, up to a point: RAD_1180 had six highly supported nodes, RAD_2317 had eight, and RAD_3831 had the most with nine. en the number of highly supported nodes slightly declined as more loci were added: RAD_6255 and RAD_25198 each had eight nodes with bootstrap support 70% (Fig. 4).Sampling loci. Analyzing dierently sized samples of loci from RAD_25198 and SSR_8 provided several crucial insights. A microsatellite dataset with seven loci sampled from SSR_8 performed better in estimating FST than a dataset with six loci, although in each case, all 100 sampled replicates fell within the 95% condence interval of Figure 2. Stacked histograms of per locus estimates of FST, FIS, and HO for each of the RAD datasets. Datasets with more loci are stacked on top of datasets with fewer loci.

PAGE 6

6 FST for the complete SSR_8 dataset (Fig. 5). For all RAD datasets, the value of FST estimated using only originally ltered data is dierent from all 100 permuted values of FST calculated from an equivalent number of loci sampled from the largest dataset (RAD_25198). For almost all datasets, FST based on sampled loci was less than FST using original loci, except for one dataset (RAD_6255), FST based on the sampled loci was greater. Strikingly, in none of the datasets did the condence intervals from the sampled loci overlap with the condence intervals of the estimated FST values from the original data (Fig. 5 ). e percentage of missing data in the largest dataset clearly had an immense impact. Even when very few loci (e.g., 239 loci) were sampled from the largest dataset, the Figure 3. Principle component analysis (PCA) for all seven data sets. Note that the scales of the axes of the SSR_8 plot are dierent than the axes of all the RAD plots.

PAGE 7

7 distribution of FST values clustered around the estimated FST using all 25,198 loci (Fig. 5), indicating that missing data, not number of loci, aected the dierences in measured FST.DiscussionInsights regarding choice of loci. Our results indicate that ltering loci using the standard cuto (i.e., 75–80% of individuals must possess data for a given locus for that locus to be retained) should not be the gold standard in RAD-Seq studies—it is possible to retain many more loci without inated statistics23– 28. FST increased as missing data increased, as predicted by simulation studies, but the relationship is more nuanced than previously assumed. FST increases as the percentage of missing data increases—up to a point—and then decreases from RAD_6255 to RAD_25198, as the percentage of missing data nearly doubles, from 41.3% to 78.1% (Table 3 ). When more loci are included, the distribution of FST across the genome is more completely sampled. However, adding loci with more missing data may cause analyses to miss low-frequency alleles in the loci with extensive missing data, which would add error to average estimates of FST. Sampling analyses conrmed that FST generally Figure 4. Trees estimated using every individual for each RAD dataset in SVDQuartets. Orange branches indicate individuals from sampling locations in the Gulf of Mexico, and blue branches represent individuals from Atlantic sampling locations. Table 3. Relevant population genetic statistics for each of the seven data sets used in this study. For each column, warmer colors indicate lower values and cooler colors show higher values. Immediately to the right of each of the four columns (FST, FIS, HO, HE) is the 95% condence interval for each statistic. Dataset % Missing FSTFISHOHERAD_239 11.7 0.046 0.009 -0.365 0.107 0.410 0.028 0.300 0.013 RAD_1180 17.1 0.057 0.004 -0.298 0.061 0.398 0.016 0.307 0.007 RAD_2317 22.2 0.057 0.002 -0.253 0.033 0.390 0.010 0.312 0.004 RAD_3831 29.4 0.066 0.002 -0.213 0.027 0.382 0.006 0.315 0.003 RAD_6255 41.3 0.108 0.002 -0.164 0.022 0.356 0.005 0.306 0.002 RAD_25198 78.1 0.080 0.003 -0.403 0.016 0.477 0.003 0.340 0.002 SSR_8 0.0 0.124 0.067 -0.110 0.694 0.431 0.187 0.388 0.111

PAGE 8

8 Figure 5. Histograms showing the distribution of the 100 samplings of loci from a larger data set. In the rst two panels, six and seven SSR loci, respectively, were randomly sampled 100 times from the SSR_8 data set, and the distribution of the 100 calculations of FST are shown. e solid blue line indicates the parameter value estimated using all eight loci, and the dashed blue lines show the 95% condence interval. In the remaining ve plots, the histogram shows parameter estimates using the number of loci (239, 1,180, 2,317, 3,831, and 6,255, respectively) in the data set randomly sampled from RAD_25198 100 times. e solid blue lines indicate the FST value estimated using all 25,198 loci, and the dashed blue lines show the 95% condence interval. e solid orange lines indicate FST estimated using the original data set (RAD_239, RAD_1180, RAD_2317, RAD_3831, and RAD_6255, respectively) and dashed orange lines show the 95% condence interval for this estimate.

PAGE 9

9 increased as missing data increased (Fig. 5 ). Heterozygosity was less aected by missing data, as there was little or no change in either observed or expected heterozygosity when the percentage of missing data ranged between 11.7% (RAD_239) and 41.3% (RAD_6255). Only the largest dataset (RAD_25198, with 78.1% missing data) showed signicantly higher heterozygosity than other datasets. Some simulation studies reported that missing data could inate FST, and would likely inate estimates of heterozygosity, leading to calls for removing loci with incomplete sampling13. However, more recent simulation studies showed that removing loci with higher mutation rates, which are more likely to have missing data, negatively impacted phylogenetic analyses22. Our study shows the importance of thoroughly exploring how loci are ltered in empirical systems. Extreme amounts of missing data yield higher estimates of FST and heterozygosity and lower estimates of FIS (Table 3 ). A large number of loci in RAD_25198 had very high values of certain statistics (e.g., hundreds of loci with FST t t 0.975 and thousands of loci with HO t t 0.975), which severely impacted average estimates of these statistics (Fig. 2 ). Notably, not all datasets have these extreme loci—dataset RAD_3831, which only requires 52.1% of individuals to have data for a given locus, and had 29.4% missing data, did not suer from extreme loci, despite liberal ltering. Although missing data caused some statistics to increase, it did not dramatically aect our conclusions. For many analyses, using datasets other than RAD_1180, especially RAD_2317 and RAD_3831, did not change the interpretation of the results. Regardless of which of the three datasets was used, FST was relatively low—between 0.057 and 0.066. Importantly, nearly doubling the amount of missing data from 17.1% (RAD_1180, the ‘gold standard’) to 29.4% (RAD_3831) resulted in a very small increase in FST and no signicant change in other statistics (FIS, HO, HE; Table 3). Furthermore, using very few loci (RAD_239) did not signicantly change any of the statistics estimated using RAD_1180 (Table 3 ). Our data indicate that the oen-used cuto of 75–80% individuals with locus data is arbitrary, and dierent cutos should be considered and evaluated on a case-by-case basis to ensure an appropriate number of loci are used. e results suggest that in many cases, only minimal ltering of loci is needed, and many more loci can be retained than typically are. Researchers who wish to maximize the number of loci in their study could likely use very low cutos (e.g., require a locus to have data for 10% of individuals). Typically, microsatellite datasets have lower FST values relative to SNPs due to a larger number of alleles, although simulation studies have shown evidence that FST can be elevated up to an order of magnitude in microsatellite datasets due to factors such as correlated allele frequencies38. Average FST ranged from 0.046 to 0.124 across all datasets—so there is not high dierentiation detected in any dataset (Table 3 ). When using any RAD dataset except RAD_239, FST calculated using RAD loci was statistically indistinguishable from the microsatellite dataset (Table 3 ). In theory, a three-to-four-fold change in FST could alter biological conclusions—possibly with deleterious results (e.g., identifying populations or management units that would be prioritized for conservation)—but no matter how the loci were ltered, there was a relatively small range of estimated FST. RAD-Seq studies where larger values of FST were detected could exhibit larger absolute changes in FST when using dierent loci ltering cutos (e.g., refs39,40). Similarly, the interpretation of FIS and HO could impact how data are considered in a biological context (e.g., identifying locations at risk for inbreeding depression). Positive values of FIS and/or low values of HO oen indicate inbreeding, which means sampling locations are more vulnerable than other sampling locations. As noted earlier, SSR_8 identied ve sampling locations with positive FIS values (Table5 ). Only one of these sampling locations (Key Largo) also had positive FIS values in any of the RAD datasets. Clearly, marker choice (microsatellite versus RAD-Seq) aected the assessment of which populations are more vulnerable based on FIS values. Agreement between these two markers types would facilitate identifying sampling locations vulnerable to inbreeding depression. However, it is understandable that dierent markers would lead to dierent results, as Table 4. Pairwise FST for each sampling location (i.e., one sampling location versus all others) for each of the seven datasets. Within each data set, lower (warmer colors) and higher (cooler colors) values of FST are shown using color-coding. SSR_8 RAD_239 RAD_1180 RAD_2317 RAD_3831 RAD_6255 RAD_25198 Average BHKFl 0.101 0.039 0.055 0.049 0.052 0.068 0.063 0.061 CpCFl 0.155 0.055 0.069 0.066 0.069 0.075 0.078 0.081 CvPFl 0.085 0.040 0.048 0.046 0.050 0.061 0.062 0.056 HwdFl 0.163 0.052 0.056 0.059 0.063 0.073 0.076 0.078 IsmFl 0.106 0.030 0.040 0.039 0.044 0.051 0.050 0.052 KyLFl 0.108 0.063 0.071 0.073 0.081 0.100 0.097 0.085 MlbFl 0.134 0.043 0.052 0.052 0.055 0.068 0.074 0.068 NPRFl 0.130 0.045 0.060 0.062 0.064 0.083 0.081 0.075 ShKFl 0.279 0.117 0.134 0.138 0.143 0.158 0.152 0.160 TCBFl 0.100 0.056 0.077 0.076 0.083 0.101 0.098 0.084 VKyFl 0.135 0.045 0.050 0.052 0.056 0.067 0.069 0.068 WPBFl 0.083 0.046 0.058 0.056 0.053 0.070 0.073 0.063

PAGE 10

10 mutation rate can aect estimation of FIS—in theory, as the mutation rate increases, FIS should decrease. e data showed opposite result though, as FIS was higher in the SSR_8 dataset for most sampling locations (Table 5). is is likely because estimates of HO and HE have larger variance when few loci are used, as in the SSR_8 dataset (Table 3 ). e relative values of HO and HE can dramatically aect the interpretation of FIS, especially when HO and HE are similar (e.g., using the equation FIS t t ( HE t t HO)/ HE would yield FIS of 0.25 when HO t t 0.4 and HE t t 0.5, but FIS t t 0.25 if HO and HE are reversed). Within RAD datasets, estimates of FIS for each sampling location were fairly consistent and FIS increased as missing data increased, but this trend was not universal (Table 3 ). Identifying vulnerable sampling locations based on HO revealed that RAD_25198 led to dierent conclusions than most other datasets. Across datasets, measures of HO within each sampling location were consistent relative to other sampling locations, except for RAD_25198 (Table 6). Missing data impacted this analysis; a large number of loci in RAD_25198 had either very high or very low HO, possibly leading to the pattern of HO in RAD_25198 that contrasted with patterns in virtually every other dataset (Table 6, Fig. 2). e PCA results show that as the number of loci increases, the denition of clusters improves, plateauing with RAD_2317 or RAD_3831. e clustering is similar in all RAD datasets with 1,180 or more loci, with Cape Canaveral individuals falling between the Gulf and Atlantic clusters. As more loci are added, the Cape Canaveral samples appear to be closer to the Atlantic cluster, especially in datasets RAD_6255 and RAD_25198 (Fig. 3 ). Taking into account the SVDQuartets results claries the clustering—all Cape Canaveral samples form a clade with all Gulf samples except one. However, this relationship is only present in datasets with 2,317 loci or greater— the putatively ‘gold standard’ dataset RAD_1180 does not show this relationship.Phylogeographic patterns in red mangroves. Based on previous studies using microsatellite data32 36, the relationship of Cape Canaveral samples to other sampling locations, as found here with both PCA and SVDQuartet analyses of RAD-Seq data, was surprising—previous studies did not nd that any of the individuals in the Cape Canaveral population clustered with any of the Gulf samples. ese new data could indicate an Atlantic-Gulf phylogeographic discontinuity, and that Cape Canaveral is an anomaly due to a lack of phylogeographic resolution, recent population founding, or human-mediated transplantation. e intermediate placement of Cape Canaveral in many of the PCAs suggests that it may actually cluster with the Atlantic samples, especially when considering datasets RAD_6255 and RAD_25198, indicating a phylogeographic break (Fig. 3 ). However, the SVDQuartets results place Cape Canaveral in a clade with the vast majority of Gulf samples, although this relationship is not highly supported in any datasets (i.e., bootstrap support is not 70% for this clade in any dataset) (Fig. 4). Assuming that Cape Canaveral is more closely related to Gulf samples, the age of the divergence between the two clades (Atlantic, Gulf t t CpCFl) comes into question. Northern Florida represents the northern limit of the range of red mangroves33. Typically, populations of these trees in northern Florida are periodically extirpated due to freezing events, and these areas are re-colonized. e lower values of HO in northern populations (CpCFl, MlbFl, ShKFl, NPRFl, TCBFl) relative to southern populations indicate that these populations were likely founded more recently from a small number of propagules. e Cape Canaveral population was likely founded by individuals from the Gulf Coast, suggesting that the divergence between the two clades (Atlantic, Gulf t t CpCFl) is very recent. Previous research indicates that gene ow is greater from the Gulf Coast to the Atlantic Coast in red mangroves; there may be ongoing gene ow from the Gulf to Cape Canaveral32. Alternatively, alleles from the Gulf Coast could have migrated into an existing Cape Canaveral population and proliferated due to other processes (e.g., dri). Another explanation for the sister relationship between the Gulf samples and Cape Canaveral is Table 5. e variation in average inbreeding coecient (FIS) among data sets and populations. Within each data set, lower (warmer colors) and higher (cooler colors) values of FIS are shown using color-coding. e average value of FIS across all data sets for each population is shown in the last column of the table. SSR_ 8R AD_239 RAD_1180 RAD_231 7R AD_383 1 RAD_6255 RAD_25198 Average BHKFl 0.049 -0.207 -0.161 -0.132 -0.126 -0.118 -0.263 -0.140 CpCFl 0.095 -0.315 -0.208 -0.166 -0.159 -0.134 -0.113 -0.155 CvPFl -0.300 -0.238 -0.201 -0.165 -0.158 -0.126 -0.109 -0.191 Hw dFl -0.206 -0.356 -0.274 -0.231 -0.206 -0.147 -0.105 -0.228 IsmFl -0.130 -0.245 -0.163 -0.130 -0.111 -0.070 0.003 -0.128 KyLFl 0.081 -0.194 -0.168 -0.146 -0.133 -0.123 0.015 -0.109 MlbFl -0.018 -0.169 -0.121 -0.096 -0.092 -0.075 -0.039 -0.092 NPRFl -0.090 -0.317 -0.251 -0.225 -0.214 -0.209 -0.307 -0.234 ShKFl 0.067 -0.554 -0.506 -0.452 -0.426 -0.400 -0.427 -0.398 TCBFl -0.041 -0.407 -0.357 -0.311 -0.300 -0.277 -0.335 -0.299 VKyFl -0.109 -0.260 -0.214 -0.185 -0.158 -0.133 -0.086 -0.172 WPBFl 0.020 -0.233 -0.236 -0.218 -0.203 -0.180 -0.152 -0.177

PAGE 11

11 human-mediated transplantation of propagules or seedlings from the Gulf Coast to Cape Canaveral. However, all available publications and information from land managers who replied to requests for information conrm that any restoration that required importation of propagules used either local propagules or seedlings from the southern Atlantic Coast (ref.41, personal communication with Rangers from Cape Canaveral National Seashore). Another possible explanation for this result is that red mangrove propagules were accidentally transported from the Gulf Coast of Florida to Cape Canaveral during construction of the Kennedy Space Center in the 1960s. Construction of the Space Center was a massive project. It is noteworthy that nearly 100,000 tons of steel was transported from the Gulf Coast to Cape Canaveral in numerous trucks; the transport of a few mangrove propagules during this process could easily have established a Gulf genotype in the Cape Canaveral area42. We conclude that, in contrast to microsatellites, RAD datasets recover a relationship between the Gulf Coast and the Atlantic Coast (excluding Cape Canaveral) that supports the presence of a maritime discontinuity in red mangroves. However, as red mangroves can disperse long distances, a population or populations that recently established in Cape Canaveral likely had a founder or founders that were predominantly of Gulf Coast origin. e fact that previous studies using SSRs did not elucidate this relationship is not surprising—both the PCA analysis and SVDQuartets analysis indicate that 1180 loci were barely sucient to infer the placement of Cape Canaveral— datasets with many more loci were needed (Figs 3 and 4). e large number of loci required to resolve such relationships highlights why liberal ltering of RAD-Seq loci is advisable.ConclusionsWe cannot overemphasize the importance of thoroughly exploring RAD-Seq datasets when performing phylogeographic analyses—it is too easy to jump to conclusions when only using one arbitrary cuto to lter loci. Our empirical data conrm that estimates of FST and/or heterozygosity may become inated as missing data increase. However, this does not happen as quickly as implied in simulation studies as loci with missing data are added— liberal ltering of loci retains loci valuable for phylogeographic or phylogenetic inference, without inating population genetic statistics. us, regardless of the cuto value used to lter loci, researchers should investigate several other cutos with both increased and decreased amounts of missing data to appreciate fully the impact of missing data on parameters in their study. We found no evidence that the 75% or 80% cuto commonly employed was optimal. In many analyses, other datasets with cutos ranging from 31.3% to 67.7% performed just as well as or better than RAD_1180. Many RAD-Seq studies aim to multiplex as many individuals as possible in a HTS run; our results show that retaining loci with more missing data is feasible and advantageous in empirical studies, and that researchers can include more samples in a single sequencing run. Our study conrmed that microsatellites were a valuable tool for inexpensively estimating population genetic statistics, such as FST, FIS, and heterozygosity. However, this study revealed that the thousands of additional loci from across the genome provided by RAD-Seq increased phylogeographic resolution. We found that red mangroves likely have a phylogeographic discontinuity at the southern tip of Florida that was not detected in previous studies using SSRs32, 36 and that a single population from the Atlantic coast of Florida arose via recent colonization by propagules (either natural or human-mediated) from the Gulf coast.Methods and Materialsƒ’Ž‡…‘ŽŽ‡…–‹‘‹•‘Žƒ–‹‘ We collected leaf tissue from plants of R. mangle from 12 locations in Florida (Fig. 1 ). At each location, we collected one leaf from 10–20 individuals that were spaced at least 15 t m apart to minimize collecting closely related individuals. For each sampling location, we randomly selected 8 individuals Table 6. e variation in observed heterozygosity (HO) among data sets and populations. Within each data set, lower (warmer colors) and higher (cooler colors) values of HO are shown using color-coding. e average value of HO across all data sets for each population is shown on the bottom row of the table. SSR_ 8R AD_239 RAD_1180 RAD_2317 RAD_3831 RAD_6255 RAD_25198 Average BHKFl 0.391 0.429 0.419 0.414 0.403 0.377 0.419 0.408 CpCFl 0.359 0.379 0.346 0.350 0.355 0.336 0.306 0.348 CvPFl 0.656 0.431 0.418 0.413 0.397 0.363 0.305 0.425 Hw dFl 0.516 0.487 0.464 0.451 0.435 0.395 0.392 0.451 IsmFl 0.531 0.443 0.414 0.406 0.392 0.364 0.387 0.420 KyLFl 0.438 0.393 0.370 0.361 0.339 0.308 0.468 0.382 MlbFl 0.328 0.412 0.395 0.394 0.386 0.360 0.311 0.372 NPRFl 0.359 0.353 0.339 0.341 0.348 0.335 0.403 0.351 ShKFl 0.313 0.318 0.306 0.305 0.317 0.312 0.394 0.320 TCBFl 0.453 0.366 0.354 0.360 0.364 0.351 0.415 0.376 VKyFl 0.422 0.449 0.431 0.424 0.408 0.371 0.332 0.408 WPBFl 0.406 0.459 0.444 0.438 0.426 0.389 0.346 0.418

PAGE 12

12 to use in genetic analyses. GPS coordinates for each sampling location were recorded (Table 1 ). Each sampled leaf was placed in a labeled bag with silica gel and stored for 1–12 months; we then extracted DNA from the dried leaf tissue using a standard CTAB protocol43.‹…”‘•ƒ–‡ŽŽ‹–‡ƒ’Ž‹¤…ƒ–‹‘ƒ†ƒƒŽ›•‹• We PCR-amplied eight nuclear microsatellite loci for R. mangle (RM 11, 19, 21, 36, 38, 41, 46, 47)37. An M13 protocol44 was used to label amplicons with four uorescent dyes (6-FAM, NED, PET, VIC). e PCR (25L reactions) contained: 5X buer (5 t L), 2.5 t mM MgCl2 (2 t L), 2.5 t mM dNTP (0.5 t L), 0.12 t M forward primer with M13 label attached (1.25 t L), 4.5 t M reverse primer (1.25 t L), 4.5 t M uorescent dye (2.5 t L), H2O (10 t L), Taq polymerase (0.5 t L), and 50 ng template DNA (2 t L). We carried out PCR in a Biometra T3 ermocycler (Whatman Biometra, Goettingen, Germany) using the following cycles: initial denaturing at 94 t C for 3 t minutes; 35 cycles of 94 t C (45 t seconds), 52 t C (45 t seconds), 72 t C (75 t seconds); nal elongation at 72 t C for 20 t minutes. We used the Applied Biosystems 3730 DNA Analyzer (Applied Biosystems, Foster City, United States) at the University of Florida Interdisciplinary Center for Biotechnology Research to detect the uorescent peaks. We determined microsatellite peaks in Geneious 6.5 (http://www.geneious.com/) using the GeneScan 600 size standard ladder for calibration45.RAD-Seq library preparation and data processing. We followed the double-digest RAD-Seq protocol developed by Peterson, et al .46. For each sample, we constructed 96 DNA libraries by digesting approximately 200 ng genomic DNA with Eco RI and MseI. We then ligated Illumina adapters and unique 8–10-nucleotide barcodes to the DNA fragments. e DNA libraries were PCR-amplied in two separate reactions and pooled to minimize early PCR bias. We size selected 250–450-bp fragments using gel electrophoresis and sequenced the DNA fragments using the 1 t t 100-bp setting on the Illumina HiSeq. 2500 platform. Raw sequence data were deposited in the NCBI Sequence Read Archive (accession numbers pending). We processed the raw Illumina reads using the FAST-X toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) to lter sequences; we required 95% of bases to be above a quality score of 30 to retain a read. We then converted the sequences from FASTQ to FASTA, demultiplexed the reads, sorted them by barcodes, and trimmed the sequences by removing the nal 2 bases to ensure that we were using only high-quality sequence data. We assembled the sequences into loci using the STACKS 1.24 pipeline47 with the following parameter settings: -n 3 -m 3 -M 2 (parameters were optimized following Mastretta-Yanes, et al .29); all other parameters were le as the default. We selected seven datasets (one microsatellite and six RAD-Seq) and used a variety of analyses to compare the results produced by each dataset (Table 2 ). We used the ‘populations’ program in STACKS to produce an unltered dataset of RAD-Seq loci using the ‘write single SNP’ command and requiring a minor allele frequency 0.05. We then removed human, fungal, and microbial contamination from the loci and ltered loci by representation across individuals using an R script to create ve smaller datasets (Data_aquisition.R; this script and all other scripts are available at https://github.com/richiehodel/ red_mangrove_RAD_SSR ). Filtered datasets were required to have locus data for a certain number of individuals for the given locus to be retained in the analysis; the number of individuals could range from 1–96 (Supplemental Fig.1). e datasets were chosen such that they encompassed a wide range of loci and missing data.Population genetic analyses. We used an R script (Basic_stats.R) and the R package ‘hierfstat’48 to cal culate average FST, the inbreeding coecient FIS, HO, and HE for each of the seven datasets. To investigate how the number of loci aected comparisons of population genetic statistics among populations, we calculated pair wise FST (one sampling location versus all others combined) for each sampling location for each dataset using GenoDive49 and an R script (Pairwise_Fst.R). Additionally, we calculated FIS and HO for each sampling location for each dataset to determine how measures that oen inform conservation practices might be aected by the number of loci and amount of missing data. We measured how missing data were partitioned across sampling locations to verify that there were not any sampling locations with unusually high or low amounts of missing data (Table 1 ). Additionally, we investigated how several population genetic statistics were distributed across loci in each of the datasets (Stat_Distribution.R; Fig. 2).Principal components and SVDQuartets. We used a principal component analysis (PCA) implemented in the R package ‘SNPRelate’50 to identify clusters of individuals in the RAD data with an R script (VCF_PCA.R) and GenoDive to run a PCA on the microsatellite data. Aer visualizing the initial results, we tested several ways of grouping sampling locations together based on geography. We used SVDQuartets51 to determine genealogical relationships among individuals. is program selects the optimal topology for a quartet of taxa, and, aer sampling millions of quartets, infers a phylogeny for all individuals based on choosing the quartets with the best scores and assembling them into a phylogenetic tree. We used an R script (Nexus_creation.R) to convert the output from the ‘populations’ program in STACKS into nexus les that could be read for the SVDQuartets analysis. For each RAD dataset, we evaluated all possible quartets and selected trees under the multispecies coalescent using QFM (Quartet Fiduccia Mattheyses) quar tet assembly52. We used non-parametric bootstrapping (100 replicates for each dataset) to assess condence in inferred genealogical relationships between individuals. e R script Tree_formatting.R was used to visualize and annotate the 50% majority-rule trees from SVDQuartets using the R packages ‘ape’53 and ‘ggtree’54.Sampling loci. To test whether the number of loci or percentage of missing data for the loci used is the more important factor impacting measures of xation, population dierentiation, and heterozygosity, we randomly sampled from RAD_25198 (the RAD-Seq dataset comprising 25,198 loci) the equivalent number of loci contained in RAD_239, RAD_1180, RAD_2317, RAD_3831, and RAD_6255, respectively, and used these ve sets of sampled loci in analyses. We used an R script (Subsample.R) to randomly sample loci without replacement from RAD_25198 and repeated the sampling 100 times for each dataset. We compared measures of FST calculated using the original datasets with results

PAGE 13

13 calculated using the sampled loci from RAD_25198 (Fig. 5 ). We used bootstrapping to calculate 95% condence inter vals around FST for the original datasets and for the sets of loci sampled from RAD_25198 (Fig. 5 ).Data availability. e datasets generated during the current study are available in the NCBI Genbank repository, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA397667 (accession numbers SRR5918296-SRR5918355).Referencest 1.t Guichoux, E. et al Current trends in microsatellite genotyping. Mol. Ecol. esour. 11, 591–611 (2011).t 2.t alia, ., ai, M. ., alia, S., Singh, & Dhawan, A. Microsatellite marers: an overview of the recent progress in plants. Euphytica 177, 309–334 (2010).t 3.t Hodel, G. J. et al e report of my death was an exaggeration: A review for researchers using microsatellites in the 21st century. Appl. Plant Sci. 4, (2016).t 4.t Seeb, J. E. et al Single-nucleotide polymorphism (SNP) discovery and applications of SNP genotyping in nonmodel organisms. Mol. Ecol. esour. 11, 1–8 (2011).t 5.t Gardner, M. G., Fitch, A. J., Bertozzi, T. & Lowe, A. J. ise of the machines–recommendations for ecologists when using next generation sequencing for microsatellite development. Mol. Ecol. esour. 11, 1093–101 (2011).t 6.t Hodel, G. J. et al A new resource for the development of SS marers: Millions of loci from a thousand plant transcriptomes. Appl. Plant Sci. 4, (2016).t 7.t Andrews, ., Good, J. M., Miller, M. ., Luiart, G. & Hohenlohe, P. A. Harnessing the power of ADseq for ecological and evolutionary genomics. Nat. ev. Genet. 17, 81–92 (2016).t 8.t Baird, N. A. et al apid SNP discovery and genetic mapping using sequenced AD marers. PLoS One 3, e3376 (2008).t 9.t Smith, A. M. et al Highly-multiplexed barcode sequencing: an ecient method for parallel analysis of pooled samples. Nucleic Acids es. 38, e142–e142 (2010).t 10.t Andolfatto, P. et al Multiplexed shotgun genotyping for rapid and ecient genetic mapping. Genome es. 21, 610–617 (2011).t 11.t Sonah, H. et al An improved genotyping by sequencing (GBS) approach oering increased versatility and eciency of SNP discovery and genotyping. PLoS One 8, e54603 (2013).t 12.t afalsi, A. Applications of single nucleotide polymorphisms in crop genetics. Curr. Opin. Plant Biol. 5, 94–100 (2002).t 13.t Arnold, B., Corbett-Detig, B., Hartl, D. & Bomblies, ADseq underestimates diversity and introduces genealogical biases due to nonrandom haplotype sampling. Mol. Ecol. 22, 3179–3190 (2013).t 14.t Etter, P. D., Bassham, S., Hohenlohe, P. A., Johnson, E. A. & Creso, W. A. SNP discovery and genotyping for evolutionary genetics using AD sequencing in Methods in Molecular Biology 157–178 (Clion, 2012).t 15.t Xu, P. et al Population genomic analyses from low-coverage AD-Seq data: a case study on the non-model cucurbit bottle gourd. Plant J. 77, 430–442 (2014).t 16.t Lowry, D. B. et al Breaing AD: an evaluation of the utility of restriction site-associated DNA sequencing for genome scans of adaptation. Mol. Ecol. esour. 17, 142–152 (2017).t 17.t Davey, J. W. et al Special features of AD sequencing data: implications for genotyping. Mol. Ecol. 22, 3151–64 (2013).t 18.t Gautier, M. et al e eect of AD allele dropout on the estimation of genetic variation within and between populations. Mol. Ecol. 22, 3165–78 (2013).t 19.t Liu, N., Chen, L., Wang, S., Oh, C. & Zhao, H. Comparison of single-nucleotide polymorphisms and microsatellites in inference of population structure. BMC Genet. 6(Suppl 1), S26 (2005).t 20.t Coates, B. S. et al Comparative performance of single nucleotide polymorphism and microsatellite marers for population genetic analysis. J. Hered. 100, 556–564 (2009).t 21.t Schopen, G. C. B., Bovenhuis, H., Viser, M. H. P. W. & van Arendon, J. A. M. Comparison of information content for microsatellites and SNPs in poultry and cattle. Anim. Genet. 39, 451–453 (2008).t 22.t Huang, H. & nowles, L. L. Unforeseen consequences of excluding missing data from Next-Generation Sequences: Simulation study of AD sequences. Syst. Biol. 65, 357–65 (2016).t 23.t Catchen, J. et al e population structure and recent colonization history of Oregon threespine sticlebac determined using restriction-site associated DNA-sequencing. Mol. Ecol. 22, 2864–2883 (2013).t 24.t Jacson, A. M. et al Population structure and phylogeography in Nassau grouper (Epinephelus striatus ), a mass-aggregating marine sh. PLoS One 9, e97508 (2014).t 25.t Bernardi, G., Azzurro, E., Golani, D. & Miller, M. Genomic signatures of rapid adaptive evolution in the bluespotted cornetsh, a Mediterranean Lessepsian invader. Mol. Ecol. 25, 3384–3396 (2016).t 26.t Blanco-Bercial, L. & Buclin, A. New view of population genetics of zooplanton: AD-seq analysis reveals population structure of the North Atlantic plantonic copepod Centropages typicus. Mol. Ecol. 25, 1566–1580 (2016).t 27.t odrguez-Ezpeleta, N. et al Population structure of Atlantic macerel inferred from AD-seq-derived SNP marers: eects of sequence clustering parameters and hierarchical SNP selection. Mol. Ecol. esour. 16, 991–1001 (2016).t 28.t Van Wyngaarden, M. et al Identifying patterns of dispersal, connectivity and selection in the sea scallop, Placopecten magellanicus using ADseq-derived SNPs. Evol. Appl. 10, 102–117 (2017).t 29.t Mastretta-Yanes, A. et al estriction site-associated DNA sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference. Mol. Ecol. esour. 15, 28–41 (2015).t 30.t Bradbury, I. et al Transatlantic secondary contact in Atlantic Salmon, comparing microsatellites, a single nucleotide polymorphism array and restriction-site associated DNA sequencing for the resolution of complex spatial structure. Mol. Ecol. 24, 5130–5144 (2015).t 31.t Jeries, D. L. et al Comparing ADseq and microsatellites to infer complex phylogeographic patterns, an empirical perspective in the Crucian carp, Carassius carassius L. Mol. Ecol. 25, 2997–3018 (2016).t 32.t Hodel, G. J., Cortez, M. B., de, S., Soltis, P. S. & Soltis, D. E. Comparative phylogeography of blac mangroves (Avicennia germinans) and red mangroves (hizophora mangle) in Florida: Testing the maritime discontinuity in coastal plants. Am. J. Bot. 103, 730–739 (2016).t 33.t Tomlinson, P. B. e Botany of Mangroves (Cambridge University Press, 2016).t 34.t Avise, J. C. Phylogeography: e History and Formation of Species (Harvard University Press, 2000).t 35.t Soltis, D., Morris, A., McLachlan, J. S., Manos, P. S. & Soltis, P. S. Comparative phylogeography of unglaciated eastern North America. Mol. Ecol. 15, 4261–4293 (2006).t 36.t ennedy, J. P. et al Contrasting genetic eects of red mangrove ( hizophora mangle L.) range expansion along West and East Florida. J. Biogeogr. 44, 335–347 (2017).t 37.t Fu, ., Dey, D. & Holsinger, E. Bayesian models for the analysis of genetic structure when populations are correlated. Bioinformatics 21, 1516–1529 (2005).t 38.t osero-Galindo, C., Gaitan-Solis, E., Crdenas-Henao, H., Tohme, J. & Toro-Perea, N. Polymorphic microsatellites in a mangrove species, hizophora mangle L.(hizophoraceae). Mol. Ecol. Notes 2, 281–283 (2002).t 39.t ang, J., Ma, X. & He, S. Population genetics analysis of the Nujiang catsh Creteuchiloglanis macropterus through a genome-wide single nucleotide polymorphisms resource generated by AD-Seq. Sci. ep. 7, 2813 (2017).

PAGE 14

14 t 40.t Manthey, J. D., Geiger, M. & Moyle, G. elationships of morphological groups in the northern icer superspecies complex (Colaptes auratus & C. chrysoides). Syst. Biodivers. 15, 183–191 (2017).t 41.t Johnson, L.. & Herren, L.W. e-establishment of fringing mangrove habitat in the Indian iver Lagoon 19–22 (Florida Department of Environmental Protection, 2008).t 42.t NASA Public Aairs. e ennedy Space Center Story (Graphic House, 1991).t 43.t Doyle, J. J. & Doyle, J. L. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 19, 11–15 (1987).t 44.t Schuele, M. An economic method for the uorescent labeling of PC fragments. Nat. Biotechnol. 18, 233–234 (2000).t 45.t earse, M. et al Geneious Basic: an integrated and extendable destop soware platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–9 (2012).t 46.t Peterson, B. ., Weber, J. N., ay, E. H., Fisher, H. S. & Hoestra, H. E. Double Digest ADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One 7, (2012).t 47.t Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Creso, W. A. Stacs: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–40 (2013).t 48.t Goudet, J. hierfstat, a pacage for to compute and test hierarchical F-statistics. Mol. Ecol. Notes 5, 184–186 (2005).t 49.t Meirmans, P. G. & Van Tienderen, P. H. Genotype and Genodive: two programs for the analysis of genetic diversity of asexual organisms. Mol. Ecol. Notes 4, 792–794 (2004).t 50.t Zheng, X. et al A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28, 3326–3328 (2012).t 51.t Chifman, J. & ubato, L. Quartet inference from SNP data under the coalescent model. Bioinformatics 30, 3317–3324 (2014).t 52.t eaz, et al Accurate phylogenetic tree reconstruction from quartets: A heuristic approach. PLoS One 9, e104008 (2014).t 53.t Paradis, E., Claude, J. & Strimmer, APE: analyses of phylogenetics and evolution in language. Bioinformatics 20, (289–290 (2004).t 54.t Yu, G., Smith, D. ., Zhu, H., Guan, Y. & Lam, T. T. Y. ggtree: an pacage for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution 8, 28–36 (2017).AcknowledgementsWe thank the following funding sources: Botanical Society of America, Florida SeaGrant, Sigma Xi, and the University of Florida Department of Biology.Publication of this article was funded in part by the University of Florida Open Access Publishing Fund.Author ContributionsR.G.J.H., S.C. and A.C.P. collected the data, R.G.J.H., S.F.M., P.S.S. and D.E.S. contributed nancially to the data collection, all authors designed analyses, R.G.J.H. wrote the manuscript, and all authors reviewed and edited the manuscript.Additional InformationSupplementary information accompanies this paper at https://doi.org/10.1038/s41598-017-16810-7. Competing Interests : e authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional aliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. e images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. e Author(s) 2017