Island biogeography and ecological modeling of the amblypygid Phrynus marginemaculatus in the Florida Keys archipelago

Abstract Aim The biogeography of terrestrial organisms across the Florida Keys archipelago is poorly understood. We used population genetics and spatioecological modeling of the Amblypygi Phrynus marginemaculatus to understand the genetic structure and metapopulation dynamics of Keys populations that are otherwise isolated by human development and ocean. Location The Florida Keys archipelago and mainland Florida. Methods We sequenced a 1,238 bp fragment of mtDNA for 103 individuals of P. marginemaculatus from 13 sites in the Florida Keys and South Florida, binned into four regions. We used population genetic analyses to understand the population structure of the species throughout its US range. Furthermore, we used ecological modeling with climate, habitat, and human development data to develop habitat suitability estimates for the species. Results We found clear genetic structure between localities. The Lower Keys, in particular, support populations separate from those in other regions studied. Ecological modeling and genetic analyses showed the highest habitat suitability and genetic isolation in the Lower Keys, but urban development across the species range has resulted in the loss of most historical habitat. Main conclusions A mainland‐metapopulation model best fits P. marginemaculatus gene flow patterns in the Florida Keys and mainland. Ocean currents likely play a role in metapopulation dynamics and gene flow for terrestrial Keys species like P. marginemaculatus, and genetic patterns also matched patterns consistent with geologic history. Suitable habitat, however, is limited and under threat of human destruction. The few remaining pockets of the most suitable habitat tend to occur in parks and protected areas. We argue that conservation efforts for this species and others in the terrestrial Florida Keys would benefit from a deeper understanding of the population genetic structure and ecology of the archipelago.


| INTRODUC TI ON
A central goal of biology is to understand how time and space shaped the evolutionary history of life (MacArthur & Wilson, 1967;Podas, Crisci, & Katinas, 2006;Warren et al., 2015). The Florida Keys is a biodiverse archipelago with high endemism (Forys & Allen, 2005;Kautz & Cox, 2001), but the biogeography of terrestrial populations along the island chain and nearby mainland Florida remains poorly understood. This is particularly surprising considering the Florida Keys was the study site used to develop some of the most iconic biogeography theory (e.g., MacArthur & Wilson, 1967;Simberloff, 1976;Simberloff & Wilson, 1969, Simberloff & Wilson 1970Wilson & Simberloff, 1969) .
The Florida Keys may support metapopulation spatial structure for some species or genetic divergence and speciation in others (Hanski, 1998;Shrestha, Wirshing, & Harasewych, 2015). In either case, the distribution of species across the Keys is important for understanding both species long-term survival and biodiversity. This is especially true considering the precariousness of Florida Keys habitats in the face of human disturbance, including human development, deforestation, nonnative species introductions, and human-induced climate change impacts like sea level rise, increased catastrophic storms, and altered fire regimes (Bancroft, Strong, & Carrington, 1995;Forys, 2005;Maschinski et al., 2011;Ross, O'Brien, & da Silveira Lobo Sternberg, 1994;Ross, O'Brien, Ford, Zhang, & Morkill, 2009). Furthermore, the Florida Keys remain a major tourist destination with over 4.5 million tourists visiting annually (McClenachan, 2013).
Our understanding of the genetic structure of Florida Keys organisms largely comes from research of marine species, where ocean currents play a major role in gene flow and migration (Apodaca, Trexler, Jue, Schrader, & Travis, 2013;DeBiasse, 2010;Kirk, Andras, Harvell, Santos, & Coffroth, 2009;Lacson & Morizot, 1991) . Genetic patterns of terrestrial species are expected to differ considerably, as the marine ecosystem acts as an uninhabitable matrix and the formation of terrestrial habitat occurred on different timescales (Hoffmeister & Multer, 1968;Shrestha et al., 2015).
Genetic research using native terrestrial species is scarce, but not absent. The mosquito Aedes aegypti showed no genetic structure along the Florida Keys, likely because they disperse via flight (Brown, Obas, Morley, & Powell, 2013). The invasive brown anole (Anolis sagrei) and greenhouse frog (Eleutherodactylus planirostris) showed introductions to the Florida Keys from Cuba but cannot inform patterns for native species across the Keys (Heinicke, Diaz, & Hedges, 2011;Kolbe et al., 2004). An allozyme electrophoretic study on the Florida Tree Snail (Liguus fasciatus) revealed low levels of genetic diversity, likely due to a recent introduction from Cuba, or the low resolution of allozyme approaches (Hillis, Dixon, & Jones, 1991). The land snail Cerion incanum showed some haplotype structure between Upper and Lower Keys, likely caused by differences in the timing of formation of the Keys. Shrestha et al. (2015) proposed that the C. incanum spread southwesterly to colonize new Keys as they formed, with Lower Key populations being the youngest. Lastly, ant gut microbiota showed genetic structure between the upper and lower keys (Hu et al. 2013 (Peterson, Lopez, Frank, Porter, & Silvy, 2004).
We aimed to quantify the genetic and ecological characteristics of P. marginemaculatus populations in the US. In particular, we aimed to uncover the genetic structure of the species across the Florida Keys archipelago, understand the evolutionary history of Keys populations in relation to the species range via phylogenetic analysis, and identify suitable habitat and locations of putative populations throughout the species' potential range. Together, these results will be the first to examine P. marginemaculatus in the wild, and will provide critical information applicable to many terrestrial island species and their conservation.

| Study site
The Florida Keys is a ca 250-km-long archipelago amounting to ca 350 km 2 of dry land, extending from Key Largo bordering mainland Florida southwest to Key West 140 km from Cuba. The Keys are made from two geologic formations that both formed during the Tarantian Pleistocene (0.126-0.0117 mya) and raising above sea level during the Wisconsin glaciation (ca. 100,000 before present; Hoffmeister & Multer, 1968;Shrestha et al., 2015). The Lower Keys (Big Pine Key to Key West) formed from cemented sand bars, resulting in oölitic limestone (termed Miami Limestone). These Keys curve west and orient laterally due to gulf stream currents. The Upper Keys (Bahia Honda Key to Key Largo), however, are constituted of fossil coral reefs (termed Key Largo Limestone) without strong lateralization. Some species show endemism to the Lower or Upper Keys, but not both because of these differences (e.g., Peck & Howden, 1985).
Florida mixed hardwood forests have persisted since the birth of terrestrial Florida 25 mya (Webb, 1990). Since then, species have had two routes to colonizing Florida habitats: a land migration down the Florida peninsula from North America, or water migration north or west from the Caribbean and Bahamas, as Florida was never connected to the Caribbean islands (Snyder, 1990). Thus, we might expect the diversity of Florida to be shaped by tropical species able to disperse over water, and temperate species only able to establish via land. This has resulted in a dominance of vertebrates from North America but flora from the Caribbean (Snyder, 1990). Exceptions include nine bird species, and two species each of bat, frog and lizard, all of which have Caribbean origins. That being said, natural migrations are only clear for a few of these species; many might have been introduced via humans, and still others have gone extinct (at least locally) since their discovery (Snyder, 1990). Invertebrate biogeographic patterns are more mixed, but still fit the model of tropical water migrators versus temperate land migrators. For example, most Florida ant species have North American origins, while Butterflies are largely Caribbean (Lenczewski, 1980). Land to support the growing human populations of Miami and the Florida Keys have been largely obtained by clearing pineland and hammock (Snyder, Herndon, & Robertson, 1990), a major conservation issue that researchers have been bringing attention to for nearly a century (Small, 1929). The first settlers in Southern Florida were concentrated in the Florida Keys. Early settlers exploited pine and hardwood trees (especially Mahogany) for lumber, fuel, and slashand-burn agricultural practices (Small, 1917;Browder, Littlejohn, & Young, 1976;Wilson & Porras, 1983) . As a result, very few stands of rockland include original forest. Industrial logging was enabled by the Florida East Coast Railroad, which reached Miami in 1896.
Rockland habitat was subject to clear-cutting for timber and fuel but made poor agricultural land due to an abundance of limestone rocks that made soil unworkable. Greatly expanding agriculture oftentimes spared rockland habitat in favor of draining glades for crops until the invention of the rock plow in the 1950s. This enabled limestone rocks to be collected and separated from soil, thereby enabling agricultural access. Limestone rocks collected via hand, plow, and mine, were used as building materials and can be seen in historical buildings, walls, and gardens of the Florida Keys today. Rocklands are nutrient and water depauperate, and thus require the heavy use of irrigation and chemical fertilizers to be agriculturally usable. As a result, abandoned rocklands show little resemblance to the original ecosystem, and are often dominated by invasive species (Loope & Dunevitz, 1981). The ability to turn rockland into space for agriculture and housing has led to the steep decline in rockland habitats that continue to the present, with practically no hope of reestablishment without human intervention (Dorn, 1956;Meyers & Ewel, 1990;Possley, Maschinski, Maguire, & Guerra, 2014;Snyder, 1990).

| Study species
We used the amblypygid species P. marginemaculatus C.L. Koch, 1840 as a model to understand the general biogeographic pattern of terrestrial Florida Keys species. Amblypygids are a small arachnid order (ca. 220 spp.) of large nocturnal predators (Chapin & Hebets, 2016). Phrynus marginemaculatus is the only amblypygid species in the US east of the Mississippi river and the most commonly studied species of amblypygid (Chapin & Hebets, 2016; Figure 1). Laboratory research has shown that P. marginemaculatus exhibit ritualized agonistic displays (Fowler-Finn & Hebets, 2006), and can learn to navigate mazes using tactile cues (Santer & Hebets, 2009a, 2009b via exceptional brain structures and sensory systems (Chapin & Hebets, 2016;Santer & Hebets, 2011). While fascinating laboratory research has been conducted on the species, no research on their habitat requirements, distribution, population ecology, or population genetics has ever been published (Chapin & Hebets, 2016;Weygoldt, 2000).
Historical records indicate that the species was found as far north as Martin County, FL, but recent records are absent. The species is also found on several Bahamian islands, Cuba, Jamaica, and Hispaniola Quintero, 1981). Research on the species, however, has only occurred in captivity, with animals collected from a single island (Big Pine Key; Hebets & Chapman, 2000;Fowler-Finn & Hebets, 2006;Spence & Hebets 2006;Santer & Hebets, 2009a, 2009b , the pet trade (Rayor & Taylor, 2006), or both (Graving, 2015).
Interestingly, P. marginemaculatus has evolved a plastron to breathe underwater, which they can do for upwards of 24 hr (Hebets & Chapman, 2000). While the function of the plastron is not well understood, it likely increases chances of survival during flooding in their terrestrial retreats. This may be particularly important in the Florida Keys, where annual hurricanes can result in flooding.
Furthermore, hurricanes, along with ocean currents, likely promote oceanic dispersal (Fleming & Murray, 2009;Gillespie et al., 2011) . A plastron allowing underwater breathing likely extends the dispersal propensities across bodies of water, and the likelihood of survival during ocean migration, in P. marginemaculatus.
Two molecular phylogenetic studies have focused on amblypygids. First, a phylogeny of the Damon variergatus group delineated two cryptic species within Damon, an African genus of amblypygid (Prendini, Weygoldt, & Wheeler, 2005). Second, phylogenetic analyses of Phrynus species in Puerto Rico revealed hidden dimensions of diversity across cave populations (Esposito et al., 2015). In particular, Esposito et al. (2015) noted high levels of diversity across mitochondrial but not nuclear, genomes. Here we focus on mitochondrial sequences to examine if similar genetic structure across small geographical scales is evident in the Florida Keys.

| Specimen collection
We collected P. marginemaculatus specimens from 13 locations in southern mainland Florida and the Florida Keys archipelago using a nonrandom sampling method (Table 1; Figure 2). We limited our survey to upland habitat types, as we assumed that P. marginemaculatus would not be found in intertidal habitats like mangrove swamps and floodplains. Additional survey areas were selected from historic records, which were all associated with these upland habitat types (Table 1). P. marginemaculatus hide under debris, especially limestone rocks, during the day (Chapin & Hebets, 2016;Hebets & Chapman, 2000). Our sampling regime included walking trails and looking for P. marginemaculatus under rocks, logs, and other larger debris. When found, we collected genetic samples and stored them in 95% ethanol on dry ice. We dissected muscle tissue from one or a few appendages, depending on the size of the specimen.

| Extraction, amplification, and sequencing
We extracted genomic DNA using QIAGEN DNeasy Blood & Tissue Kits. We followed the standard kit protocol but used chilled ethanol (−20°C) and a 50-mL final elution. We amplified a 1,238 nucleotide sequence of the mitochondrial gene cytochrome c oxidase subunit 1 (COI) by performing 34 iterations of the following cycle on a thermal cycler: 30 s at 94°C, 35 s at 48°C, and 90 s at 65°C, beginning with an initial cycle of 2 min at 94°C and ending with 10 min at 72°C. Using illustra PuReTaq Ready-To-Go PCR beads and 400-nM forward and reverse primers, the long or short read of COI was sequenced for each sample. LCOI1490 was used as the forward primer for both the long and short reads, while HCOI2198 was used as the reverse primer for the short reads and C1-N-2776 was used as the reverse primer for the long reads (LCOI1490 GGTCAACAAATCATAAAGATATTGG, HCOI2198 TAAACTTCAGGGTGACCAAAAAATCA, andC1-N-2776 GGATAATCAGAATATCGTCGAGG;Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994). We chose a mtDNA sequence as mtDNA is much more informative than nuclear DNA among Amblypygi (Esposito et al., 2015).

| Ecological modeling
We used locality data from the 103 field-collected specimens to es-  (Miller & Glantz, 1988). There is no evidence to indicate that P. marginemaculatus occurs beyond this line (Quintero, 1981 (Phillips et al., 2006), well beyond the ideal minimum sample size to obtain reliable results (Proosdij et al., 2016). Each model was assessed with the area under the receiver operating characteristic curve (AUC; Hanley & McNeil, 1982). AUC values represent a measure of the MaxEnt model's ability to discriminate between suitable and unsuitable areas in the modeled distribution (Anderson & Gonzalez, 2011). AUC values range from zero to one, with one indicating a perfect differentiation of suitable and unsuitable habitat. We compared predicted models against distribution literature for P. marginemaculatus (Quintero, 1981

| Population genetics
We sequenced a 1,238 bp region of the mitochondrial COI gene of 103 individuals. Sequences had an overall base composition of 24.6% adenosine, 25.2% cytosine, 16.5% guanine, and 33.7% thymine. Key Largo exhibited the lowest genetic diversity among regions, with the Lower Keys and mainland regions exhibiting considerably higher genetic diversity (Table 1). Additionally, a range-wide mantel test for IBD of geographic coordinates and Nei's distance was nonsignificant (Mantel r = −0.07, P = 0.637). A hierarchical AMOVA revealed population structure among populations, but not regions ( Table 2). The AMOVA indicated significant genetic structure between localities, suggesting that oceans limit dispersal. Stratified cross-validation of DAPC resulted in a mean successful assignment of 0.92% and 88.2% conserved variance with 20 principal components. All regions separated into distinct clusters, further evidencing regional genetic structure ( Figure 3). Furthermore, pairwise G ST showed relatively high divergence of Key Largo localities from the rest of the range, and low divergence between mainland and island sites, indicating ongoing gene flow (Figure 2).

| Ecological modeling
Species distribution modeling using MaxEnt found good model fit for climate-only models (mean AUC = 0.978 ± 0.02, n = 100 models; Figure 4a-c). Habitat suitability was highest in the Florida Keys ( Figure 4), but also extended to the southeastern end of mainland Florida. Areas of predicted suitable habitat on the southeastern mainland corresponded to the geologic features of the peninsula, which had a permutation importance of 5.8%. However, suitable habitat was identified primarily by environmental variables that contributed most to the model: precipitation of the coldest quarter (74.8% permutation importance) and mean diurnal temperature range (14.7%). Altitude also had predictive power with 2.4% permutation importance in the first set of models. Jackknife tests of variables in isolation from all others revealed that annual mean temperature had the highest training gain for models, followed by mean temperature of the driest quarter, mean diurnal temperature range, and mean temperature of warmest quarter.
Models including land use categories performed slightly worse (mean AUC = 0.873 ± 0.06, n = 100 models) than those without land use but appear to have refined the habitat suitability of P. marginemaculatus (Figure 4d-f). Land use categories had the highest permutation importance (50.5%) followed by mean temperature of the driest quarter (18%), mean diurnal temperature range (12.7%),  We saw an alarming 22%-34% human-induced decline in suitable habitat under the best fit model ( Figure 5; Tables 3 and 4).

TA B L E 2 Hierarchical analysis of molecular variance (AMOVA) using cytochrome c oxidase I (COI) sequences of
Models that included land use showed steeper declines of 29%-48% of suitable habitat due to human development (Tables 3 and 4). This indicates that, while climatically identified habitat shows considerable decline, human impacts particularly target land use habitat types important for the species.

| D ISCUSS I ON
Pine rocklands succeed to tropical hardwood hammocks after two or three decades of fire suppression (Loope & Dunevitz, 1981;Robertson, 1981), but it remains unclear if P. marginemaculatus occurs in both habitats. Both periodic seawater flooding and fire occur naturally in rockland habitats (Snyder et al., 1990) and P. marginemaculatus, like much of the fauna and flora species in these habitats (Hofstetter, 1975;Robertson, 1953Robertson, , 1962, have evolved to survive these stochastic events.
The mainland and Lower Keys show the greatest genetic diversity. This is likely because these two regions include the largest expanses of area, and both are represented by small patches of suitable habitat within an inhospitable matrix. Key Largo showed relatively low diversity, despite the island's size. However, only a small portion of Key Largo's land area is suitable habitat-most of the island is highly developed.
Genetic variation was significant between localities but not larger regions (Table 2), suggesting a more complicated genetic structure than our a priori regional delineations. Population-level genetic structure is aligned with intuition, considering all populations in our study occur on Keys or islands of habitat surrounded by human disturbance (Figure 2). Regional structure, however, is somewhat surprising. We selected regions a priori based on geologic features: the Key Largo region for its size and proximity to the mainland; the Upper and Lower Keys for their geologic variation in formation, if not timing thereof; and the mainland, as an obvious delineator from island populations. This generally aligns with the population structure of native Cerion land snails, which showed isolation between the Upper and Lower Keys (Shrestha et al., 2015). Species on the Lower Keys likely have a unique evolutionary history separate from the rest of the Keys.
Our DAPC analysis shows clear genetic structure across regions, with key populations all being closest related to mainland localities.
This pattern suggests that a mainland-metapopulation model described the landscape genetics of the Keys. As mentioned, genetic analyses of land snails in the Florida Keys also show an Upper-Lower Key division (Shrestha et al., 2015). But research on marine species like bicolor damselfish (Eupomacentrus partitus) and common reef sponge (Callyspongia vaginalis) found much lower levels of divergence between regions than we found for P. marginemaculatus (DeBiasse, Richards, & Shivji, 2010;Lacson et al., 1989) . This is likely because P. marginemaculatus is much more dispersal limited than a marine fish. Interestingly, the least divergent locality pairs of theses marine species included Key Largo, which is the same pattern we find in our study. This matching pattern supports the idea that currents play a major role in P. marginemaculatus genetic structure.
Pairwise G ST showed that the greatest divergence was between Key Largo and the other regions ( Figure 2). This could be caused by migrants from the Bahamas, the closest islands of which are ca. 100 km from Key Largo. Surprisingly, the mainland had low G st estimates with the Upper and Lower Keys. We posit that this is due to ongoing gene flow between these regions. Ocean currents likely push rafting P. marginemaculatus to the Lower Keys as might major weather events including hurricanes (Fleming & Murray, 2009).
Experiments with GPS-equipped buoys show that this is a major current pathway (Lee & Smith, 2002) and P. marginemaculatus, with the ability to breathe underwater, are aptly suited to survive the voyage (Hebets & Chapman, 2000). Other research has emphasized the importance of ocean currents in the Florida Keys in genetic structure, but this has been limited to marine species. For example, a study of three marine invertebrates showed high gene flow and connectivity across the Keys, with a pattern of southern migration (Richards et al., 2007). The salt marsh snake (Nerodia clarkii), which is somewhat restricted to shallow waters, showed genetic structure between Upper N and Lower Keys that was explained by IBD (Jansen et al., 2007). In general, our SDM predicted habitat suitability for P. marginemaculatus in parts of Monroe, Miami-Dade, Collier, Lee, and Hendry counties. This limited range generally matches museum records, with occasional records from counties as far North as Charlotte county (Quintero, 1981 (Hu et al., 2014). The mosquito A. aegypti showed practically no genetic structure among the Florida Keys, likely because they disperse via flight (Brown et al., 2013). Shrestha et al. (2015) proposed that the C. incanum spread southwesterly to colonize new Keys as they formed, with Lower Key populations being the youngest. Lastly, ant gut microbiota showed genetic structure between the Upper and Lower Keys (Hu et al., 2013).
Our models identified areas of tropical hardwood hammocks and pine rocklands as the most suitable habitat types for P. marginemaculatus. Indeed, this is where almost all of our observations occurred and is corroborated by published collection sites of the species for laboratory research (Hebets & Chapman, 2000;Weygoldt, 2000).
The two habitat types, being at relatively high elevation, are the primary targets for human development in Florida (Noss, LaRoe, & Scott, 1995;Snyder et al., 1990). Thus, our modeling results show that P. marginemaculatus suitable habitats are also areas where human disturbance has been, and continues to be, an imminent threat to the habitat and species.
Pine rockland forests, once common throughout southeastern Florida, are now one of the most threatened habitats globally, with at least 98% global loss (Noss et al., 1995). This includes ca. 8,000 ha in Everglades National Park and a mere 920 ha outside the park's boundary (Bradley, 2005). Pine rockland habitat was identified as one of the most suitable habitats based not only on our models that included land use categories but also those with only climate variables and geology. Pine rockland habitats occur on exposed limestone substrates where limestone rock outcroppings are common and provide important microhabitat for P. marginemaculatus (Chapin & Hebets, 2016;Hebets & Chapman, 2000). Human development, fire suppression, and climate change have altered or entirely removed many areas once dominated by this community (Kautz & Cox, 2001;Possley, Woodmansee, & Maschinski, 2008;Ross et al., 1994).  (Maschinski et al., 2011;Ross et al., 2009). In this sense, increases in the frequency and intensity of hurricane storm surges reshape pine rockland and similar vegetation communities (Ross et al., 1994 Caicos (Quintero, 1981). In particular, gene flow from the Bahamas and Cuba to Key Largo and the Lower Keys could be ongoing, but this remains unexamined. Temporally, we excluded historical samples with the concern that they would not inform modern biogeographic patterns. While historical material would increase sample sizes, our modern collections span the spatial range of Florida historical samples bar northern limits, where populations may no longer occur.
Sampling was uneven across sites, which could bias results. Lastly, we used only one sequence in our analysis. We chose a mtDNA marker because nDNA appears highly conserved in Amblypygi (Esposito et al., 2015), and mtDNA is thus, more informative. Caution should be taken, however, in interpreting results from only one mtDNA marker, and future research using more thorough genomic sequencing could reveal more high-resolution biogeographic patterns.

| Conservation
Our results point to a primary threat to the population health of P. marginemaculatus in the wild: habitat loss by human development.
Approximately 77,000 people permanently reside in the Florida Keys (US Census Bureau). Given the small land area of the archipelago, this accounts for an average density of 205.54 people per square kilometer and leaves approximately 150 km 2 uninhabited by humans, much of which is seasonally or permanently flooded habitat unsuitable for many terrestrial species. Furthermore, this does not include the impact from commercial, transportation, and utilities developments. Fortunately, much of the remaining habitat for P. marginemaculatus is within protected areas managed by federal and state agencies. This, however, does not protect habitats from all threats, including the impacts of increased hurricanes, sea-level rise, poaching, and microhabitat alterations.
Secondarily, P. marginemaculatus is collected from the wild for sale in the pet trade. While we do not have data on the number of individuals collected for sale in the pet industry, personal observations lead us to believe that it must be in the hundreds to thousands.
Wild populations of P. marginemaculatus would benefit from captive breeding that allows these fascinating animals to be kept as pets without reducing numbers in the wild. More research on wild populations is needed to assess population health, and we encourage researchers to conduct both field and laboratory studies on these fascinating organisms.