Unexpectedly high genetic diversity in a rare and endangered seabird in the Hawaiian Archipelago

Seabirds in the order of Procellariiformes have one of the highest proportions of threatened species of any avian order. Species undergoing recovery may be predicted to have a genetic signature of a bottleneck, low genetic diversity, or higher rates of inbreeding. The Hawaiian Band-rumped Storm Petrel (‘Akē‘akē; Hydrobates castro), a long-lived philopatric seabird, suffered massive population declines resulting in its listing under the Endangered Species Act in 2016 as federally Endangered. We used high-throughput sequencing to assess patterns of genetic diversity and potential for inbreeding in remaining populations in the Hawaiian Islands. We compared a total of 24 individuals, including both historical and modern samples, collected from breeding colonies or downed individuals found on the islands of Kaua‘i, O‘ahu, Maui, and the Big Island of Hawai‘i. Genetic analyses revealed little differentiation between breeding colonies on Kaua‘i and the Big Island colonies. Although small sample sizes limit inferences regarding other island colonies, downed individuals from O‘ahu and Maui did not assign to known breeding colonies, suggesting the existence of an additional distinct breeding population. The maintenance of genetic diversity in future generations is an important consideration for conservation management. This study provides a baseline of population structure for the remaining nesting colonies that could inform potential translocations of the Endangered H. castro.


INTRODUCTION
Many seabird species in the order Procellariiformes have experienced historical bottlenecks due to vulnerability to anthropogenic disturbances (Weimerskirch, Brothers & Jouventin, 1997;Milot et al., 2007;Paleczny et al., 2015). Current management, including predator control, laws protecting seabirds, and habitat restoration are helping to stabilize and restore some of these seabird populations (Young et al., 2009;Jones & Kress, 2012;Young et al., 2013). However, despite this active management, many recovering populations still have relatively low genetic variation (Barrowclough, Corbin & Zink, 1981; pairs suspected to remain (Pyle & Pyle, 2017), there is concern that the Hawaiian populations may have problems normally associated with small population size, such as demographic stochasticity and inbreeding (Brook et al., 2002;Frankham, 2005;Kennedy, 2009). Furthermore, as a member of the Procellariiformes order they are likely highly philopatric (Milot et al., 2007;Jones & Kress, 2012;Antaky et al., in press), with a high rate of return to the natal colony for breeding, potentially contributing to limited gene flow.
Given these pressing management concerns, as a first step in assessing the vulnerability of the remaining populations to a changing environment, in this study we evaluated patterns of genetic variation in Band-rumped Storm Petrels nesting in the Hawaiian Islands. High-throughput sequencing delivers high yields of genetic data across the genome and allows for accurate estimates of genetic relatedness, which is particularly useful for small, potentially inbred populations of non-model organisms. Using high-throughput sequencing, we assessed inbreeding, genetic diversity, genetic structure, and demographic history.

Sample collection
We obtained 18 H. castro samples from across the Hawaiian Islands from 4 years of field effort to add to six existing museum specimens, for a total sample size of 24 individuals. Colony-sourced samples came from birds found near suspected or known breeding areas on Kaua'i, Maui, and the Big Island, while non-colony sourced samples came from the island of O'ahu ( Fig. 1

Laboratory analyses
We individually extracted DNA from the blood and feather samples using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's protocol. We quantified the extracted DNA with the AccuClear TM Ultra High Sensitivity dsDNA Quantitation Kit (Biotium, Hayward, CA, USA). Due to low DNA yield, we performed whole genome amplification on individual samples with the REPLI-g UltraFast Mini-kit (Qiagen, Valencia, CA, USA) which effectively and accurately increases yields of high-fidelity DNA (Ahsanuddin et al., 2017). We prepared the replicated whole genomic DNA from all 24 individuals for reduced representation genomic sequencing using the ezRAD protocol version 3.2 (Toonen et al., 2013;Knapp et al., 2016).

Genetic data analyses
We used the dDocent pipeline (Puritz, Hollenbeck & Gold, 2014) to assemble loci and call single nucleotide polymorphisms (SNPs) within the aligned sequences. In dDocent the following settings were used: 90% similarity to cluster reads, match score of one, mismatch score of four, gap penalty of six, minimum coverage of four within individuals, and minimum coverage of three between individuals. We filtered the resulting Variant Call Format (VCF) file using vcftools (Danecek et al., 2011), retaining 13,708 SNPs found in 90% of all individuals with a minimum quality value of 30, and with 20-200× read coverage. We separated the VCF into mitochondrial SNPs (mtDNA) and nuclear SNPs (nDNA) by aligning sequences to the H. castro mitochondrial genome (Antaky et al., 2019) using bwa-mem (Li, 2013), samtools (Li et al., 2009), and vcftools (Danecek et al., 2011). We calculated fixation indices (F st ) for both nDNA and mtDNA by using the Weir & Cockerham (1984) unbiased calculation in vcftools (Danecek et al., 2011). We analyzed the nDNA and mtDNA separately in R (R Core Team, 2013) using the package 'PCAdapt' to run a principal components analysis. We transformed the trimmed alignments in PGDSpider (Lischer & Excoffier, 2012) and ran them in STRUCTURE (Pritchard, Stephens & Donnelly, 2000) to identify the likely number of populations from which the samples came and infer proportion of ancestry for each individual. We ran STRUCTURE for one, two, three, and four populations (K), with 10 iterations for each K at a burn-in period of 10,000 steps and 10,000 steps after burn-in and 10 iterations at a burn-in period of 100,000 steps and 100,000 steps after burn-in. We performed multiple runs to obtain better estimates of the posterior probability of each K value. We input results from both STRUCTURE runs into the program STRUCTURE HARVESTER (Earl & VonHoldt, 2012), to calculate the ad-hoc statistic (ΔK) suggested by Evanno, Regnaut & Goudet (2005) that takes into account the change in the log probability of the data between increasing numbers of clusters. We used TASSEL (Bradbury et al., 2007) to determine the nucleotide diversity (π), Watterson estimator of diversity (θ), and Tajima's D (D T ). For population statistics, individuals were grouped by island from which they were sourced from. We constructed a Bayesian Skyline Plot, a coalescent-based graphical method, in BEAST v.2.5.0 (Bouckaert et al., 2014) and TRACER v.1.6 (Rambaut et al., 2016) using mtDNA from all individuals to infer potential historical fluctuations in effective population size (N e ). The Bayesian Skyline Plot framework takes into account genealogy, demographic history, and substitution-model parameters in a single analysis (Ho & Shapiro, 2011). We ran the Bayesian Skyline Plot analysis using strict clock models as it is considered a good approximation for intra-population level analyses and implemented in other studies on seabird evolution (Younger et al., 2015;Iglesias-Vasquez et al., 2017). We used the molecular clock rate of 1.94 × 10 −8 substitutions/site/year which was calibrated from Hydrobatidae (Storm Petrel) mtDNA using only Hydrobates spp. (Weir & Schluter, 2008).

RESULTS
From our 24 individuals (nine samples from Kaua'i, two from O'ahu, two from Maui, and 11 from the Big Island) we obtained a total of 650,048,040 sequences. We calculated genetic diversity statistics using 13,708 genomic SNPs shared among all individuals. The average nucleotide diversity (π) was similar across individuals found on Kaua'i (0.159), O'ahu (0.197), Maui (0.172), and Big Island (0.176). The Watterson estimator of expected genetic diversity at equilibrium (θ) was higher for Maui (0.376) and O'ahu (0.333) with the caveat of low sample size, but lower for Kaua'i (0.193) and the Big Island (0.225) where the sample size was higher. Tajima's D (D T ) was negative for all island colonies ( Table 2). The inbreeding coefficient (F IS ) ranged from −0.354 to 0.097 across islands.
We calculated pairwise F st values between island-scale groupings of individuals using both nDNA and mtDNA. We identified 13,641 nuclear and 67 mitochondrial shared SNPs among all individuals. Based on nDNA among islands, the highest differentiation was found between two groupings of islands (Supplemental Material S1): individuals on Maui differed from those on Kaua'i and Big Island, and individuals on O'ahu differed from those on Kaua'i and Big Island. The lowest F st value was between the Maui and O'ahu individuals, although these locations were only represented by two individuals and must be interpreted cautiously. F st values based on mtDNA among islands followed a similar pattern, showing differentiation between the same two groupings of islands (Supplemental Material S1). The Maui and O'ahu island samples showed no differentiation within nDNA (F st = 0), and therefore we performed the same analyses combining the Maui and O'ahu individuals into a single population, which did not produce qualitatively different results or alter conclusions (Supplemental Material S2). We investigated population structure using both nDNA and mtDNA with a Bayesian clustering approach. Using the criteria of Evanno, Regnaut & Goudet (2005), the grouping K = 2 received the highest support for nDNA and the grouping K = 3 received the highest support for mtDNA (Supplemental Material S3). Although the Evanno criteria cannot calculate the likelihood for K = 1, due to the separation of the Maui/O'ahu group and Big Island/Kaua'i in the K = 2 structure plot and PCA plot, it is unlikely that all individuals tested fall into a single genetic population. Furthermore, the Maui/O'ahu group had a 0.977 average probability of assigning together while the Big Island/Kaua'i group showed a 0.904 average probability of assigning to the same cluster. The structure analysis from the nDNA shows some support for separation of populations by island ( Fig. 2A), which is consistent with philopatry of breeding individuals. We evaluated relationships among individuals using a principal component analysis (PCA). The PCA based on nDNA resulted in at least two groupings, with some indication of separation by island, the same pattern as the structure analysis (Fig. 3A). The PCA and structure analysis based on mtDNA did not show the same separation of population by island as seen in the nDNA (Figs. 2D and 3B).
Reconstruction of the population size history by means of the coalescent Bayesian skyline plot using mtDNA data suggests a likely continuous population decline over the last 500 years for H. castro with indication of a small expansion in the last 25 years (Supplemental Material S4). Although the general trend in the skyline plot shows a decline, the median effective population number stays relatively stable with broad confidence intervals indicating that a constant population size over the whole time period is also possible. The current effective population size (N e ) ranged from 166 to 2327 individuals with a mean of 414. As the Bayesian skyline plot assumes that there is panmixia within the population tested, the four individuals found on Maui and O'ahu were removed from this analysis as they may belong to a different genetic group. Furthermore, results from this analysis should be interpreted with caution as the SNP data may differ from the overall mutation rate for mtDNA, as the more slowly evolving sites within mtDNA are likely not fully represented.

DISCUSSION
In this study, we found that H. castro in the Hawaiian Islands had relatively low inbreeding estimates and high genetic diversity, despite a relatively small population size and an Plots are based on analysis of 13,641 nuclear SNPs for K = 2 (A) and K = 3 (B), and 67 mitochondrial SNPs for K = 2 (C) and K = 3 (D) using STRUCTURE. Each bar represents an individual bird and the color represents the assignment probability to a particular genetic group. The most likely clustering for K, denoted with an asterisk, was determined by STRUCTURE HARVESTER.
Full-size  DOI: 10.7717/peerj.8463/ fig-2 assumed high degree of philopatry based on taxonomic order (Antaky et al., in press). Individuals from the presumed breeding colonies on the Big Island and Kaua'i show little differentiation, but individuals recovered from Maui and O'ahu do not assign to breeding colonies on either the Big Island or Kaua'i island, suggesting the presence of another distinct population in the region. Analysis of Hawaiian populations of H. castro using SNP data indicates an excess of rare alleles (mean D T = −3.425 ± 2.76), low rates of inbreeding (mean F IS = −0.164 ± 0.221), and high nucleotide diversity (mean π = 0.176 ± 0.016). Although the species has undergone a decline in population size, there is no evidence of inbreeding, potentially due to promiscuous mating behavior or sex-biased dispersal (Greenwood, 1980;Amos et al., 2001;Milot et al., 2007). Nucleotide diversities of H. castro in the Hawaiian Islands were higher than those found in studies of some non-endangered seabird species using RADseq methods (Dierickx et al., 2015;Tigano et al., 2017) but not all (Clark, 2018). The Bayesian skyline plot suggests a steady population decline after the introduction of mammalian predators to the Hawaiian Islands, with a small expansion in the last 25 years, but due to small sample size and assumptions of the Bayesian skyline analysis, evidence of recent population size change is relatively weak (Supplemental Material S4). Past studies estimated the breeding population of H. castro to be 660 individuals based on radar surveys, observational records, and song meter data (Pyle & Pyle, 2017), which falls within the 95% confidence intervals for effective population size found in this study 559). While these results are consistent, comparisons of census and effective population size are notoriously difficult (Frankham, 1995;Turner, Wares & Gold, 2002), and further studies investigating population size are necessary as skyline plots are not appropriate for the inference of complete historical demography (Grant, 2015).
Due to the nature of cryptic Endangered species with low population numbers, there are potential introduced biases from small sample size and opportunistic sampling. We acknowledge that with reduced sampling, we risk not capturing the complete genetic diversity in the populations sampled. Thus, additional studies will be useful to validate the results suggested from this study. Whole genome amplification was used to increase DNA yield in samples. Although whole genome amplified DNA for next-generation sequencing and array applications have been debated, studies using Qiagen's REPLI-g kits to amplify DNA yield have proven to be effective in limiting amplification bias across alleles and producing comparable results to the non-amplified sequence (Cichon et al., 2008;Ahsanuddin et al., 2017). Due to inaccessibility and lack of knowledge of nest site locations, sampling cannot be performed on confirmed breeding adults in the field. It is possible that carcasses found at or near colonies may be juveniles or adults visiting from other islands or populations. To limit this bias, we sequenced each individual and performed analyses that do not assume population assignment (PCA, STRUCTURE). However, for population diversity statistics and F st analysis (Table 2; Supplemental Material S1), which rely on population allele frequencies, we note that results should be interpreted with caution due to small sample sizes.
Based on the structure analysis and PCA, we found evidence for at least two distinct groups, with individuals from Kaua'i and the Big Island grouping together and individuals from Maui and O'ahu islands not assigning to that same population ( Figs. 2A and 3A). These genetic patterns do not match the island chain's geography, and instead may be due to genetic drift in small fragmented colonies, ocean regime around the islands (Friesen, 2015), or that the individuals collected on Maui and O'ahu were visiting birds that may belong to another distinct breeding population outside Hawai'i. In sum, the Maui and O'ahu island individuals differ from the individuals breeding on Kaua'i and the Big Island. This is further supported by the same analyses performed with Maui and O'ahu sourced individuals combined together (Supplemental Material S2). Due to small sample size, however, our results should be interpreted with caution. The PCA and structure analysis based on mtDNA did not show the same pattern of separation as nDNA, possibly because mtDNA does not account for male-mediated dispersal, or because the mtDNA dataset had fewer SNPs included in the analysis.
With some indication of differences among islands in the structure analysis and PCA plots, genetic data are consistent with the expectation that H. castro is a highly philopatric species. While Procellariiformes are more likely to return to their natal colony than disperse, most species are not completely philopatric, and a small percentage of individuals are likely to disperse to new colonies (Antaky et al., in press). It only takes a small amount of dispersal, that is, less than ten migrants per generation, to homogenize genetic structure (Mills & Allendorf, 1996). Seabirds also mate while on visiting forays to neighboring colonies, increasing gene flow beyond that expected based on dispersal from the natal colony (Young, 2010). The individual sampled in 1994 on the Big Island (Table 1), which clustered with the Maui/O'ahu group in both the PCA and structure analysis, is most likely a migrant from the Maui/O'ahu group that was visiting or migrated to breed on the Big Island. Thus, complex population structure must be taken into account when interpreting population genetics in highly mobile species (Bowen et al., 2005).
Relatively high genetic diversity despite population declines has been observed in other long-lived endangered seabird species (e.g., the Hawaiian Petrel Pterodroma sandwichensis, Welch et al., 2012; the Balearic Shearwater Puffinus mauretanicus, Genovart et al., 2007; the Magenta Petrel Pterodroma magenta, Lawrence et al., 2008), and may be explained by evolutionary history. An ancient large population of H. castro may lead to retained genetic diversity (Goossens et al., 2005). Despite a likely population decline since the introduction of nonnative mammalian predators to the Hawaiian Islands within the last 1,100 years (Pyle & Pyle, 2017; this manuscript), relatively high genetic diversity in H. castro is not completely unexpected as only a few hundred individuals may be needed to maintain a majority of genetic diversity (Gaither et al., 2010;Tison, 2014).

CONCLUSIONS
This study found little population structure between Kaua'i and the Big Island, and no inbreeding within the Hawaiian populations of H. castro, indicating that at least some individuals are dispersing among the breeding colonies on these islands to maintain gene flow. However, the fact that bycatch birds from Maui and O'ahu do not assign to the same breeding colony as Kaua'i and the Big Island also supports the existence of a second discrete population in the Hawaiian Islands or possibly outside of Hawai'i (e.g., Japan).
Although a lack of detection at the suspected Maui Nui breeding colony precludes direct testing, this island may host a breeding colony distinct from the others (United States Fish and Wildlife Service, 2016), or there may be unknown temporal separation of nesting populations in the Hawaiian Islands that have yet been tested (Raine et al., 2017) similar to that reported for H. castro in Cape Verde (Monteiro & Furness, 1998;Deane, 2013). Continued efforts to find active colonies in the Hawaiian Islands are essential to assess population connectivity and for species recovery (Young et al., 2019).
Populations of H. castro currently do not appear to be in any danger of a genetically induced extinction vortex (Gilpin & Soulé, 1986). However, they remain vulnerable to other threats (Jones et al., 2008;Croxall et al., 2012;Spatz et al., 2014). Reduced fledging success and adult mortality due to invasive predators continue to impact population growth (Galase, 2019). Predator control, translocation, and related management efforts to increase chick survival, attract conspecifics, help expand colony range, minimize adult mortality, and increase nesting success will be crucial in achieving recovery in this species (Raine et al., 2017;Antaky, Galase & Price, 2019).
In summary, despite a historical population decline, continued small population size, and separation of hundreds of miles among islands, this study finds no evidence that populations of H. castro in the Hawaiian Islands are inbred. H. castro colonies in Hawai'i appear to have escaped any severe genetic bottleneck, and the populations do not seem at risk for an extinction vortex associated with loss of genetic diversity (Gilpin & Soulé, 1986). Nevertheless, the small population size of H. castro warrants continued conservation programs to achieve recovery, as seabirds play an important role in food webs in both marine and terrestrial ecosystems in the Pacific (Hobson, Piatt & Pitocchelli, 1994;Fukami et al., 2006) and hold cultural significance to Hawaiian communities (Kamakau, 1987;Rose, Conant & Kjellgren, 1993;United States Fish and Wildlife Service, 2005;National Park Service, 2006).