Using genetics to inform restoration and predict resilience in declining populations of a keystone marine sponge

Genetic tools can have a key role in informing conservation management of declining populations. Genetic diversity is an important determinant of population fitness and resilience, and can require careful management to ensure sufficient variation is present. In addition, population genetics data reveal patterns of connectivity and gene flow between locations, enabling mangers to predict recovery and resilience, identify areas of local adaptation, and generate restoration plans. Here, we demonstrate a conservation genetics approach to inform restoration and management of the loggerhead sponge (Spheciospongia vesparium) in the Florida Keys, USA. This species is a dominant, habitat-forming component of marine ecosystems in the Caribbean region, but in Florida has suffered numerous mass mortality events. We developed microsatellite markers and used them to genotype sponges from 14 locations in Florida and a site each in The Bahamas, Belize and Barbuda. We found that genetic diversity levels were similar across all sites, but inbreeding and bottleneck signatures were present in Florida. Populations are highly structured at the regional scale, whilst within Florida connectivity is present in a weak isolation by distance pattern, coupled with chaotic genetic patchiness. Evidence of a weak barrier to gene flow was found in Florida among sites situated on opposite sides of the islands in the Middle Keys. Loggerhead sponge populations in Florida are vulnerable in the face of mass mortalities due to low connectivity with other areas in the region, as well as distance-limited and unpredictable local connectivity patterns. However, our discovery of Florida’s high genetic diversity increases hope for resilience to future perturbations. These results provide valuable insight for sponge restoration practice in Florida.


Introduction
Population declines in keystone species have a number of negative impacts on associated communities, ecosystem functioning, and the provision of ecosystem services (Sweeney et al. 2004;Hicke et al. 2012; Thomson et al. 2015;Sorte et al. 2017). Genetic factors are a significant determinant of population health and fitness, and can influence both the longevity of populations and the success of conservation strategies (Frankham 2005). However, genetic information is unavailable for the vast majority of species, and thus the application of conservation genetics theory to practice has been limited (Shafer et al. 2015;Taylor et al. 2017).
Declining populations are vulnerable to low genetic diversity due to the effects of genetic drift, in which rare alleles have a higher probability of being lost due to random chance in smaller populations. These effects are amplified considerably in populations that experience a rapid decline, or bottleneck, through which substantial genetic variation is randomly eliminated in a short space of time (Sbordoni et al. 1986;Bellinger et al. 2003;Bristol et al. 2013). This threatens population survival, as genetic diversity is an important determinant of long-term population persistence (Frankham 2005). Indeed, high genetic diversity bolsters the resilience of populations, because they harbour a higher adaptive capacity with which to respond to perturbations such as disease, environmental change, or declining environmental conditions (Hughes and Stachowicz 2004;Ehlers et al. 2008;Evans et al. 2017). Low genetic diversity is also related to inbreeding depression, where recessive deleterious alleles are more likely to combine within individuals and reduce fitness (Whitlock 2000;Reed and Frankham 2003;Charlesworth et al. 2009), further compromising the long-term prospects for survival of the population.
Connectivity-the movement of individuals or propagules among populations-is an important counterforce against declining population size, low genetic diversity and local extinction. A well-connected population receives a regular supply of immigrants, thus boosting population size. Crucially, if these migrants successfully reproduce, they can help replenish the gene pool with new alleles, thus countering the effects of genetic drift through gene flow (Garant et al. 2007;Saenz-Agudelo et al. 2011;Frankham 2015). Conversely, isolated populations with little connectivity are more vulnerable to extinction due to limited immigration and gene pool restriction (van der Meer et al. 2013). Assessing levels of genetic connectivity among geographical sites is therefore another key step in managing vulnerable populations.
An interesting case study to apply such genetic information to conservation practice exists among sponge populations in Florida Bay and the Florida Keys (USA). In nearshore hard-bottom habitats in this area, sponges form a dominant component of benthic communities (Chiappone and Sullivan 1994;Tellier and Bertelsen 2008), and perform a number of vital functional roles and ecosystem services. Given their high relative biomass, they provide the majority of architectural complexity and habitat structure in the area (Herrnkind et al. 1997). This is especially important given that Florida Bay is a nursery area for a number of economically important fish and invertebrate species, including snapper (Lutjanus spp.), stone crabs (Menippe mercenaria), and Caribbean spiny lobsters (Panulirus argus). Several species of sponge are themselves 1 3 the target of commercial fisheries in the region . Moreover, sponge endosymbionts are important in creating soundscapes that form an acoustic cue for larval settlement in a variety of taxa (Butler et al. 2016). As filter feeders, sponges drive nutrient cycling dynamics in the area (Fiore et al. 2017;Hoer et al. 2018;Valentine and Butler 2019), and contribute to the maintenance of water quality (Peterson et al. 2006;Butler et al. 2018).
Sponge communities in the Florida Keys have suffered a number of mass mortality events (Butler et al. 1995;Stevely et al. 2010;Wall et al. 2012) associated with recurring blooms of the cyanobacteria Synechococcus spp. (Fourqurean and Robblee 1999;Berry et al. 2015), as well as stochastic cold weather events (Colella et al. 2012) and storm damage (Stevely et al. 2010). These mass mortalities have had dramatic consequences for the ecosystem, including declines in local juvenile lobster populations (Butler et al. 1995;Herrnkind et al. 1997), increased susceptibility to further cyanobacterial blooms (Peterson et al. 2006;Wall et al. 2012) and diminished underwater soundscapes predicted to impact larval recruitment from a variety of taxa (Butler et al. 2016). Furthermore, sponge population recovery is potentially forestalled by limited dispersal, as adults are sessile, and sponge larvae are generally short-lived, with larval durations of a few hours to a few days before settlement (Maldonado 2006;Maldonado and Riesgo 2008).
Due to their keystone role in the ecosystem and the impacts of their decline, sponge restoration work has been undertaken in the area for a number of years, where healthy sponges have been fragmented and translocated to areas that have suffered mortalities (Butler et al. 2016;Valentine and Butler 2019). However, cyanobacterial blooms and sponge mass mortalities continue to recur across different areas of the Florida Keys and Florida Bay. Coupled with work to identify the proximal causes of sponge mortality and the implementation of habitat improvement measures, understanding the genetic status of the populations is imperative for future restoration and management planning. In addition, investigating connectivity patterns will aid understanding of source-sink interactions across the Bay, and identify priority areas for restoration.
In this study, we investigated these topics in the loggerhead sponge, Spheciospongia vesparium (Lamarck, 1815). Spheciospongia vesparium is common throughout Florida Bay and has the largest biomass of all sponge species in the Bay (Tellier and Bertelsen 2008). It is also found on reefs and in lagoons throughout the Greater Caribbean region. Reproduction and larval biology have not yet been studied in S. vesparium, therefore limiting our ability to predict dispersal and population genetic pattners. However, studies of other Clionaidae species suggest that varied reproductive characteristics exist within the family: sexual and asexual reproduction have both been observed (Rosell and Uriz 2002;Schönberg 2002;Maldonado and Riesgo 2008), and similarly, gonochorism and hermaphroditism are also both found within the family (Piscitelli et al. 2011;González-Rivero et al. 2013). The Clionaidae are oviparous (i.e., broadcast spawning of both the sperm and eggs) (Ereskovsky 2018), and fertilization and larval development are mainly external, although in Cliona vermifera eggs are fertilized internally and the zygote released (Bautista- Guerrero et al. 2014). Larvae are lecithotrophic (i.e., do not feed), and larval duration is shortin Cliona viridis, it was estimated at < 10 days (Mariani et al. 2000). Clionaidae larvae have so far been observed to show weak swimming ability, with crawling behaviour common (Mariani et al. 2000(Mariani et al. , 2001.
Here, we aimed to describe patterns of genetic diversity and genetic connectivity in S. vesparium at hard bottom sites across the Florida Keys. In addition, we sampled three other locations in the Greater Caribbean to act as comparative populations, and to observe drivers of population structure at the regional scale.

Sample collection and preservation
We collected S. vesparium samples from shallow water sites (< 2 m depth) in four main localities: the Florida Keys/ Florida Bay (USA), Abaco Island (The Bahamas), Barbuda, and Caye Caulker (Belize) ( Fig. 1; Table 1). We sampled a number of sites across the Upper, Middle and Lower Florida Keys: 12 sites on the Florida Bay side of the Keys and 2 collection sites on the Atlantic side (Table 1; Fig. 1). Our sites in Florida included both those that have previously been affected by cyanobacterial blooms and mortalities, and those that have not. At each site in Florida we sampled between 10 and 32 individuals (average of 18.6 ± 1.2 SEM), and in Abaco, Barbuda, and Caye Caulker we sampled 12, 20, and 10 individuals, respectively (Table 1). We avoided sites where restoration work had taken place in order to observe the natural patterns of population structure and genetic diversity as far as possible. We collected small tissue fragments (~ 2 cm 3 ) and immediately transferred the samples into 95% ethanol, which was renewed after 24 h.

Microsatellite development
For this study, we characterised twelve new tri-and tetra-nucleotide microsatellite loci (see Supplementary Material for full details of the methods). In brief, DNA from a single S. vesparium sample collected from Long Key (Bay-side) was sequenced using Illumina MiSeq 2 × 250 base pair technology. We then processed the sequence reads using the Palfinder Galaxy bioinformatics pipeline (Griffiths et al. 2016) to quality filter the data, screen for microsatellites and design primers. We tested 36 loci, of which 12 could be successfully amplified and scored, and were subsequently used in this study.

DNA extraction and genotyping
We checked sponge tissue samples under a dissecting microscope to remove any visible endosymbiotic invertebrates, and then extracted total DNA using the DNeasy® Blood and Tissue Kit (Qiagen). We combined 10 of the 12 microsatellite primer pairs in two multiplex (5-plex) PCRs using the fluorophores 6-FAM and HEX (Table S1), and ran two primer pairs, Vesp36 and Vesp9 in singleplex PCRs due to problems encountered in multiplexing these loci. We utilized a three-primer universal tail approach for fluorescent labelling PCR products, as described in Blacket et al. (2012) and Culley et al. (2013). We carried out PCRs using the Type-it® Microsatellite PCR Kit (Qiagen) with the following cycling conditions: 95 °C for 5 min, 28 × (95 °C for 30 s, 60 °C for 90 s, 72 °C for 30 s), 60 °C for 30 min. For any amplification failures, PCRs were repeated in singleplex reactions with lowered (50-59 °C) annealing temperatures.
We sized PCR products by capillary electrophoresis using a 3730 DNA Analyzer (Thermo Fisher Scientific) with GeneScan™ 500, 600 or 1200 LIZ® size standard (Thermo Fisher Scientific), or a homemade ROX-based size standard. On all plates, we included both positive and negative controls. We scored alleles using Genemapper® v3.7 software (Thermo Fisher Scientific), and corrected allele sizes according to the positive controls to account for differences in allele length based on the machine or size standard used. We then binned alleles using the R package 'MsatAllele' v1.02 (Alberto 2009) in RStudio v3.0.3 (R Core Team 2014).

Quality control and summary statistics
We calculated the probability of linkage disequilibrium between pairs of loci using Genepop on the Web v4.2 (Raymond and Rousset 1995;Rousset 2008), with p values corrected for multiple tests using the false discovery rate procedure of Benjamini and Yekutieli (2001), as calculated using the R function p.adjust. We estimated null allele frequency at each locus using the EM algorithm (Dempster et al. 1977) in FreeNA (Chapuis and Estoup 2007

Genetic diversity, inbreeding and bottlenecks
We used Genodive v2.032b (Meirmans and Van Tiendener 2004) to calculate observed heterozygosity (H O ) and gene diversity/ expected heterozygosity (H S ; Nei 1987). We also tested for probability of departure from Hardy-Weinberg Equilibrium (HWE) in Genodive using the AMOVA (least squares) method and 50,000 permutations (p values corrected for multiple tests using the Benjamini and Yekutieli method, calculated as previously). We calculated average allelic richness and private allele richness rarefied to the lowest sample size (maximum g = 10) in ADZE v1.0 (Szpiech et al. 2008). We repeated these analyses with all the Florida sites grouped as one population and each separately. We estimated inbreeding coefficients (Avg Fi) in INEst v2.1 (Chybicki and Burczyk 2009), correcting for the presence of null alleles. The program includes three possible parameters that can affect inbreeding coefficient estimation: null alleles ('n'), inbreeding ('f') and genotyping failure ('b'). We ran the individual inbreeding model (IMM) for all combinations of these parameters and calculated the Deviance Information Criteria (DIC) for each run to determine the best model fit for the data. We ran the model using 500,000 Monte Carlo Markov Chain (MCMC) cycles with 50,000 burnin cycles.
We used INEst to find evidence of genetic signatures of recent population bottleneck events. The program implements two tests; the first identifies heterozygosity excesses in respect to allelic richness (Cornuet and Luikart 1996), and the second identifies M-ratio (mean ratio of allelic richness to allelic size range) deficiencies (Garza and Williamson 2001). Both phenomena have been observed when populations experience rapid reductions in size. We used the two-phase mutation model, and tested significance using a Wilcoxon signed-rank test with 1000 permutations.

Genetic connectivity patterns
We estimated genetic differentiation among sites by calculating pairwise F ST (Wright 1943(Wright , 1949 and D (Jost 2008) in Genodive v2.032b. For F ST values, we tested their significance in Genodive using 50,000 permutations, and corrected significance for multiple tests using Benjamini-Yekutieli correction as described above.
We used two different approaches to infer the number of population clusters ('K') in the data. Firstly, we used the Bayesian individual-based assignment model implemented in the 'Geneland' package v.4.0.8 (Guillot et al. 2005;Guillot et al. 2008) in RStudio, which uses spatial and genetic data to infer K and calculate the probability of individual assignment. Due to the assumptions of this model, we used only the seven loci that did not deviate from HWE in the majority of the sites, and deleted samples in which missing data was present in the majority of the HWE loci (n = 285). We first ran the no-admixture model to obtain estimates of cluster membership and allele frequencies. We used the uncorrelated allele frequencies, spatial and null allele models, and ran the program with 1,000,000 MCMC iterations, 100 thinning and 1000 burnin, and uncertainty on coordinates set to 0.0005. We set the maximum number of nuclei to 855, and the maximum rate of the Poisson process to 285. We repeated this with K ranging from 1 to 17, with 10 independent runs for each value of K. We then ran the admixture model using the estimates obtained from the noadmixture run with the highest average posterior probability. For the admixture model, we used 1,000,000 MCMC iterations, a thinning of 100, and a burnin of 1000. We extracted the q-matrix of estimated individual membership proportions to each of the detected clusters, and used Distruct v1.1 (Rosenberg 2004) to graphically display the results.
We used Flock v3.1 (Duchesne and Turgeon 2012) as an alternative method to infer membership to population clusters. This method estimates K and partitions samples into K clusters based on iterated reallocation, uses no a priori information on sampling location, and does not assume populations are in HWE. We tested K from 1 to 17 in 50 independent runs per value of K, and ran each model with 20 iterations (i.e. 20 rounds of reallocation). We used plateau analysis based on log likelihood difference (LLOD) scores, as described by Duchesne and Turgeon (2012), to infer the most likely value of K. We carried out hierarchical clustering approaches for both the Geneland and Flock analyses by first running the models using all sites, and then repeating the process on any multi-site clusters identified.
We used Discriminant Analysis of Principle Components (DAPC) (Jombart et al. 2010) as implemented in the package 'adegenet' v.2.1.1 (Jombart 2008) in RStudio to examine genetic variation among the sites based on allele frequencies. We used the function optim.a.score to calculate the optimum number of principle components (PCs) to retain in the analysis to prevent over-fitting of the model, whilst preserving the maximum discriminability. We included all sites in the first instance, and then conducted a further analysis on the Florida sites alone to examine the presence of fine-scale structure.

3
We looked for evidence of barriers to gene flow among the Florida and The Bahamas sites using the software Barrier v2.2 (Manni et al. 2004). We excluded the Barbuda and Belize sites from this analysis because of the large geographic distances separating them from the other sites, as this does not offer an appropriate theoretical framework to search for oceanographic barriers. Barrier uses the spatial coordinates of the sampling sites and Delauney triangulation to partition geographic space into polygons, creating a Voronoï tessellation map with each site contained within a single polygon whose edges border neighbouring adjacent sites. Monmonier's (1973) maximum difference algorithm then uses this tessellation map along with a genetic distance matrix (Jost's D) to detect genetic discontinuities among neighbouring sites. We assessed the robustness of the computed barriers by repeating analysis on 100 resampled bootstrap D matrices. We created the resampled bootstrap matrices in the R package 'diveRsity' v1.9.90 (Keenan et al. 2013). We computed increasing numbers of barriers until bootstrap support fell below 50%, reaching a maximum of three barriers. Following computation of barriers, we used AMOVAs to examine the partitioning of genetic variation across barriers.
To test the presence of genetic isolation by distance (IBD) within Florida, we performed a Mantel test to detect association between matrices of linearised pairwise genetic distances (F ST /[1 − F ST ]) and the logarithm of geographic distances. We calculated least-cost oceanographic distances between sites (i.e., the shortest distance possible, excluding landmasses) using 'marmap' v0.9.5 (Pante and Simon-Bouhet 2013) in RStudio, and carried out the Mantel tests in 'ade4' v1.7-10 ( Dray et al. 2007) in RStudio, with 9999 permutations to calculate significance.
We used Geneclass2 v2 (Piry et al. 2004) to detect first generation migrants among the sampling locations, and their putative population origins. We used the Bayesian criteria of Rannala and Mountain (1997) for likelihood estimation, and the Monte Carlo method of Paetkau et al. (2004) for probability computation, with the L home criterion, as source populations for all individuals were unlikely to have been sampled. We used a significance threshold of p < 0.01 and carried out simulations with 10,000 individuals.

Quality control and summary statistics
In total, we collected samples from 326 individuals across 17 sites (Table 1). Twentytwo samples were removed from the final dataset due to amplification failure in over 50% of the loci, leaving 304 individuals. Two individuals from the Lakes Passage had identical genotypes, one of which was removed from the dataset for analysis, yielding 303 individuals. Following correction for multiple tests, no significant linkage disequilibrium was found between pairs of loci. Null allele frequency was high in some markers (Tables S1, S2); however, post hoc analysis showed that the null allele-corrected global F ST value was only marginally higher (+ 0.002) than the uncorrected value when all loci were included in the analysis (Table S2). Furthermore, the difference between the uncorrected and corrected F ST did not increase as more loci were added (r 2 = − 0.03608, p = 0.4504), and therefore all loci were retained for the population genetics analyses. The number of alleles per locus over all sites ranged from 4 (Vesp23) to 27 (Vesp30) ( Table S1).

Genetic diversity, inbreeding and bottlenecks
Genetic diversity (allelic richness and gene diversity) was slightly lower at the Florida and Barbuda sites compared to The Bahamas and Belize. However, overlapping error bars among many of the sites indicate that this is only significant for a few Florida sites (Table 2; Fig. 2; Table S3). Broadly, genetic diversity can therefore be considered to be the same across Florida and non-Florida sites. Average rarefied allelic richness ranged from 3.408 (Pigeon Key) to 4.399 (Long Key-Bay-side) and gene diversity (H S ) ranged from 0.569 (in Pigeon Key) to 0.735 (in Belize) ( Fig. 2; Table 2; Table S3).
Observed heterozygosity over all loci varied from 0.251 (Craig Key-Atlantic) to 0.504 (The Bahamas) ( Table 2; Table S3). All sites had lower than expected levels of heterozygosity (Table 2; Table S3), and significant departures from HWE were found in a number of loci and populations (Table S4). The DIC analysis in INEst determined either the 'nfb' (null allele, inbreeding and genotyping failure) or 'nb' (null allele and genotyping failure) models to be the best fit for the sites in this study (Table 2). This indicates that null alleles and genotyping failure would affect inbreeding coefficient estimations in all the sites, but in ten of the sites, inbreeding was also an influential component of the model. The null allele-corrected inbreeding coefficients were positive in all locations, ranging from 0.036 (The Bahamas) to 0.343 (Craig Key-Atlantic). However, the posterior 95% probability intervals included zeros at all sites, and therefore F IS cannot be considered to be significantly above zero. When the Florida sites were grouped together as a single population, however, the posterior 95% probability interval was above zero, which may indicate significant inbreeding across the area. We found deficiencies in M-ratios at four sites, indicating the presence of recent bottleneck events (Boca Chica Channel, p = 0.0385; Little Crane Key, p < 0.001; Craig Key (Atlantic, p < 0.001; Long Key (Atlantic), p < 0.001). However, none of the sites showed significant heterozygote excess in comparison to allelic richness.

Genetic connectivity patterns
Pairwise F ST ranged from -0.019 (no differentiation) between Craig Key (Atlantic) and Long Key (Atlantic), to 0.273 (great differentiation) between Pigeon Key and Barbuda (Table 3). Among the four regional locations (Florida, The Bahamas, Barbuda, Belize), F ST values were large and significant, showing strong differentiation. Among sites within Florida, F ST values were lower (≤ 0.116), but significant differentiation was present between many pairs of sites. In general, higher differentiation could be observed among Fig. 2 Average allelic richness and private allelic richness per site (rarefied to maximum sample size g = 10). Error bars ± 1 SE Table 3 Pairwise Upper and Lower Keys sites than comparisons involving the Middle Keys sites, but patchiness can be observed throughout. Patterns of D were similar, and ranged from -0.035 (between the Craig Key and Long Key Atlantic sites, as previously) to 0.668 (between The Bahamas and Waltz Key) (Table 3). Private alleles were present at many sites (Table 2;  Table S3), and average private allelic richness was higher among the non-Florida sites (Table 2; Fig. 2).
Using Geneland, K = 4 was found for each independent run, with each regional location forming a separate population cluster (Fig. 3). In contrast, Flock showed strong evidence for K = 2. Samples were broadly partitioned into a Florida cluster and a cluster comprising individuals from Belize, Barbuda and The Bahamas. Two individuals from Florida (Craig Key Atlantic and Lakes Passage) fell into the Belize, Barbuda and The Bahamas cluster; otherwise, clustering was concurrent with sampling locations. When a second Flock analysis was carried out on the Barbuda, Belize and The Bahamas cluster, the samples were partitioned into K = 3 concurrent with sampling locations. When repeating the models for only the Florida samples, the Geneland model was unable to converge, indicating that K = 1 or the presence of strong isolation by distance in the data. Similarly, no plateau was obtained in Flock, indicating K = 1.
The DAPC showed clear separation of the Barbuda and The Bahamas sites from all other sites (Fig. 4a). All the Florida sites clustered together, with inertia ellipses showing substantial overlap among sites. The Belize site clustered closely to the Florida sites, with some Belize samples showing overlap with the Florida point cloud. In the Florida-only DAPC analysis, no clustering patterns were present, but points from sites more closely situated geographically tended to be closer together in the DAPC plot (Fig. 4b).
In the PCoA carried out on all sites (Fig. 5a), the first axis separated Florida from The Bahamas, Belize and Barbuda, and the second separated the Upper Keys from the Lower Keys and Atlantic side of the Middle Keys; the Bay side Middle Keys were distributed among both. In the Florida-only PCoA (Fig. 5b), points were distributed in a loose isolation by distance fashion, but notably the sites on the Atlantic side of the Middle Keys (Long Key and Craig Key) were clustered with Waltz Key, and separated from the sites on the Bay side of the Middle Keys. When the analysis was replicated with Jost's D instead of F ST , the patterns observed were very similar (data not shown).
The AMOVA showed that 18.2% of the total variation was found among the four main locations, while 2.4% was found among the sites within Florida. 30.9% of variation was found among individuals within sites, while 48.5% was within individuals (Table 4). Barrier software suggested the presence of two barriers with high bootstrap support: the first was a barrier between Florida and The Bahamas, with a bootstrap score of 100%. A second barrier separated the Atlantic sites from their adjacent Bay-side sites in the Middle Keys in Florida (Fig. 6). As barriers are computed based on the tessellation map, this barrier comprised a number of polygon edges, which showed bootstrap support ranging 39-99% (Fig. 6). A further barrier was estimated to separate the Middle and Lower Keys sites, however, this had low bootstrap support (10-54%) (Fig. 6). AMOVA analysis of sites separated by barriers showed that a large proportion of genetic variation was present between Florida and The Bahamas (20.1% of genetic variation, F ST = 0.201), while only a small amount of genetic variation was found between the Atlantic and adjacent Bay-side sites in the Middle Keys (3.3% of genetic variation, F ST = 0.033) ( Table 4). In both cases, more variation was found across the barrier than among sites on the same side of the barrier (Table 4).
Isolation by distance within Florida was significant, but the effect size was relatively small (r = 0.229, p = 0.031) (Fig. 7). We obtained comparable results when repeating the analysis with Jost's D (r = 0.225, p = 0.033).
Three putative first generation migrants were detected. All potential migrants were found within Florida sites, and all originated from other Florida sites. Two migrants were found at Long Key (Bay), with origins from Waltz Key (p = 0.0007; distance 89 km) and Fiesta Key (p = 0.0019, distance 5 km), and the third migrant was found at Little Crane Key with inferred origins of Kemp Channel (p < 0.0001, distance 12 km).

Genetic diversity and bottlenecks
Genetic diversity was similar throughout all of the sites sampled in Florida and the Caribbean, and was comparable to levels observed in other demosponge species (Chaves-Fonnegra et al. 2015;Riesgo et al. 2019). This implies that recurring mass mortality events have not significantly reduced genetic diversity in Florida, however, pre-mortality data is not available to confirm this hypothesis. Nevertheless, this study does provide a baseline with which future assessments of genetic diversity can be compared. We did not find signatures of genetic bottlenecks in sites that have been affected by cyanobacterial blooms. However, bottleneck signatures were present in four Florida sites that have not been affected by cyanobacterial blooms. These sites may have suffered unrecorded mortality events due to a different cause, such as disease, climatic variation, or hurricane disturbance. As sponges rapidly disappear once dead, and leave no visible skeleton, mass mortalities in sponges can be overlooked unless specific, regular monitoring is undertaken (Wulff 2006).
Levels of genetic diversity (allelic richness, gene diversity) in other sponge populations that have experienced mass mortalities vary by species. Spongia officinalis has high genetic diversity with no bottleneck signatures (Dailianis et al. 2011), whilst the opposite was found for a congener, S. lamella (Pérez-Portela et al. 2015). In Ircinia fasciculata, evidence of bottlenecks have been found at many (but not all) sites known to have suffered mortalities (Riesgo et al. 2016). In other species within the Florida reef tract, such as the coral Acropora cervicornis and sea urchin Diadema antillarum, genetic diversity was similar to other Caribbean sites tested (Chandler et al. 2017;Drury et al. 2017), despite mass mortality events. High genetic diversity, despite recent mass mortalities, may be due to high levels of connectivity with other sites. This would provide a pathway for re-colonisation and would increase the effective population size (N e ), protecting the population against the effects of genetic drift (Dailianis et al. 2011;Riesgo et al. 2016). However, high variance in reproductive success can occur in broadcast spawning marine invertebrates, reducing N e , and thus increasing vulnerability to bottlenecks (Hedgecock 1994).
Similarity in genetic diversity across all sampling sites implies that S. vesparium in Florida may still have sufficient genetic variation for resilience against some future stressors. In the Florida Keys, those include anthropogenic effects on water quality and global climate change (Wall et al. 2012;Kearney et al. 2015;Butler and Dolan 2017), as well as further cyanobacteria blooms. However, bottleneck signatures in some sites suggest that genetic diversity may have been previously lost due to unknown causes, and therefore caution should be exercised in management to prevent possible reductions in genetic variation.

Inbreeding and null alleles
Observed heterozygosity was lower in Florida and Barbuda than Belize and The Bahamas. However, all sites showed excesses in homozygosity, and departures from Hardy Weinberg equilibrium were present across loci and sites. This phenomenon can be caused by inbreeding, but it can also be attributed to the presence of null alleles, which were found in a number of loci and across all sampling sites in our study. Null alleles are caused by mutations in primer binding regions that prevent primers from binding, and subsequently cause amplification failure in PCR, either in both alleles (resulting in missing data), or for only one allele (resulting in false homozygotes). High null allele frequencies are commonly found in sponge microsatellite studies (Dailianis et al. 2011;Guardiola et al. 2012Guardiola et al. , 2016Chaves-Fonnegra et al. 2015;Pérez-Portela et al. 2015;Richards et al. 2016), suggesting that the problem may be common in the phylum, and is a known issue in other taxa (e.g., molluscs and insects; Chapuis and Estoup 2007). To reduce the impact of null allele bias on our estimates of inbreeding, we corrected F IS values for null alleles.
We found positive F IS values in all populations when corrected for null alleles, but this was not statistically significant in any of the sites when tested individually, potentially due to small sample sizes. Our genetic clustering analyses concluded that Florida was a single population, which enabled us to group the Florida sites for more statistical power, resulting in a significant positive mean F IS value. This suggests the presence of inbreeding in Florida S. vesparium populations, although the large 95% posterior probability intervals at the individual site level preclude a more fine-scale spatial assessment. Inbreeding has negative implications for fitness, thus our results highlight a potential concern for the health, reproductive success and longevity of S. vesparium in Florida.
Inbreeding is often characteristic of populations that have experienced declines. Hence, the positive F IS values we observed for S. vesparium in Florida may be due to mass mortality events, coupled with limited regional-scale connectivity to replenish the gene pool. However, high inbreeding coefficients are widespread in the Porifera (Guardiola et al. 2012;Bell et al. 2014;Chaves-Fonnegra et al. 2015;Pérez-Portela et al. 2015;Giles et al. 2015;Padua et al. 2017). This suggests that the positive F IS values for Porifera may be general characteristics of the phylum, perhaps associated with high philopatry due to limited larval dispersal. Additionally, sponges in the Clionaidae family can be simultaneously hermaphroditic (contain both eggs and sperm at the same time) (Piscitelli et al. 2011), so self-fertilisation, or selfing, could theoretically be possible; however there is currently no recorded evidence of selfing in the Porifera phylum. Positive F IS values can also be caused by excess homozygosity driven by Wahlund effects. These effects can occur when there is population structure within a site or group, and can be caused by reproductive asynchronicity or recruitment of different genetic cohort (Duran et al. 2004;Chaves-Fonnegra et al. 2015;Riesgo et al. 2016). With this in mind, it is difficult to fully gauge the implications of positive F IS for population health.

Genetic connectivity patterns
Spheciospongia vesparium exhibited strong population structure at the regional (Caribbean) spatial scale, indicating that connectivity among sponge populations in the four countries we sampled is low. These results are congruent with those of other sponge species, which exhibit high differentiation at large spatial scales in the Caribbean (López-Legentil and Pawlik 2009;Chaves-Fonnegra et al. 2015;de Bakker et al. 2016;Richards et al. 2016;DeBiasse et al. 2016), but also in other regions (Duran et al. 2004;Xavier et al. 2010;Pérez-Portela et al. 2015;Riesgo et al. 2016Riesgo et al. , 2019Brown et al. 2017;Padua et al. 2017;Taboada et al. 2018). Dispersal in marine species is affected by a number of factors and the complex interactions between them (Cowen et al. 2006;Cowen and Sponaugle 2009), including ocean current patterns and life history characteristics such as pelagic larval duration, larval behaviour, and reproductive strategies (Butler et al. 2011;Selkoe and Toonen 2011;Coelho and Lasker 2016). Although reproductive and larval traits for S. vesparium are not known, sponge larvae generally have short pelagic larval durations, limiting their dispersal capacity. This includes previously studied members of the Clionaidae family (Warburton 1966;Mariani et al. 2000Mariani et al. , 2001, to which S. vesparium belongs. Furthermore, Clionaidae larvae have been found to exhibit weak swimming ability, and commonly crawl (Mariani et al. 2000(Mariani et al. , 2001, further minimizing dispersal capacity. Our results on the population structure of S. vesparium are consistent with expectations of connectivity as determined by regional ocean current patterns. Our analyses indicated the presence of a barrier to gene flow between the Florida sites and Abaco in The Bahamas, which concurs with patterns found in genetic studies of other sponges (López-Legentil and Pawlik 2009;Richards et al. 2016;DeBiasse et al. 2016) and corals (Brazeau et al. 2005;Baums et al. 2010), and biophysical modelling predictions of fish and lobster larvae (Cowen et al. 2006;Truelove et al. 2017). This break is likely due to the strong Florida Current, which runs between The Bahamas and Florida, and can act as a strong barrier to dispersal.
Genetic differentiation was much larger between Florida and Abaco than between Florida and Belize, despite the geographic distance being much larger for the latter pair, as shown by genetic distance calculations and the DAPC analysis. Connectivity between Florida and Belize could be aided by the Caribbean Current and Loop Current, which can support larval transport from Belize towards Florida in as little as 7 to 10 days (Muhling et al. 2013). This is likely to be higher than the larval duration of most sponge species, but locations 'upstream' from Florida, such as the Yucatán Peninsula, could act as intermediate 'stepping stones' to aid gene flow between these areas, as appears the case for some marine diseases , thus reducing genetic differentiation.

Connectivity across the Florida Keys
Florida formed a single genetic cluster in our analyses (Geneland, Flock, DAPC), and the AMOVA showed that only 2.4% of the total genetic variation in the dataset was among sites in Florida. These results indicate that some level of connectivity is present across the Keys. We also found evidence of recent migration between Florida sites in our first generation migration analysis. According to genetic distance and PCoA analyses, sites such as Long Key (Bay-side) in the Middle Keys and Boca Chica Channel in the Lower Keys appeared well-connected to sites throughout the Florida Keys range. The complex currents found across the Florida Keys are likely to aid in connectivity among disparate sites. Although the main current dominating the area is the north-easterly running Florida Current, there are many local oceanographic processes that can affect larval dispersal patterns. Westerly running counter currents arise as a result of downwelling winds and offshore eddies and gyres (Lee and Williams 1999;Yeung et al. 2001;Kourafalou and Kang 2012), and eddies themselves also drive connectivity in the area (Sponaugle et al. 2005). Connectivity is also influenced by a species' life history. Reproduction of S. vesparium has not been described, however, oviparity occurs in some members of the Clionaidae family (Maldonado and Riesgo 2008;González-Rivero et al. 2013); if S. vesparium is also oviparous, additional dispersal of the gametes before fertilization may increase connectivity over longer distances compared to viviparous sponges. However, in other oviparous sponges, egg masses have been observed to stick to the substrate close to the mother sponge due to their envelopment in an adhesive material (Mariani et al. 2001). Furthermore, fertilization rates generally decrease over increasing gametic dispersal distances in broadcast spawners (Levitan 1991;Lauzon-Guay and Scheibling 2007). However, even a relatively small proportion of far-dispersing eggs that get successfully fertilized could increase the genetic connectivity between populations (Trakhtenbrotl et al. 2005).

3
Despite evidence of connectivity, there was still population structure among the Florida sites, demonstrating that the area does not form a completely panmictic population. Isolation by distance accounted for some of the structure across sites: genetic similarity decreases with geographic distance. This suggests that distance-limited dispersal influences population structure on smaller (< 160 km) spatial scales and is again likely to be due to the short pelagic larval duration found in sponges.
The Barrier analysis suggested a barrier to gene flow between adjacent Atlantic and Bay-side sites in the Middle Keys. The AMOVA confirmed that more genetic variation was found across the barrier than among sites on the same side of the barrier. In addition, the Florida-only PCoA showed the Atlantic sites (along with Waltz Key) separated from the rest of the sites by the second axis. These results suggest that dispersal through the channels between the islands of the Keys archipelago is limited, at least in the Middle Keys where we sampled. This is somewhat surprising considering the strong tidal flux through these channels (Smith 1994;Smith and Lee 2003), however, weakly-swimming larvae caught in the tidal flow may struggle to settle in areas close to the channels before being transported offshore or into the Bay. Furthermore, larval exchange could be limited by spatially and temporally variable inflow and outflow through the channels (Smith 1994;Yeung et al. 2001;Lee and Smith 2002). That being said, despite moderate statistical support for a barrier, the F ST value across the barrier was only 0.033, showing low genetic differentiation. Furthermore, they did not form separate populations in Geneland and FLOCK analyses. This indicates that although genetic differentiation is higher than would be expected due to distance alone, it is only a weak barrier to gene flow. More substantial population structure was found in the seagrass Syringodium filiforme between the Bay and Atlantic sides of the Keys in the same area (Bijak et al. 2018). This is likely due to vegetative propagation of S. filiforme through the sediments compared to larval propagation of S. vesparium through the water column.
The Barrier analysis also showed a putative barrier occurring between the Middle Keys and Lower Keys sites. However, this had low bootstrap support, and is likely to be an artefact of the isolation by distance pattern in the area, rather than a physical or oceanographic barrier to dispersal (Meirmans 2012).
Other patterns of population structure within Florida did not correlate with known physical or oceanographic features. For example, sponges near Waltz Key (a semi-isolated lagoon) were genetically different than those at many other sites in the Lower Keys, but not those in the Middle Keys. Although counterintuitive, this is not uncommon. Unexpected patterns of fine-scale genetic structure have also been observed in other sponges found along the Florida Keys reef tract (DeBiasse et al. 2010;Chaves-Fonnegra et al. 2015). Furthermore, a dispersal model based on water circulation patterns and larval characteristics did not accurately predict genetic connectivity patterns for A. cervicornis across the Florida reef tract, with the genetic data revealing more complex connections than the model predicted (Drury et al. 2018). Such patterns of chaotic genetic patchiness in the marine environment can be caused by 'sweepstakes reproductive success', the random survival of certain larval cohorts due to oceanographic conditions (Hedgecock 1982(Hedgecock , 1994Hedgecock and Pudovkin 2011). These effects can be found in species with high fecundity and high larval mortality. Alternatively, variable current regimes, as found in the Florida Keys (Lee et al. 1992), can result in temporally variable dispersal pathways. Both of these situations could lead to spatially heterogeneous genetic structure through genetic drift.

Sponge restoration implications
Our results have important implications for sponge restoration practice. Genetic diversity in S. vesparium is naturally high and, in addition, clonality is low, with only two identical genotypes found in our dataset. To maintain these genetic diversity levels, restoration should be carried out though the selection of genetically-diverse donor sponges. Donor sponges should not be extensively fragmented to produce a number of genetically identical transplants in a single location; instead, minimal fragmentation of many individuals and transplantation of whole sponges should be used. By maintaining high genetic diversity, restored populations can uphold evolutionary potential and resilience against future stressors, as well as avoid the negative fitness consequences of inbreeding. As our results indicate an absence of population clusters within the Keys, strong local adaptation does not appear to be present. This indicates that outbreeding depression is not a concern, and sourcing donor sponges does not have to be restricted to certain sites or environmental conditions.
Our findings highlight the importance of restoration work in Florida. Connectivity on the regional scale was low in our study, suggesting that immigration and gene flow into Florida may be limited. Populations in Cuba or the Gulf of Mexico may be more connected to Florida that the sites sampled here. However, patterns observed in this study suggest that migration is likely to be limited due to the oceanographic distances and limited pelagic larval duration. Active management on the local scale is therefore likely to be of vital importance to ensure that population numbers are maintained.
We show that connectivity is present over the range of the Keys, and we did not observe genetically isolated sites that would need to be prioritised for restoration action. However, our results also imply that connectivity in Florida is unpredictable, as we observed unexplained fine-scale structure. Furthermore, isolation by distance suggests that dispersal is distance-limited. These results show that natural repopulation of barren areas may be slow, especially if healthy populations are moderately distant. This may be compounded by the loss of acoustic larval recruitment cues in the area, itself caused by loss of sponge-associated endosymbionts (Butler et al. 2016). Given the crucial role of this important keystone species, sponge restoration is an important strategy in facilitating a more rapid return to ecosystem function following mass mortality events. However, this approach must be coupled with thorough investigation into the causes of the ongoing mass mortalities and ecosystem instability in the Florida Keys, and the implementation of measures to mitigate these issues. Furthermore, genetic diversity and its distribution among sites should be monitored regularly to ensure that genetic variation is maintained throughout the restoration program. This can now be accomplished relatively quickly using the molecular tools described in this study.