The paradox of retained genetic diversity of Hippocampus guttulatus in the face of demographic decline

Genetic diversity is the raw foundation for evolutionary potential. When genetic diversity is significantly reduced, the risk of extinction is heightened considerably. The long-snouted seahorse (Hippocampus guttulatus) is one of two seahorse species occurring in the North-East Atlantic. The population living in the Ria Formosa (South Portugal) declined dramatically between 2001 and 2008, prompting fears of greatly reduced genetic diversity and reduced effective population size, hallmarks of a genetic bottleneck. This study tests these hypotheses using samples from eight microsatellite loci taken from 2001 and 2013, on either side of the 2008 decline. The data suggest that the population has not lost its genetic diversity, and a genetic bottleneck was not detectable. However, overall relatedness increased between 2001 to 2013, leading to questions of future inbreeding. The effective population size has seemingly increased close to the threshold necessary for the population to retain its evolutionary potential, but whether these results have been affected by sample size is not clear. Several explanations are discussed for these unexpected results, such as gene flow, local decline due to dispersal to other areas of the Ria Formosa, and the potential that the duration of the demographic decline too short to record changes in the genetic diversity. Given the results presented here and recent evidence of a second population decline, the precise estimation of both gene flow and effective population size via more extensive genetic screening will be critical to effective population management.

The loss of genetic diversity is a primary concern in conservation biology and should be at the forefront of strategy to promote biodiversity conservation management 1,2 . It was once thought that marine fishes were too abundant and fecund to be faced with threats from overexploitation 3 . Fishing was a right, oceans were open access, and baseline data were, more or less, non-existent. More than a century of literature and observations later, it is now recognized that the paradigm of inexhaustible marine fishes and exploitation is unequivocally false [4][5][6] . Unfortunately, while the burden of proof often falls to scientists to demonstrate that anthropogenic activity is the root cause of environmental degradation, the lack of baseline data for comparison often hampers inferences on impact 7 . As a result, many marine environments have been subjected to degradation in the form of pollution, overexploitation and habitat loss. While the paradigm of inexhaustible ocean resources has eroded, many ongoing discussions continue in the field of conservation biology. One of which, relevant here, was proposed by Lande 8 , who suggested that species are often driven to extinction before genetic consequences have time to take effect. However, a bulk of empirical evidence indicates the contrary 9 , and more recently, reduced genetic diversity has been documented in threatened species 10 . Quantifying gene flow, metapopulation structure and demographic and stochastic events are necessary when assessing whether a reduction in population size is likely to lead to a reduction in genetic diversity 11 .
Interestingly, just as there has been contention about the effects of genetic changes on threatened populations, a similar debate has been held over inbreeding in the wild. Despite many questioning the likelihood that inbreeding depression can lead to declines in wild populations 12-14 that ca. 90% of studies have shown reduced fitness in inbred populations compared with non-inbred individuals 15 . For example, inbreeding depression was reported in 99 species of birds in which failed hatching rate increased with genetically similar parents 16 , as well as in fish, where inbred populations of Poeciliopsis monacha suffered from spinal curvature, deformation and reduced resistance to low oxygen 17   DNA extraction, amplification, fragment analysis, genotyping and pre-processing. DNA from filament samples was extracted by standard phenol-chloroform protocol, according to Sambrook 42 . Following DNA extraction, amplification of 12 highly polymorphic microsatellite loci (Electronic supplement S1) isolated for H. guttulatus (or obtained by cross-amplification in H. hippocampus) was performed as described in Pardo et al. 43 and van de Vliet et al. 44 . PCR reactions with labelled primers were performed using standard procedures, and amplified products were run on an ABI 3130XL (Applied Biosystems) automated sequencer. Genotype scoring was performed using STRand (https:// vgl. ucdav is. edu/ infor matics/ strand. php). It was assumed that a single panmictic population was being assessed 31 , and downstream analyses of the 2013 sample pooled all locations into a single year population. To filter for missing data, the missingno function of the poppr 45 R-package 46 was used, removing loci with more than 25%, and individuals with more than 10% missing data, respectively. Genotyping errors and null alleles were estimated using calculations by Summers  Relatedness analysis. Relationships between individuals from the two sampling years were estimated with the program ML-RELATE 54 . ML-RELATE calculates coefficients of relatedness (r) and putative relationships among individuals (e.g., unrelated, siblings, parent/offspring) using a maximum likelihood approach. Random genotype simulations were performed 1000 times for likelihood ratio tests and identified plausible relationships (full-sibling, half-sibling, parent-offspring, unrelated) based on a 99% confidence interval criterion 54 . Estimates of the percentage of sample pairs identified as full-siblings and half-siblings were calculated. ML-RELATE can accommodate the presence of null alleles in the estimates via a Monte-Carlo randomization test 55 and the U-statistic of Rousset 56 . The test is one-tailed, i.e. it estimates the probability of obtaining the observed U statistic or greater value under Hardy-Weinberg conditions. If null alleles are present, ML-RELATE will use a maximum likelihood estimate of their frequency in all calculations which is considered to be more accurate than other estimators 47,49,50,54 . Inbreeding. Estimates of inbreeding depression often rely on limited datasets that can be affected by the presence of null alleles, which have been shown to induce an upward bias in inbreeding coefficient (F is ) estimates 57 . This is especially true when sampling in small geographical areas in which inbred or closely related individuals occur 58,59 , as is the case for H. guttulatus from Ria Formosa. Therefore, in this study, simultaneous estimation of the inbreeding coefficient, null allele frequencies and random genotyping failures were conducted using the software INEST v2.2 59 . This software implements a Bayesian approach comparing different combinations of parameters (f = inbreeding, n = null alleles, b = genotype failures) and chooses the model with the lowest deviance information criterion (DIC) as the best fit. MCMC (Monte Carlo Markov Chain) used 500,000 cycles and a burn-in of 50,000, retaining parameters every 100th cycle.
Genetic demography. To explore the possibility of a demographic bottleneck in Ria Formosa seahorses, two tests were performed. First, we used a test for heterozygosity excess 60 implemented in INEST v2.2 59 , under the null expectation that a population under a mutation-drift equilibrium should have an equal probability of exhibiting a heterozygosity excess or deficiency at a given locus. Populations that have experienced a recent reduction in size typically exhibit fewer alleles and lower heterozygosity. However, allelic diversity declines at a faster rate than heterozygosity due to random genetic drift, which eliminates lower frequency alleles. Thus, this test contrasts the expected heterozygosity assuming Hardy-Weinberg equilibrium with expected heterozygosity under mutation-drift. The latter is more sensitive to declines in low-frequency alleles that are expected when populations become reduced in size. A bottleneck is inferred when heterozygosity under Hardy-Weinberg equilibrium exceeds heterozygosity under mutation-drift 60,61 . This analysis was run using a two-phase mutation model (TPM) with the proportions of multi-step mutations (p g ) set as 0.22 and the mean size of multi-step mutations (D g ) as 3.1 as recommended by Peery 61 . Stepwise mutation (SMM) and infinite allele models (IAM) were www.nature.com/scientificreports/ also carried out for comparison, as these represent the extremes of mutation models. The Wilcoxon signed-rank test was applied to test for the significance of heterozygosity excess for both mutation models using 10 6 permutations to approximate the exact value, as well as combined z-scores 62 . Second, a mode shift test was used to determine whether a shift has occurred in allele frequency classes. This test assumes that stable populations have a peak in allele numbers at the lowest frequency class, displaying an L-shaped distribution. Population bottlenecks tend to distort this distribution to the right, as the loss of alleles results in a rightward skew in the allelic distribution 62 . To test if the data conformed to the L-shaped null hypothesis of a stable population, Bottleneck v2.2.01 63 was used.
Effective population size. Estimates of contemporary effective population size (N e ) may be obtained from single or temporal samples 64 . In this study, the former approach was used, because the samples from 2001 and 2013 covered a restricted time interval and there were no age estimates from the sampled individuals 65,66 . N e was assessed from levels of linkage-disequilibrium (LD) 67 with NeEstimator v2.1 68 . Monogamy was assumed 69 and 95% confidence intervals of Ne were estimated using the jackknife-across samples method 70 . Sample sizes are an important consideration when estimating N e, and small sample sizes can result in considerable biases 66 . Due to the sample size disparity between the 2001 and 2013 samples, genetic data for both years were rarefied by generating 1000 replicate files with 30 individuals in each. Alleles with low frequencies (< 0.02) were omitted to avoid bias because their inclusion can bias estimates of N e 68 . Population-specific N e was calculated by taking the harmonic mean of LDNe estimates from each year the population was sampled. As 95% confidence intervals of Ne based on the LD often included infinity, only 95% lower bounds are reported, as these are considered valuable indicators of changes in population size 64 .

Results
A total of 434 individuals were screened across 12 loci (https:// tinyu rl. com/ ya5tg k6a; Annex B). Four loci (Hgut9, USC1, USC6 and USC7) and 135 genotypes were removed after filtering at 25% and 10% of missing data, respectively, resulting in 291 informative individuals across eight loci, 33 (Table 1), and average PAr increased from 4.1 to 6.5 during this period (Table 1), despite a decrease at the Hgut6 locus from 8 to 4. The average inbreeding coefficient across loci declined from 0.16 in 2001 to 0.05 in 2013 (Table 1).     (Fig. 4a,b) for both tested methods. Both approaches (single point estimate with true sample size values, and point estimates using randomly generated files with rarefied sample sizes) provided larger estimates for the 2013 sample. However, the 95% confidence intervals for population-specific N e values overlapped, and upper confidence limits were indeterminate (Fig. 4a).

Discussion
The initial hypothesis for this investigation was based on the premise that the drastic decline in census population numbers reported between 2001 and 2008 would be reflected in reduced genetic diversity in subsequent generations. Based on this prediction, lower genetic diversity was expected in the 2013 sample. In contrast to this expectation, the H. guttulatus population of the Ria Formosa showed an increasing level of genetic diversity between 2001 and 2013. Only the number of private alleles at one microsatellite locus (Hgut6) showed any substantial decrease. Over this same period, although N e has risen and F is has seemingly declined, no evidence was found to support the occurrence of a genetic bottleneck, apart from the observation that the overall relatedness has increased. Before discussing the implications of these findings, it is important to highlight some of the caveats of this study. Firstly, it must be stressed that the sample size disparity between both years may have, to some extent, influenced the results presented here. Some diversity estimates are more prone to biases due to sample size than www.nature.com/scientificreports/ others, making them more useful in the context of a comparison between years. For instance, while the number of alleles per locus in a small sample is greatly biased relative to an estimate based on a large sample size, this can be overcome with rarefaction, which improves the ability to compare between samples differing in size 71 . Secondly, the genetic diversity of the 2001 sample was derived from a higher number of sites (15) than that sampled in 2013 (6) 26,36 . Thirdly, the presence of null alleles in the dataset may have affected the study, particularly as cross-species amplification was used 72 . While the frequency of null alleles in the 2001 dataset was higher than that inferred for 2013, likely due to poor preservation, no evidence of null alleles in previous studies has been detected using these markers 43,44 . Thus, we found no reason to alter allele frequencies or exclude loci, except when using INEST, in cases in which the inclusion of null alleles or genotyping errors was found to upwardly bias F is estimates. Finally, it is important to note that neutral microsatellites, particularly a limited number of loci, do not represent genome-wide diversity, and may fail to detect changes in the genetic diversity of loci experiencing selection 73 .
Genetic diversity. Genetic diversity indices were generally high and showed an overall increase between 2001 and 2013, rejecting our initial hypothesis of a reduction in genetic diversity due to census population declines over the study period.  74 . Although these differences may be due in part to the use of different markers, Hgut4 was used in both studies, and as both studies employed similarly large sample sizes, this is unlikely to be a sampling artefact. When all available demographic and genetic diversity results are taken into consideration (Table 3), evidence for an association between population reductions based on census data and reductions in genetic diversity as estimated by microsatellites is equivocal.
In other similar lagoonal environments in the Mediterranean, Taranto and Thau, expected heterozygosities at eight microsatellites loci are very low (N: 40, H e : 0.31 and N: 8, He: 0.39, respectively) 34 . Other coastal locations in the Mediterranean typically display low expected heterozygosities ranging from 0.30 to 0.54 with sample sizes from 2 to 17 see Table 3 23 . A possible scenario is that gene flow into the Ria Formosa has enabled some level of genetic diversity to permeate back into the population, a pattern supported by a more extensive temporal survey of seahorses in the area (Wilson et al., unpublished data). This hypothesis is also supported by studies that have shown that South Iberian populations display low levels of microsatellite genetic differentiation 74 and are panmictic with respect to 286 SNPs 31 . Despite suggestions that seahorses are capable of assisted long-distance translocations, likely through rafting on floating debris, carried by oceanic and coastal currents 75 , that this is implausible in the Ria Formosa, where H. guttulatus have high site-fidelity and very limited dispersal capacity 38 , and the lack of intermediate required habitat does not facilitate seahorse dispersal. The panmictic nature of the South Iberian populations remains to be fully explained even though a commonly cited rule of thumb is that one migrant per generation is needed to sufficiently minimise the loss of polymorphism and heterozygosity within subpopulations 18,76,77 . Relatedness and inbreeding. While the present results are inconsistent with a genetic bottleneck related to the inferred population decline, heterozygote deficiency was noted at six out of the eight loci, which could be a result of recent population expansion 60 , corroborating other findings 36,37 . It is possible that the doubling of related individuals between 2001 and 2013 is an artefact of the demographic decline inferred in 2008. With fewer individuals in the period after the population reduction, the likelihood of sampling related individuals would be expected to be higher, due to the remaining individuals contributing to the recovery reported by 2012. Neverthe- Table 3. Summary results from previously published genetic and demographic data from Ria Formosa. N number of individuals analysed, genetic diversity: observed heterozygosity, for microsatellite data, and haplotype diversity for mtDNA data. www.nature.com/scientificreports/ less, both our investigation and Woodall's 74 study used the heterozygosity excess method 60 , which is expected to have more power to detect recent population bottlenecks relative to the M-ratio statistic 78 . With this in mind, we suggest that H. guttulatus in the Ria Formosa may be naturally prone to population fluctuations that may be exacerbated by anthropogenic activity, e.g. clam farming, fisheries, dredging 79,80 . Seahorse populations are vulnerable to fluctuations [81][82][83] , which makes the conservation of their preferred habitats critical 26,41 and of great urgency if disturbances and habitat degradation are identified.
Effective population size. Both 2001 and 2013 N e point estimates were negative, indicating that the results can be entirely explained by sampling error rather than drift 68 . This situation usually occurs when N e is large, which may suggest that the Ria Formosa population was substantial during both sampling periods. Linkage disequilibrium due to drift is low in large populations and difficult to detect, rendering the estimation of N e in large populations (i.e., N e > 1,000) challenging, particularly with the small number of markers employed in this study 64 . More comprehensive genetic screening via high-throughput sequencing may offer an opportunity to overcome this limitation, and more precisely estimate N e 64 . Reporting the lower limits of these estimates remains valuable as an indicator of changes in population size 64 . Here, lower bound N e means remained relatively stable, with a slight increase, less than the increases inferred from harmonic means and point estimates. The interpretation of these N e estimates must be made with caution 84 because there is only weak genetic differentiation between the Ria Formosa and other South Iberian subpopulations 31,74 . This implies a level of historical and/or contemporary gene flow, which seems improbable given the species' life-history traits and habitat preferences. Nevertheless, N e estimates may be inflated and/or reflect Ne of a larger metapopulation rather than the local population. While this cannot be ruled out, even with an apparent increase in N e , the estimates of harmonic means reported here centre around 2000 individuals, and 50% of the lower bounds are below 500 (Fig. 4b). A commonly used guideline employed in a conservation context is the 50/500 N e rule of thumb 85 , although more recently, this has been revised upward to 100/1000 N e 86,87 . N e ≥ 100 thought to be required to avoid inbreeding depression and limit the loss of total fitness to 10% and N e ≥ 1000 is suggested in order to retain evolutionary potential in perpetuity. Our lower bound N e estimates are worryingly close to levels thought to be required for the maintenance of long-term evolutionary potential (500/1000 N e ), particularly in light of unpublished data that has shown that another severe decline has taken place, indicating an overall decline over the last 20 years.
Overall, based on the results presented here and in previous literature, it seems that the H. guttulatus population of the Ria Formosa has been able to retain relatively high genetic diversity despite evidence of a recent demographic decline. Perhaps the decline reported in 2008 was not sustained for long enough for drift to have had an effect and a sufficient number of individuals were left to retain high genetic diversity. According to empirical observations e.g. kit foxes 88 and simulations 89 , changes in genetic diversity following a decline in population size may take a number of generations to become obvious, and that the short time span of our study was not enough to capture the full range of genetic parameter fluctuations.
Future studies using an increased number of loci compared to those used here (such as Single Nucleotide Polymorphisms, SNPs) could help to test this prediction. Alternatively, a recently published survey of available seahorse microsatellites identified 18 highly polymorphic loci for H. guttulatus, which would provide a more extensive panel of markers for the analysis of neutral diversity 90 . If SNPs or more microsatellite loci reveal declines in genome-wide genetic diversity, F is will become extremely critical, as small differences in F is values can significantly alter extinction risk 18,91 . Furthermore, the estimation of N e for the broader region may highlight if the N e value reported here is inflated through gene flow. Considering the lower bound estimates of N e , the precautionary principle must be taken. We recommend that active management actions for this population should be continued and enhanced to protect the diversity and numbers indicated here.
Coastal lagoons such Ria Formosa rank among the most productive ecosystems on Earth, and provide a wide range of ecosystem services and resources. Anthropogenic impacts are escalating in many coastal lagoons worldwide because of increasing population growth and associated land-use alteration in adjoining coastal watersheds 92 . Among the main stressors affecting coastal lagoons are habitat loss and alteration; eutrophication; sewage and organic wastes; fisheries overexploitation; sea-level rise; chemical contaminants, sediment input/ turbidity and floatables and debris 93 . Despite Ria Formosa being a semi-protected lagoon, many activities that alter the environment are permitted, including legal fisheries, anchoring and dredging, which contribute to habitat loss in the lagoon 79,94 . Habitat loss was identified as the leading probable cause for seahorse decline 38 , which was to some extent confirmed by other studies 36,37 . The observed extreme fluctuations in seahorse populations are correlated to the availability of holdfasts 95 . However, we cannot discard the impact that other undetermined abiotic or biotic factors such as warmer temperatures, food availability, or abundance of predators 36 or overfishing may have on local populations. There is empirical evidence of a rapid and somewhat uncontrollable unreported and unregulated fishing (IUU) activity which may be a significant factor contributing to population decline (J. Palma, person.comm.).
It is difficult to explain the maintenance of the high levels of genetic diversity observed in Ria Formosa individuals when census data suggest such a dramatic decline. While we cannot discard the possibility that the individuals have moved to other areas, the lagoon has been visited regularly by divers experienced in seahorse detection. During 2020, the number of visited sites was increased, and although an overall decline in population size was observed, local increases were found at the microhabitat level. Expansions of individual home ranges in response to low densities of potential mates have been observed in previous years 35 , suggesting that short distance migrations outside the survey area may have influenced census results. Seahorses have been found to colonize previously degraded habitats after being enriched with artificial structures 95 . We discard the possibility that the seahorses might have moved outside Ria Formosa in search of partners, because of the lack of suitable habitats (e.g., availability of holdfasts), both to the west and the east of the lagoon. It is perplexing not to find a genetic

Conservation implications.
This study showed that genetic diversity and effective population size of the Ria Formosa long-snouted seahorse was retained through a severe recent demographic decline as inferred from census data. However, this population has suffered an even steeper decline in recent years (2018-2020) (J. Palma, pers. comm.). Therefore, despite the genetic results, conservation measures and demographic and genetic monitoring should continue to be implemented. Enhanced protection and restoration of H. guttulatus habitats would help to retain the genetic diversity in the Ria Formosa. The use of artificial holdfasts has been shown to be a useful tool for repopulating areas in which habitats have been damaged or disappeared 95 . However, this does not resolve the issue of degrading habitat, which should be at the forefront of conservation in the Ria Formosa. Clam farming poses a significant threat to the seagrass Zostera noltii 80 , and other species of seagrass have been degraded by coastal construction and dredging 79 . H. guttulatus have a preference for complex habitats 26,35,41 , which are likely to be subject to human impacts, fluctuations and disturbances in the dynamic barrier island system that characterizes the Ria Formosa 97 . Identifying temporally stable areas and enhancing habitat may help to ensure the long-term viability of seahorses in the Ria Formosa 98 , in addition to the enhanced protection of important refuges. Finally, future studies will need to clarify the potential impact of gene flow between subpopulations of the South Iberian Peninsula, as previously suggested by other authors 31,74 in order to determine the relevant conservation unit for seahorses in the region.

Conclusion
Molecular data suggest that, contrary to our initial hypothesis, a severe demographic decline reported in 2008 did not reduce the genetic diversity of H. guttulatus in the Ria Formosa, nor cause a genetic bottleneck, and that the effective size of this population has instead modestly increased. A more in-depth genomic study could help to illuminate genetic changes at the whole genome level and identify genetic regions experiencing differential selection, which may be especially prone to population fluctuations. Finally, although the goal of this study was to assess genetic changes in diversity, it must be acknowledged that the effective conservation of essential habitat is critically important 99,100 . Bearing in mind the observed pattern of population fluctuations, the use of artificial holdfasts may alleviate pressures from environmental variation. Active and enhanced monitoring and conservation of key habitats should remain at the forefront of conservation efforts in nearshore environments such as the Ria Formosa.