Genetic structure as a response to anthropogenic and extreme weather disturbances of a coastal dune dwelling spider, Arctosa sanctaerosae

Abstract The continued increase in the number of tourists visiting the Northern Gulf Coast (NGC), USA, in the last century, and the resulting sprawl of large cities along the coast, has degraded and fragmented the available habitat of Arctosa sanctaerosae, a wolf spider endemic to the secondary dunes of the white sandy beaches of the NGC. In addition to anthropogenic disturbance to this coastal region, hurricanes are an additional and natural perturbation to the ecosystem. The data presented here explore the status of populations of this species spanning the entire known range and the factors influencing population demography including anthropogenic disturbance and severe tropical storms. Using microsatellite markers, we were able to document the genetic structure of A. sanctaerosae, including current and historic patterns of migration. These results combined with ecological and census data reveal the characteristics that have influenced population persistence: ecological variables affecting the recovery of the population clusters after severe tropical storms, genetic fragmentation due to anthropogenic disturbance, and their interaction. These findings demonstrate the significance that the high traffic beach communities of the NGC and their impact on the once intact contiguous dune ecosystem have on recovery after severe tropical storms. Contemporary modeling methods that compare current and historic levels of gene flow suggest A. sanctaerosae has experienced a single, contiguous population subdivision, and the isolates reduced in size since the onset of commercial development of the NGC. These results point to the need for monitoring of the species and increased protection for this endangered habitat.

The interaction of the hurricanes and human-induced pressures on coastal taxa has not been extensively studied. Human modifications are known to severely limit the dynamic ability of an ecosystem for vary naturally (Nordstrom, 2000) or can amplify the impacts of naturally occurring stochastic disturbances (Jonzen et al., 2004;Schrott et al., 2005). These modifications include alteration to the supply and transport of the sand as well as climate change-induced sea level and surface temperature rise which is predicted to increase the severity of tropical storms (Komar, 1998;Slott et al., 2006).

The ecosystem along the Northern Gulf of Mexico Coast (NGC)
is experiencing growth in both of these classes of disturbance. The NGC including Northern Florida, Alabama, Mississippi, and Louisiana has had 35 hurricanes of category 3 or higher make landfall in the last 100 years and in that same period, the human population of the major cities along NGC has increased approximately 15-fold from 10,000 individuals to 150,000 (U. S. Census Bureau, 1990). This increased habitat fragmentation due to human encroachment is hypothesized to have reduced population size of flora and fauna through the creation of barriers to migration and subsequent subdivision of larger populations into a series of smaller populations. Small populations experience reduced population viability and persistence (Palstra & Ruzzante, 2008;Reed, 2010;Reed, Lowe, et al., 2003;Reed, O'Grady, Brook, et al., 2003;Saccheri et al., 1998). Small populations are also more susceptible to loss of genetic diversity due to random genetic drift, they maintain lower levels of fitness (Reed, 2005;, and have reduced adaptive potential (Blows & Hoffman, 2005;Reed, 2005;Reed, Lowe, et al., 2003;Reed, O'Grady, Brook, et al., 2003) when compared to larger populations.
The spatial heterogeneity in hurricane impacts suggests that spatial autocorrelations in population fluctuations (Burgman et al., 1993;McCarthy & Lindenmayer, 2000;Reed, 2004) might be especially important to metapopulation persistence in this system.
Because hurricanes reduce population size via the direct destruction of habitat, results from one habitat-specific species with similarities in vulnerability to storm-driven mortality should be relevant to the persistence of all species limited to that habitat.
Using the nocturnal burrow-dwelling wolf spider, Arctosa sanctaerosae, Gertsch and Wallace, 1935 (Araneae: Lycosidae) (Figure 1), we are able to investigate the effects of habitat fragmentation (human encroachment) and a catastrophe regime (severe tropical storms) as well as the interaction of these natural and anthropogenic disturbances. This species is an ideal subject to explore these questions, as it is entirely restricted to the secondary dunes in the coastal dune system of the NGC (McNatt et al., 2000). Given the general lack of invertebrate conservation work (Skerl & Gillespie, 1999) and the discrete generation length that spiders have, this taxon will provide insight into threats faced by other invertebrates and small vertebrate species of interest in the region (e.g., several species of the beach mouse Peromyscus polionotus, Osgood, 1907).
The purpose of this study was threefold: (a) explore the effects of severe tropical storms on spider density, and relate these to the physical attributes of the dune system, extent of disturbance, and the distance from the point of landfall of the tropical storm; (b) describe the genetic diversity and structure of these population clusters; and finally (c) to investigate gene flow, past and present, among these. The comparison of historic versus recent gene flow can potentially be used to differentiate isolation of populations caused by hurricanes and other forms of prehuman isolation (historically constant) versus contemporary causes of reduced connectivity including those caused by conversion of habitat by humans within the last 100 years. Assuming that hurricanes have occurred since the evolution of this species and that human disturbance on a large scale in the region started only a hundred years ago, the dominant force that is shaping the current status of the species and its populations should emerge. Gaining insight into the effects of environmental F I G U R E 1 Arctosa sanctaerosae, Gertsch and Wallace, 1935 (Araneae: Lycosidae) perturbations, anthropogenic encroachment, and their interaction would be invaluable informing the long-term conservation goals species who share the same endangered coastal dune habitat along the NGC.

| Population densities and growth rates
Densities were measured by hand collecting individuals inside three independent 12 m by 12 m quadrats (144 m 2 ) randomly placed within the secondary dunes of ten sites across the NGC during the summer months between 2003 and 2007 ( Figure 2). Counts were made of spiders at or near their burrows one hour after nightfall on three consecutive clear nights. The quadrats were located each year using GPS data and resampled. Growth rate (R) was calculated as R = ln (N 1 /N 0 ) for each site. All density estimates of zero were adjusted to 0.5 to ease statistical analysis based on the assumption that these population numbers were most likely not zero but too small to be detected.
These rates were compared before and after landfall of major Physical measures, dune height (dh) and the density of vegetation, were also quantified within the quadrats. These two measures were included due to their provision and maintenance of habitat for A. sanctaerosae as well as prey. Vegetation cover of dunes was quantified by counting the number of stems of sea oats, Uniola paniculata L. (Liliopsida:Poaceae) and similar vegetation within three independent, randomly chosen one square meter quadrats within the larger quadrats on each of the three sampling nights per site. Dune height was measured from the height of the apparent high tide mark on the active beach to the highest point on the secondary due. Model selection for explaining population growth rates was accomplished using an information-theoretic approach. Stepwise multiple regression was subsequently used to explore the relationships among these ecological factors and changes in them, as a function of hurricane landfall, associated with each site in order to identify those that best explained the variation in population decline and recovery observed after each of these two hurricanes.

| Census population size
Census population sizes were calculated for each of the five clusters (recovered by genetic clustering algorithms detailed below) by multiplying population densities across the clusters by the total habitat area. The total extent of habitat was estimated by calculating the area of secondary dunes within each cluster. This was done using GOOGLE EARTH by creating polygons and using the area tool to calculate the total area. The density per quadrat was then multiplied against total area to estimate census population size.

| Genetic sampling
Tissue samples from 20 individuals collected from each of the sites were collected between 1st June and 11th June 2007 to be used in genetic analyses. All individuals were stored at −80°C in 100% ethanol, and approximately 1 mg of tissue acquired from the legs was used for subsequent DNA extraction (Qiagen DNeasy kit). Target microsatellite sequences were amplified using 11 microsatellite primers, ten developed for the study species  and an addition developed for a sister taxa. Fluorescence-labeled fragments were visualized on an ABI 3130, and allele sizes were determined through comparison with a known size standard (GeneScan -500 ROX) using GENEMAPPER version 3.7. All scores were checked manually, and ambiguous fragments were reanalyzed.

F I G U R E 2
The five population clusters recovered from genetic data using GENELAND and STRUCTURE in circles with major human developments labeled: (a) Gulf Shores, AL, (b) Pensacola Beach, FL, (c and d) Destin, FL

| Population structure
Population differentiation was accomplished using the Bayesian clustering algorithms in the programs STRUCTURE (Pritchard et al., 2000) and GENELAND ver. 3.2.4 (Guillot et al., 2005;Guillot & Santos, 2009). Both of these programs have idiosyncrasies and can be used in concert to not only suggest a current number of clusters and assign individuals to them, but also to find support for a suggested model. GENELAND was run and included the correlated model, which assigns individuals to geographic clusters without prior knowledge of the site where that individual was sampled. This accounts for several factors relevant to this system. Spatially, we expect areas of intense human impact to create a barrier between populations and so we add into the model the set of georeferenced coordinates, in this case the location of our collection sites (Guillot & Santos, 2009). GENELAND and STRUCTURE both recovered five clusters and these assignments were used for all subsequent demographic analyses.

| Genetic diversity
FSTAT ver. 2.9.3.2 (Goudet, 2002) was used: (a) to test for Hardy- to calculate allelic richness (A r ) and (d) to calculate gene diversity (H s ). GENEPOP (Raymond & Rousset, 1995;Rousset, 2008) was used to test population differentiation and to estimate the number of null alleles. ARELEQUIN ver. 3.5.1.2 (Excoffier, 2010;Excoffier et al., 2005) was used to calculate expected (H e ) and observed (H o ) heterozygosities and pairwise F st values between populations. Other statistics have been suggested for estimating gene flow in microsatellites (Goldstein & Pollock, 1997;Goldstein et al., 1995;Weir & Cockerham, 1984;Zhivotovsky, 1999), but F st is the most commonly used metric and was employed here. Expected heterozygosity was computed among and within populations using Levene's method (1949) in the software package POPGEN (Yeh et al., 1999). Isolation by distance was tested using a Mantel test of log transformed linearized F st and geographic distance values was conducted among the clusters.

| Historic and recent geneflow estimation
Historic mutation scaled migration rates (M) were calculated using MIGRATE ver. 3.2.1 (Beerli, 2010). The estimates given are longterm averages heavily influenced by the recent past and are calculated over the last 4N e generations. The rate of recent gene flow was calculated using BAYESASS ver. 1.3 (Wilson & Rannala, 2003), which estimates migration rates within the last two to three generations using a Bayesian inference framework and gametic disequilibrium among immigrants and their descendants. Individuals from each of the five clusters were grouped suggested by GENELAND for both tests.

| Historic and recent bottleneck detection
To test for evidence of bottlenecks, M ratios were calculated (Garza & Williamson, 2001) in ARLEQUIN version 3.1 (Excoffier, 2010;Excoffier et al., 2005). The effectiveness of this method has been shown to be maximized when the bottleneck is older and lasted several generations before recovery (Williamson-Natesan, 2005).
BOTTLENECK (Cornuet and Luikart 1996;Luikart & Cornuet, 1998;Piry et al., 1999) was used for estimations of bottlenecks within the previous 4Ne generations, for this species approximately the last 100-500 years. The method employed in this program has been shown to have the ability to detect less severe and more recent reductions in population size.

| Historic and recent effective population size estimation
Historic effective populations sizes were measured using three methods. First, using the equation Ө = 4Neμ where Ө is the mutation scaled effective population size and μ is the mutation rate (Gaggiotti & Excoffier, 2000). Mutation rate was held constant and at a rate of 5x10 -4, and a coalescent approach was used to estimate of Ө over the last 4Ne generations in the program MIGRATE (Beerli, 2010;Beerli & Felsenstein, 1999).
Historic effective population size was also calculated using the methods of Hartl and Clark (1989) and Ohta and Kimura (1973).
These methods assume an infinite allele model (IAM) and a stepwise mutational model (SMM), respectively. Both methods hold N e as a function of H e . The mutational models estimate the upper and lower extremes of mutation, and so the true Ne is likely to be found between the two estimates (Busch et al., 2007). A paired t test was then carried out between the historic and recent estimates of Ne.
Recent effective population sizes were estimated using ONeSamp (Tallmon et al., 2008) as well as LDNe (Waples & Do, 2008). LDNe uses linkage disequilibrium to estimate N e . The major issue affecting its usefulness for this study is that linkage disequilibrium can be caused by inbreeding, substructure, or immigration. The first is known to be an issue in this species, and so the results of these tests bear careful scrutiny. ONeSAMP on the other hand uses eight different genetic parameters including: "the number of alleles divided by allele length range (Garza & Williamson, 2001), the difference of the natural logarithms of variance in allele length and heterozygosity (King et al., 2000), expected heterozygosity (Nei, 1987), number of alleles per locus, Wright's F IS (Nei, 1987), the mean and variance of multilocus homozygosity, and the square of the correlation of alleles at different loci (Hill 1981)" (Tallmon et al., 2008). This is expected to provide more accurate results, although high inbreeding levels (F IS ) may confound results. The results for both methodologies including estimates using both 0.01 and 0.02 for the lowest allele frequencies for the analysis in LDNe are reported.

| Population densities and growth rates
The mean population growth rate (R) across sites during the period between 2003 and 2004, which included the effects of Hurricane Ivan, was −2.703 (SE = 0.517). In the following one-year period (2004)(2005), the growth rate across sites was 1.738 (SE of 0.416).
This demonstrates a pattern of reduction and recovery of densities, assuming that the sites had been stable for some time before this.
This assumption is based on the extended period of time prior to Hurricane Ivan (8 years

| Census population size
Census population sizes were large and range from 71,000 to 315,000 individuals spread across the available habitat per cluster (Table 5). There was a significant difference between recent N e and estimates of N c for each of the clusters.

| Genetic diversity
The number of alleles per locus varied from 3 to 11, and null alleles were estimated to be <0.16% across all loci. Each of the eleven loci tested were found to be in gametic equilibrium. Tests of linkage disequilibrium, within and across populations, were all insignificant, satisfying the assumption that all loci are unlinked. The measures of diversity (allelic richness and gene diversity) decreased longitudinally, with the lowest levels of observed diversity in the westernmost populations. Heterozygosity, gene diversity, and allelic richness all followed this east -west pattern (Table 6). F is scores for each of the clusters loosely followed this pattern and ranged from 0.18 to 0.06. Pairwise F st scores ranged from 0.05 to 0.30, and all were significant ( Table 7). The Mantel test showed Pearson's r of .68 that was highly significant (p = .001) meaning that individuals are more likely to find mates from populations geographically close to themselves rather than at random across all populations.

| Historic and recent geneflow estimation
Using MIGRATE, a mean historic migration rate was calculated as the number of effective migrants (Table 8). BAYESASS estimates of recent migration were significantly lower, suggesting these clusters showed higher levels of isolation in the recent past. A paired t test of the means found the difference between historic and recent migration rates of A. sanctaerosae to be statistically significant (df = 4, p < .002).

| Historic and recent bottleneck detection
Tests of heterozygosity using the statistical package BOTTLENECK showed no significant excess in the five clusters. However, the M TA B L E 5 Estimates of historic and recent effective population size (N e ) and census population size (N c ) for each of the five population clusters recovered from genetic data using GENELAND and STRUCTURE

| Historic and recent effective population size estimation
Historic estimates of N e (Table 5) across methods varied. The mean of the three methods was used for subsequent analyses and was significantly higher than the estimate of recent effective population size and the lower C.I. did not cross the upper C.I. of the estimates of recent N e .
Estimates of recent effective population size from LDNe varied widely and included negative values as well as upper estimates that reached infinity. Varying the lowest allele frequency from 0.01 to 0.02 did not narrow the 95% confidence intervals (Table 5).
ONeSAMP gave results that appeared to be more biologically significant. Each of the clusters has a N e of less than one hundred individuals with the exception of cluster 5. Mean estimates range from 32 to 153 individuals with less variation around the estimates using this method compared to LDNe. Only the ONeSAMP estimates were used for comparisons against historic N e .

| D ISCUSS I ON
The impact of hurricanes on the density of A. sanctaerosae was direct, density-independent, and the best predictor of the magnitude of population reduction was the distance from the point of landfall.
The effects of Hurricane Katrina showed a gradient of severity ranging from the most severe in the western sites to no measurable effect in sites beyond 283 km from the site of landfall. Hurricane Ivan, by contrast, was centered directly over the western sites and was of sufficient severity to affect every site.
More specifically, population clusters 2 and 4 are receiving more effective migrants than the other three clusters. The rate of emigration into Cluster 4 is best explained by its geographic isolation and small size in relation to the other populations, suggesting this locality may serve as a sink population that requires immigration to persist. The severity of effects experienced by Cluster 2 due to hurricanes in the recent past could explain the high number of migrants from Cluster 1. Additionally, the spatial heterogeneity of hurricane effects is not limited to the distance from landfall.
The hurricane's leading edge has higher wind speeds relative to the trailing edge and is formed on the eastern side of the storm in the northern hemisphere. This frequently leads to an increased number of tornadoes spawned by the storm and more destructive power to the east of the storms point of landfall (Williams & Sheets, 2001). Hurricane Ivan made landfall due west of Cluster 2, and thus, the severity of effects may have been higher than that at Cluster 1, which with its high migration rate, provided a rescue effect. This migration appears to be elevated after major storm events.
The model that best explained population recovery in the year following a storm included both distance from landfall of the hurricane and the recovery of dune height. This dune habitat that is required for population recovery is often removed or highly partitioned by commercial development along the NGC. Five genetically distinct population clusters were recovered, and the barriers separating them were the high traffic beach communities of the NGC including the following: Gulf Shores, AL; Orange Beach, AL; Pensacola Beach, FL; Destin, FL; and Panama City Beach, FL; (Figure 1). This suggests that the commercial development of the NGC and the storm-related population reductions work in concert to reduce the size of population clusters, reduce the cluster's ability to recover from storms, and reduce the population cluster's connectivity.
This type of habitat fragmentation shapes population genetic structure by reducing population cluster connectivity. In addition to low rates of gene flow, reductions in population size and loss of available high-quality habitat can lead to genetic isolation and eventual extinction of populations (Reed 2008). Previous studies of spider population demography in fragmented habitats have found levels of structure and inbreeding similar to those found in the current study suggesting a decreased role of ballooning . The genetic isolation seen in the population clusters of Arctosa sanctaerosaeis is likely due to a decline in effective population size and effective migration in the recent past with F is scores ranging from 0.18 to 0.06, which suggests low amounts of dispersal across available habitat.
Historically, these population clusters experience severe reduction in size due to hurricanes or other catastrophic events; however, they had the ability to recover within one to two generations through migration to recently vacated habitat patches. This migration erased enough of the genetic signature of the reduction that it is undetectable using the current methodology. The genetic data suggest that there has been a significant decrease in the amount of migration between the five clusters within the last 100 years. However, this lack of migration in the recent past combined with the correlation of the patterns of genetic structure and the timing of the development of high traffic beach communities supports placing the timing of the subdivision to within the last 100 years. Once subdivided, the migration between clusters was reduced, and genetic drift, direct removal of diversity due to severe tropical storms, and commercial development led to declining population sizes and reduced diversity.
The observed genetic variation decreases along a geographic gradient from east to west with the alleles found in western populations as subsets of the alleles found in the eastern populations.
Allelic richness and gene diversity are highest in the east and decrease as you move westward. The best three possible explanations for the higher genetic diversity seen in the eastern clusters are (a) the relatively few numbers of hurricanes to make landfall in the east in recent years, (b) significantly fewer high traffic beach communities, or (c) the possibility that A. sanctaerosae has expanded its distribution westward over time from the site of its divergence from its sister taxa. Hurricane Michael made landfall in October of 2018 in the eastern most cluster, which exhibits the highest richness and diversity, at Mexico Beach, Florida, and directly on several collection sites. The effects this will have on the species' persistence would not be known for decades.
All of this data combined paints a picture of a metapopulation that has evolved in the presence of severe tropical storms and has been recently impacted through commercial development that has created barriers to repopulation, gene flow, and directly lowered N e . With the severity of tropical storms predicted to increase over time due to elevations in surface temperature of the Gulf of Mexico caused by climate change (Goldenberg et al., 2001;Slott et al., 2006;Webster et al., 2005), population connectivity will become increasingly important for population and even species persistence.
Moving forward, conservation measures must recognize the need for a contiguous dune system, not only for the physical preservation of the coastal dune ecosystem, but also for maintaining population structure and connectivity of this and presumably other species. Importantly, increased continuity and protection of these dunes will also lead to a healthier dune system that prevents erosion and inland destruction in the face of tropical storms. It must be recognized that a contiguous dune system, while beneficial to the focal species, also serves the broader goals of ecological and commercial interests along the NGC. By removing the dunes for commercial development, we have increased erosion and must reclaim the beach sands regularly as well as replace structures damaged during major tropical storms at great cost. It makes ecological as well as financial sense to restore these dune systems to their original state and allow natural processes to maintain them. Only the halting and/or reversal of the current developmental trends, including the complete removal of barriers interrupting the corridors for inter-population migration, will result in the long-term persistence of A. sanctaerosae.

ACK N OWLED G M ENTS
I would like to thank all the colleagues who graciously reviewed the manuscript as well as Brice Noonan for laboratory assistance.

CO N FLI C T O F I NTE R E S T
There are no sources or any potential sources of conflict of interest.
There is no interest or relationship, financial or otherwise, that might be perceived as influencing the objectivity of this work.