Genetic diversity, structure, and effective population size of an endangered, endemic hoary bat, ʻōpeʻapeʻa, across the Hawaiian Islands

Tissue sampling

Previous genetic studies of ʻōpeʻapeʻa consisted of small sample sizes with limited geographic coverage on each island (Baird et al., 2015; Russell et al., 2015; Baird et al., 2017; Pinzari et al., 2020). In order to provide robust analyses capable of addressing population genetics questions, population genetic data from a wider variety of resources was necessary. We assembled a biological sample collection using wing membrane and/or muscle sourced from 339 individuals across four Hawaiian Islands (Fig. 1, Table S1). The mtDNA analyses included 321 individuals (180 males; 141 females) and microsatellite analyses included 298 individuals (176 males; 119 females; three unknown). A total of 280 bats were examined for both mtDNA and microsatellite DNA. Tissues were obtained between 2005 and 2020 from live bats (n = 197) captured in mist-nets according to Kunz & Parsons (2009), and from carcasses (n = 142, 2007–2020) submitted from wind energy facilities, the Hawaiʻi State Division of Forestry and Wildlife, the U.S. Geological Survey’s National Wildlife Health Center Honolulu Field Station, the Hawaiʻi Wildlife Center (2007–2020), and from four museum skins (accessions: BPBM177155, BPBM178067, BPBM184307, BPBM184308) during 1988–1999, vouchered in the Bernice Pauahi Bishop Museum (Honolulu, Hawaiʻi). Carcasses were refrigerated or frozen upon discovery. Muscle tissues from necropsies were stored at −20 °C. A sterile 3-mm biopsy punch was used to obtain tissue from each wing of carcasses and live bats. Wing tissues were stored in NaCl-saturated 20% DMSO or on silica gel desiccant beads in the field, then frozen at −20 °C until extraction. USGS procedures are consistent with the guidelines for capture, handling, and care of mammals according to the Institutional Animal Care and Use Committee (IACUC) through the University of Hawaiʻi (#04-039-17) and the American Society of Mammalogists (Sikes, Care & Mammalogists, 2016). Biological samples from bats were collected under Federal and State of Hawaiʻi endangered species collection permits: USFWS TE003483; DLNR-DOFAW WL05-03 through WL19-52.

Molecular methods

DNA was extracted from wing and muscle tissue using a DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol for purification of total DNA from animal tissues, or a QIAmp DNA Micro Kit (Qiagen, Valencia, CA, USA) with protocol modifications for wing punches (Corthals et al., 2015). Isolated DNA was amplified for 1 mtDNA marker, 20 nuclear microsatellite loci, and two genes on the X and Y sex chromosomes.

To assess population structure, female-biased gene flow, and maternal effective population size, a 657-bp region of the mtDNA gene cytochrome c oxidase I (COI) was amplified using the forward primer HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) and the reverse primer LCOI490 (5′-GGTCAACAAATCATAAAGATATTG-3′) (Folmer et al., 1994). Polymerase chain reaction (PCR) was conducted using Illustra Hot Start mix PCR beads (GE Healthcare, Chicago, IL, USA) in 25-µL reaction volumes, containing 20.5 µL of sterile water, 0.5 µL of each primer (10 µM), and 2.5 µL of genomic DNA template. PCR cycling conditions followed Hebert et al. (2003). PCR products were cleaned using Exo-Sap (Affymetrix, Santa Clara, CA, USA). Sanger sequencing of PCR products using forward and reverse primers was performed via ABI Prism 3500 Genetic Analyzer (Applied Biosystems, Carlsbad, CA, USA) at the University of Hawaiʻi at Hilo Core Genomics Facility (UHH CGF) or ABI 3730xl DNA Analyzer (Applied Biosystems, Carlsbad, CA, USA) by Sequetech DNA sequencing services (Mountain View, CA, USA). Sequence chromatograms were manually trimmed, edited, and contigs formed in Geneious Prime v.2020.2.4. All sequences were checked against GenBank using a BLAST search to confirm a match to ʻōpeʻapeʻa. Sequences produced in this study (n = 262) were combined with ʻōpeʻapeʻa COI sequences from (Russell et al. (2015); n = 59), available on GenBank (accession numbers KR350020 through KR350078). The combined dataset was aligned and trimmed by eye using MEGA-X v10.2.0 software (Kumar et al., 2018). The total COI dataset for this study consisted of 321 individuals from Hawaiʻi (H, n = 166), Maui (M, n = 92), Oʻahu (O, n = 47), and Kauaʻi (K, n = 16), spanning collection years 2005–2019 (Table S1).

A suite of 20 nuclear microsatellite loci developed for the conspecific northern hoary bat (L. cinereus) were amplified in ʻōpeʻapeʻa samples (Table S2) following multiplex methods for PCR and genotyping described in Korstian et al. (2013); Korstian, Hale & Williams (2015); Korstian, Hale & Williams (2014) and Pylant et al. (2016). Loci were amplified in five multiplex groups (E, F, MP1, MP2, and MP5) using a Qiagen Multiplex Reaction Kit (Qiagen, Valencia, CA, USA). Amplified microsatellite multiplexes were diluted in sample loading solution with size standard (CEQ 8000 DNA standard kit 500 Beckman Coulter or ABI DS33/LIZ500 Ladder ABI Systems) and electrophoresed on a Beckman Coulter CEQ 8000 (at UHH CGF) or an Applied Biosystems 3,730 × l Genetic Analyzer (GENEWIZ, South Plainfield, NJ, USA). Genotypes were sized and scored using Beckmann CEQ 8000 fragment analyzer package or Geneious Prime v.2020.2.4. Because fragment size assessment was conducted on two different platforms, the allele calls were standardized and checked for allelic dropout by running 8% (n = 24) of samples on both platforms and all microsatellite data were analyzed using Tandem (Matschiner & Salzburger, 2009), a program designed for aligning fragment sizes and merging cross platform allele bins. Genotyping errors from allelic dropout, stuttering error, and null alleles in microsatellite loci were screened using Micro-Checker 2.2.3 (Van Oosterhout et al., 2004).

Genotyping was attempted at 20 microsatellite loci for a total of 298 samples collected by island as follows: H, n = 131; M, n = 102; O, n = 49; and K, n = 16 (Tables S1, S3). We discarded the CotoG12 locus early on, as it showed presence of null alleles. Exact tests for deviations between observed (H_o) and expected (H_e) heterozygosity, based on Hardy-Weinberg Equilibrium (HWE), were performed in Arlequin v3.5.2.2 (Excoffier & Lischer, 2010), whereas likelihood ratio tests for genotypic linkage disequilibrium (LD) were conducted in GenePop v4.7 (Raymond & Rousset, 1995; Rousset, 2008). For both HWE and LD, Bonferonni corrections for multiple comparisons were applied to p-values with the p.adjust function in the R “stats” package v4.0.5 (R Core Development Team, 2020). The program FreeNA, with the ENA method (Chapuis & Estoup, 2007) corrected for positive bias in estimates of microsatellite pairwise F_ST values due to the presence of null alleles.

Mitochondrial sequences, microsatellite genotypes, and detailed tissue collection information produced through this study are available from https://doi.org/10.5066/P9COQ3ZK (Pinzari et al., 2022). The mtDNA COI sequences were also placed in the NCBI GenBank repository under accession numbers OL894241 – OL894502.

Genetic diversity

Patterns of genetic diversity across the four islands were inferred from ʻōpeʻapeʻa COI sequences. Per-island, average pairwise genetic distances (Dxy) and percent sequence divergence were calculated using MEGA-X v10.2.0 (Kumar et al., 2018), whereas diversity statistics number of polymorphic sites (S), haplotype diversity (Hd), per site nucleotide diversity (π), and Watterson’s theta (Θ) (per site and per sequence) were calculated using DnaSP v6.12.03 (Rozas et al., 2017). Those diversity statistics were calculated for individuals island and partitioned across collection years with a minimum of nine samples per time category. To account for uneven sample sizes across islands, the COI sequences were also analyzed by the minimum number of sequences available for each island (with sequence sub-sets chosen randomly). Accordingly, those datasets consisted of sequences numbering 16 (K, O, M, H), 47 (O, M, H), and 92 (M, H). Patterns of haplotype diversity were also assessed by examining relationships among collapsed unique COI haplotypes using a maximum parsimony TCS network (Clement, Posada & Crandall, 2000) and visualized using PopArt Population Analysis with Reticulate Trees (Leigh & Bryant, 2015).

Genetic diversity was assessed for nuclear microsatellites, with the number of polymorphic loci, number of alleles, expected and observed heterozygosity, gene diversity, and inbreeding coefficients calculated in Arlequin v3.5.2.2 (Excoffier & Lischer, 2010). To estimate allelic richness while accounting for uneven sample sizes, a rarefaction analysis was conducted using HP-Rare v1.1 (Kalinowski, 2005), with a minimum gene size of n = 32 gene copies (based on the smallest population sample size - Kauaʻi).

Population structure and isolation by distance

The hypotheses that ʻōpeʻapeʻa populations are genetically structured by island, and that regional population samples may show genetic subdivision within islands, was tested using model and non-model based approaches on both mtDNA haplotypes and microsatellite allele frequencies. Population genetic differentiation at island-scale was assessed using pairwise F_ST and analysis of molecular variance (AMOVA) of mtDNA and microsatellite data, calculations were conducted in Arlequin v3.5.2.2 (Excoffier & Lischer, 2010). The microsatellite analysis was conducted using two datasets, observed alleles (uncorrected), and alleles corrected for null alleles using FreeNA with the ENA method (Chapuis & Estoup, 2007). Genetic structure was also assessed using two multivariate clustering methods, principal components analysis (PCA) and discriminant analysis of principal components (DAPC), implemented in the packages “ade4” v1.8 (Dray & Dufour, 2007) and “adegenet” v2.1.3 (Jombart, 2008) using R v4.1.3 (R Core Development Team, 2022). The PCA identifies components of variation that contribute the greatest percentage differences observed within the genotype data, whereas DAPC constructs linear combinations of alleles having the largest between-group variance and the smallest within-group variance (Jombart, 2008).

Genetic structure and potential for admixture was also modeled by calculating ancestry coefficients and individual sample cluster membership (K) with the program Structure v2.3.4 (Pritchard, Stephens & Donnelly, 2000). Because our total dataset of multilocus microsatellite genotypes contained uneven numbers of individuals per island, and unbalanced sampling underestimates latent population structure when using the program Structure (Puechmaille, 2016; Wang, 2017), we ran analyses for several different scenarios to alleviate potential biases: (1) by island using the total dataset while allowing for uneven sample sizes (298 individuals; H = 131, M = 101, O = 49, K = 16), and (2) with mostly even sample sizes (106 individuals; H = 30, M = 30, O = 30, K = 16). Structure parameters were set to assume admixture, correlated allele frequencies within populations (Falush, Stephens & Pritchard, 2003), and the prior for admixture, α, was set to vary between clusters and decreased from 1.0 to 0.1 as suggested in Wang (2017). We used a range of 1-10 for values of K and 20 iterations per K, with a burn in length of 1 ×10⁵ and 2 ×10⁵ Markov chain Monte Carlo simulations (MCMC). The interpretation of Structure output was aided by Structure Selector (Li & Liu, 2018), which computes the following measures: L(K) the maximal estimate of posterior probability for a given K; ΔK, an ad hoc estimate of the most probable number of K-clusters based on the rate of change in the log probability of the data between consecutive K-values (Evanno, Regnaut & Goudet, 2005); and supervised measures (MedMedK, MedMeanK, MaxMedK, and MaxMeanK) under three varying thresholds for ancestry coefficients (0.5, 0.6, and 0.7) as recommended in Puechmaille (2016). The ancestry coefficients were visualized across clusters using Clumpak (Kopelman et al., 2015).

Isolation by distance as a potential explanation for restricted gene flow or population structure between and within islands, was tested using a Monte-Carlo style Mantel test (with 999 replicates) for correlations between microsatellite-based matrices of Edward’s genetic distance and Euclidean geographic distances using package “adegenet” v2.1.3 in R (Jombart, 2008); and (2) regression of genetic distance ((F_ST/(1- F_ST)) and log transformed geographic distance in GenePop (Rousset, 2008). Isolation by distance tests were calculated at group level, represented by bats collected at 59 different collection locations across the Hawaiian Islands and within islands: H = 27, M = 21, O = 5, and K = 6.

Historical and contemporary gene flow

The historical migration patterns of ʻōpeʻapeʻa among the islands of Hawaiʻi, Maui, and Oʻahu could include unidirectional or bidirectional movement across islands, or archipelago-wide panmixia. Using Migrate-n v.3.7.2 (Beerli, 2006; Beerli & Palczewski, 2010), we considered six different modeled migration scenarios: (1) uneven versus even migration rates between all pairs of islands; (2) archipelago-wide panmixia; (3) Maui as the source population of Oʻahu and Hawaiʻi populations (with no self-recruitment); (4) directional migration from Hawaiʻi to Maui and Oʻahu; (5) directional migration from Maui to Oʻahu and Hawaiʻi; and (6) directional migration from Oʻahu to Maui and Hawaiʻi. Migrate-n recommends limiting population modeling to two or three populations, thus we excluded Kauaʻi from this analysis. Microsatellite data were recoded using the amplicon length divided by repeat number (standardized in Tandem), and a Brownian motion model was used to estimate the stepwise mutation of microsatellites. Uniform prior distributions of theta were bound at 0 and 150, window size 15, and migration parameter (M) bound at 0 and 800, with window size set to 80 and the mutation rate option set relative to the data. Model runs consisted of 25 million generations, sampling every 500 generations, and using a burn-in of 6.25 million. Four variably heated chains were used following the methods of (Beerli & Palczewski, 2010). Run convergence was determined from posterior probability distributions. The Bayes factors and model probabilities were calculated using log marginal likelihoods from Bezier approximations with a custom python script provided in Beerli et al. (2019).

To provide measures of contemporary migration rates (0–2 generations previous) that may have occurred recently between Hawaiʻi, Maui, Oʻahu, and Kauaʻi, microsatellite data were analyzed using BayesAss v3.0.3 (Wilson & Rannala, 2003). Each island was treated as an individual population with posterior probability of migrant history estimated from individual ancestry measures. The program was run for 50 million generations with a burn-in of 10 million, sampling every 1,000 generations, and parameter settings as follows: step sizes of migration rate m = 0.1, allele frequencies α = 0.4, and inbreeding coefficients f = 0.7. Several pilot runs were conducted using varying parameters to ensure the acceptance rate of final parameters remained between 40 and 60% throughout the duration of the run. Run convergence was visualized in Tracer v 1.6 (Rambaut et al., 2018) to ensure adequate mixing through generations. The posterior probability of migrant history was inspected for each individual.

To complement the BayesAss analysis, recent migration was assessed with the gensback option (set at 2) in Structure v2.3.4 (Pritchard, Stephens & Donnelly, 2000). Sampling locations were used as priors to identify likely non-residents among a population sample or individuals that might have descended from a recent migrant. For this analysis, the number of clusters was set to K = 4, based on model (K) selection. The migration rate priors were varied (0.01, 0.05, 0.1) to assess sensitivity (Pritchard & Wen, 2004) and we used a burn-in of 200,000 and run length of 500,000 MCMC chains.

Sex-biased dispersal testing

For carcasses that could not be reliably sexed from morphological features, genetic methods were employed to identify sex for tests of sex-biased gene flow. Of 113 samples that lacked morphological sex data, sex genotyping was performed by amplifying regions of introns on the Zfx and Zfy sex chromosomes using PCR and visualizing fragments by electrophoresis (Korstian et al., 2013; Pinzari & Bonaccorso, 2018).

Models to detect sex-biased dispersal have several assumptions, including nonoverlapping generations, dispersal of juveniles prior to reproduction, and that sampling occurred post-dispersal. Like many species, bats have overlapping generations, yet sex-biased dispersal tests have been used to explore dispersal in bats (Petit, Balloux & Goudet, 2001; Kerth, Mayer & Petit, 2002; Salgueiro et al., 2008). Because detecting sex-biased dispersal depends on the proportion of population sampled and thus greater sample size is preferred, we performed tests using microsatellite loci combined from all islands (n = 295, 4 groups) and for regional population samples within each island; H = 130, 4 groups; M = 99, 4 groups; O = 49, 2 groups; and K16, 2 groups. To compare dispersal probabilities between sexes of adults in our data, we evaluated the mean and variance assignment indices (mAI_c, vAI_c), F_ST and F_IS, relatedness, observed heterozygosity (H_o) and gene diversity (H_s) measures as described in Goudet, Perrin & Waser (2002). We conducted sex bias testing with the p one-sided hypothesis testing option using FSTAT v2.9.4 (Goudet, 2003).

Historical and contemporary effective population sizes

The long-term historical female (N_Ef) and long-term (N_e) effective population sizes were estimated using complementary mitochondrial and microsatellite DNA techniques. The mtDNA COI sequences represent the long-term historical female effective population size from tens to thousands of generations before present. These estimates were calculated with the equation θ=2N_eu, with the Watterson’s estimator (θ) value determined using Arlequin v3.5.2.2 (Excoffier & Lischer, 2010). The per sequence per generation mutation rate, u, (Tajima, 1993; Schenekar & Weiss, 2011) for COI has not been published for the genus Lasiurus, therefore, following Korstian, Hale & Williams (2015), two different substitution rates based on cytochrome b in the family Vespertilionidae were applied as proxies (Nabholz, Glémin & Galtier, 2008). Adjusted to our COI sequence length (657 bp), u was set to 8.061 ×10⁻⁶ and 1.089 ×10⁻⁴. Long-term N_e (i.e., global population) size was estimated from the average expected heterozygosity (H_e) of the microsatellite loci using two mutation models, the infinite allele model (IAM) and the stepwise mutation model (SMM). Here, N_e was calculated using the respective equations: SMM N_e = { (1/1-H)²-1)/8µ, and IAM N_e = H/(4µ(1-H) } (Nei, 1987), with mutation rates for microsatellite loci set to 10⁻³ or 10⁻⁵, consistent with Korstian, Hale & Williams (2015). The N_e values were calculated by averaging the IAM and SMM model results for each mutation rate.

The null hypothesis that per-island long-term N_Ef and N_e estimates do not differ from contemporary (recent generations) effective population size (N_eC) was tested using the microsatellite linkage disequilibrium method in NeEstimator V2 (Do et al., 2014). The linkage disequilibrium single-sample method relates a decrease of N_e with an increase in genetic drift. Estimates were generated for each island with all samples over all collection years, and individual island datasets were also subdivided into temporal groupings in cases where collection year spanned a minimum two-year period sample size. Random mating was assumed and alleles with frequencies below 0.05 were removed. Estimates of N_e are reported with ± 95% confidence intervals determined with both jackknife and parametric re-sampling methods.

Demographic patterns—bottlenecks and expansions

To assess whether island population demographic signatures indicate past expansions and contractions, microsatellite data were examined using a modified Garza-Williamson M-ratio test and COI sequences were subjected to tests of neutrality. The M-ratio method compares the per-locus ratio of the number of size classes with at least one observed allele to all possible allele categories calculated from the full set of allelic size ranges, with the assumption that a bottleneck event will result in loss of alleles from allele size categories. M-ratio values below a critical value of 0.68 generally indicate a recent population bottleneck (Garza & Williamson, 2001). For each island, the M-ratios were empirically calculated in Arlequin v3.5.2.2 (Excoffier & Lischer, 2010) and compared to critical values (M_c) simulated for each island using the software Critical M (Garza & Williamson , 2001). This analysis was run by setting the probability of changes greater than one step (ρg) to 0.2, the size of one-step changes (Δg) to 3.5, and varying θ from 0.01 to 10, which encompassed all estimates of θ produced by mitochondrial sequences (Busch, Waser & Dewoody, 2007.

The neutrality tests on the COI dataset were conducted using DnaSP v6.12.03 (Rozas et al., 2017). Two summary statistics, Tajima’s D (Tajima, 1989) and Fu’s F (Fu, 1997), were applied to test the null hypothesis that the number of rare and common mutations are equal, as is expected for a stable population, while the alternative is an increase in low frequency mutations consequent to a population undergoing expansion or rapid growth.

Genetic diversity

Inter-island differences in molecular population genetic diversity were prominent in the analysis of the mtDNA COI gene, (Tables 1 and 2). The haplotype and nucleotide diversity values were variable between islands, and between years (Table 1, Table S4). Among sampled islands, Hawaiʻi showed the greatest number of unique haplotypes and polymorphic sites, while Kauaʻi showed the lowest number of unique haplotypes and fewest polymorphic sites (Table 1). Haplotype diversity ranged from 0.669 (M) to 0.125 (K) while nucleotide diversity ranged from 0.012 (M) to 0.0005 (K). The low mitochondrial diversity observed on Kauaʻi does not appear to be an artifact of sample size because results for Hawaiʻi and Maui were mostly robust across the number of rarefied sequence sets (Fig. S1A); however, nucleotide diversity observed for Oʻahu dropped considerably when sample sizes were increased (Table 1, Fig. S1B).

The average pairwise genetic distance values ranged from 4.271 to 0.125, with differences generally greater between than within island populations, and Maui and Hawaiʻi showing the greatest spread (Table 2). The percent sequence divergence ranged from 2.63% to 0.06%, with Kauaʻi- Maui and Hawaiʻi-Maui, showing greatest differences (>2.4%).

Consistent with other studies the mtDNA haplotype network indicated the presence of two clades (Fig. 2). Among the 22 unique haplotypes detected in 316 individuals, the network indicated three haplotypes were shared across islands while most were restricted to individual islands. Sequences of 5 individuals from Oʻahu (O8, O9, O14, O17, and O27) displayed heteroplasmy. None of these sequences showed evidence of contamination and all sequences were translatable into amino acids with no stop codons that would indicate they were nuclear copies of mitochondrial DNA (Numts). Most mutations were transitions and the heterozygous sites occurred across several individuals from the same collection site. Removal of these Oʻahu individuals resulted in a haplotype network with a much greater number of base pair changes, emphasizing distinctiveness of two mitochondrial clades (Fig. 2). However, when Oʻahu individuals were included, network results display a closer relationship between the two mtDNA clades (Fig. S3).

The majority of samples (n = 278, 93.2%) were genotyped at all 19 microsatellite loci, with a minority of samples (n = 20, 6.8%) genotyped at only 9–18 loci (Table S3). The number of alleles per locus ranged from 2–24 (Table S3). No evidence for consistent or significant linkage disequilibrium was found between any two loci across the four islands following Bonferroni corrections, thus all loci were considered independent, however the locus LAS8539AC was discarded, because it returned strong evidence of deviation from HWE across islands and presence of null alleles. Across the remaining 18 microsatellite loci, the per-island number of polymorphic loci ranged from 15 (O and K) to 18 (H), the per-locus and per-island mean numbers of alleles ranged from 4.93 (K) to 9.22 (H), and the observed heterozygosity between bats on different islands ranged from 0.56 (O and H) to 0.59 (M). The inbreeding coefficients (F_IS) ranged from 0.013 (H) to −0.035 (O) and were not significant (p-value threshold 0.05), consistent with absence of the Wahlund effect within individual island datasets (Table 3). The rarified private allelic richness varied by island; 1.29 (H), 0.83 (M), and 0.25 (O and K). Partitioning data by years showed the average microsatellite diversity measure were mostly stable between time periods (Table S5).

Genetic structure

An AMOVA test for mitochondrial haplotype differences showed significant differentiation among and within islands, with variance among islands at 67% compared to within islands at 33%, F (3, 320) = 0.669, p < 0.001. In contrast, the AMOVA conducted for microsatellite loci showed a very small percent of variance due to differences among islands (7.60%) or within islands (1.76%), with 90.6% of variance due to the differences within individual bats (F = 0.076, p < 0.001). Consistent with population structure, global F-statistics across all microsatellite loci were significant (p <0.001: F_IS = 0.019 [95% confidence interval (CI) [−0.005–0.04], F_ST = 0.076 [95% CI [0.05–0.09]], and F_IT = 0.093 95% CI [0.06–0.12]).

Pairwise F_ST of nuclear microsatellites were lower than those calculated using mtDNA COI, a pattern that is typical of the two marker types. However, both datasets revealed significant genetic differentiation between all island pairs (p < 0.05), with the exception of Hawaiʻi and Kauaʻi (Table 4) at the COI marker. The highest significant pairwise differences for COI were between Maui and Hawaiʻi and between Maui and Kauaʻi. Although F_ST values slightly decreased when corrected for null alleles using FreeNA with the ENA method, they still showed the same pattern as the uncorrected microsatellite data (Table 4, Table S6). The lowest F_ST values for the microsatellite data occurred in pairwise comparisons between Maui and Hawaiʻi, and Oʻahu and Kauaʻi.

Population structure and isolation by distance

Genetic subdivision by island was evident from PCA and DAPC scatterplots, which showed mostly non-overlapping scatter for each island, except between Oʻahu and Kauaʻi (Fig. 3). For PCA, axis 1 (7.85% total variation) separated Oʻahu and Kauaʻi from Hawaiʻi and Maui, axis 2 (5.49% total variation) separated Maui from the other islands (Fig. 3A), and axis 1 with axis 3 (3.56% total variation) further separated Kauaʻi and Oʻahu from Hawaiʻi and Maui (Fig. 3B). For DAPC, a plot of components 1 and 2 indicated three distinct clusters with overlap between Oʻahu and Kauaʻi populations (Fig. 3C), while the plot of components 1 and 3 indicate four genetic clusters defined by island with further separation of Kauaʻi (Fig. 3D).

Results from the program Structure were concordant with PCA and DAPC scatterplots, revealing strong genetic differentiation among nearly all island populations. The number of genetic groups varied from three to four, contingent on sample datasets and based on post hoc estimators (Fig. 4). Structure results from the whole dataset supported differentiation between Hawaiʻi and Maui clusters, and a third cluster composed of individuals from both O′ahu and Kauaʻi, with Evanno and Puechmaille post hoc estimators supporting K = 3 (Fig. 4A). When the numbers of individuals across islands was roughly even, individuals clustered into island groups at K = 5 and above (Fig. 4B); however, post hoc metrics resulted in K = 3 for both Evanno and Puechmaille estimators. Post hoc metrics for even sampling resulted in K = 3 for the Evanno method but K = 4 for the Puechmaille estimator. When testing for potential population substructure within larger islands of Hawaiʻi and Maui, clustering patterns within Maui resulted in K = 2 for the Evanno method and K = 3 for the Puechmaille estimator. We did not find overwhelming evidence for distinct genetic clusters using these methods within the Island of Hawaiʻi. Members of the two mitochondrial clades were dispersed across the clusters identified by the Structure analysis (Fig. S4).

Genetic structure across all island bat populations was further corroborated by the Mantel test correlation between matrices of Edward’s genetic distance and Euclidean geographic distances across collection sites on four islands (r = 0.224, p = 0.001; Fig. 5A). Results of isolation by distance (IBD) Mantel tests conducted with adegenet (Jombart, 2008) within islands supported a relationship between genetic and geographic distance across collection sites within Hawaiʻi (r = 0.322, p = 0.002; Fig. 5B), collection sites within Maui (r = 0.265, p = 0.056; Fig. 5C), but not across collection sites within Kauaʻi (r = 0.260, p = 0.075) or O’ahu (r = −0.156, p = 0.597). Regression of genetic distance [(F_ST/(1- F_ST)] and log transformed geographic distance supported the alternative hypothesis that distance had an effect on genetic differentiation among collection sites across four islands (b = 0.04, 95% CI [0.02–0.04]), and within Maui (b = 0.05, 95% CI [0.03–0.10]). However, the null hypothesis could not be rejected for collection sites within Hawaiʻi, Oʻahu, or Kauaʻi.

Dispersal and gene flow

Evidence for sex biased dispersal was inconsistent across test metrics and groups (Table 5). Across all islands, assignment indices (mAIc and vAIc), and measures of inbreeding coefficients (F_IS) measures had a pattern of weak male dispersal; however, this was not significant after p one-sided hypothesis testing. Only samples from Maui showed significance for male biased dispersal in mAIc, vAIc, and H_S. Within Hawaiʻi, mAIc, F_ST, and F_IS measures had a pattern of female-biased dispersal, but vAIc supported male dispersal. Within Oʻahu, mAIc supported female dispersal, while vAIc and F_ST indicated male dispersal; however, F_IS supported no dispersal of either sex. Within Kauaʻi, mAIc, F_ST, and F_IS supported female dispersal. The p one-sided hypothesis testing was not significant for any measures within Hawaiʻi, Oʻahu, or Kauaʻi. Within Maui, male biased dispersal was supported with mAIc, vAIc, and F_IS measures, and for mAIc, vAIc, and H_S was significant after hypothesis testing.

Results from BayesAss produced contemporary migration proportions, with confidence intervals including zero, thus the proportion of migrants was not significantly different from zero (Fig. 6A). Rare dispersal events may have occurred more than two generations ago, as BayesAss identified a handful of individuals on Hawaiʻi, Maui, and Oʻahu with migrant ancestry posterior probabilities, but in most cases less than 0.30 probability of being a second-generation migrant. Recent migrant ancestry results from Structure using several migration prior settings, also resulted in little to no evidence for contemporary migration (Fig. S4). Only two individuals from Oʻahu showed probability of assignment to a different island; however, the assignment probabilities were mixed across generations and populations from Hawaiʻi and Maui and decreased as migration priors increased.

Comparison of Bayes factors and model probabilities across six dispersal models from Migrate-n resulted in the probability of 1.00 for a model placing Maui as the “source” population of Oʻahu and Hawaiʻi migrants . In this scenario, Maui’s historical theta {Θ (4N_eµ)} was ∼3 times larger than O′ahu and ∼1.4 times larger than Hawaiʻi (Fig. 6B). The mutation-scaled immigration rates (M = m/µ, m = proportion of immigrants) from Maui to Oʻahu was less than the rate from Maui to Hawaiʻi. With the mutation rate for mammalian microsatellites (µ= 1.2 ×10⁻³) from Weber & Wong (1993) applied, the rough estimate of historical migration rate from Maui to Oʻahu was 0.45 and to Hawaiʻi was 0.57 individuals per generation.

Estimates of effective population size

Long-term historical female effective population size (N_Ef) estimates from COI mitchondrial sequence data were similar across Hawaiʻi, Maui, and Oʻahu (min 20,796 and max 286,658, ranges contingent on mutation rate) with the Kauaʻi sample showing a smaller range (1,384–18,194; Table 6). Using methods based on heterozygosity, the long-term N_e estimates from microsatellite data were similar across the four islands with wide ranges and degree of overlap. However, the SMM estimates for Oʻahu and Kauaʻi were half that of Hawaiʻi and Maui (Table 6). Long term Ne was generally lower than N_Ef based on COI (except for the low N_e based on COI for Kauaʻi). The long-term N_e trended higher for the SMM model compared to the IAM, with N_e bracketing a wide range of 248 to 62,162. Estimates from averaged results between the model types for each mutation rate had a pattern of larger N_e for Hawaiʻi and Maui (∼490 to 49,000), and smaller N_e for Oʻahu and Kauaʻi (∼310 to 31,000). The average estimates using the lower mutation rate produced values that fell within the range of historical N_Ef estimates from mitochondrial sequences.

Linkage disequilibrium estimates of contemporary effective population sizes (N_eC) of ʻōpeʻapeʻa across all islands and collection years were far smaller than historical long-term N_Ef and N_e, but confidence intervals included infinity for Hawaiʻi and Kauaʻi, which precludes making strong inferences and indicates potentially larger N_eC populations (Table 7, Table S7). Estimates of N_eC were largest for Hawaiʻi, similar between Maui and Kauaʻi, and smallest for Oʻahu (Table 7). The Maui and Oʻahu N_eC estimates resulted in bounded confidence intervals. Hawai′i Island N_eC for annual groups in 2009–2012 and 2018–2019 were 391 and 350 respectively, while the 2018–2019 parametric confidence interval was the only metric that produced an upper bound (4,262 in Table S7). Maui’s N_eC for annual groups 2012–2014 and 2016–2019 were 65 and 95, respectively, and each had confidence intervals with upper bounds <500. Oʻahu’s N_eC for annual groups 2013–2015 and 2017–2020 were 21 and 11, respectively, with confidence interval upper bounds <50. In some cases, the sample size exceeded the N_eC estimate.

Demographic patterns–bottlenecks and expansions

The modified Garza-Williamson M-ratios for each island population resulted in M-ratio mean values lower than the simulated range of critical M (M_c), indicating either past population size reductions or possible bottleneck events across Maui (0.43 ± 0.20, M_c range [0.69–0.80]); Oʻahu (0.32 ± 0.20), M_c range [0.65–0.80]); and Kauaʻi (0.25 ± 0.15), M_c range [0.58–0.80]). For Hawaiʻi the M-ratio standard deviation encompassed critical M range (0.56 ± 0.29, M_c range [0.70–0.81]) The Tajima’s D neutrality test indicated a past expansion for the Hawaiʻi population, showing significant negative values across most rarefied groups (Table 8), however Maui and Oʻahu populations showed no evidence of expansions. Tajima’s D tests conducted on subsets of samples by collection year for each island, continued to support past expansion signals for Hawaiʻi for both collection periods(Table S9). The Fu’s Fs testing did not significantly support past population expansions (Tables S8, S9).

Genetic diversity

Insular patterns from mitochondrial and nuclear loci generally mirrored heterozygosity values reported for single nucleotide polymorphisms in (Pinzari et al., 2020), with higher diversity values reported for Maui, and lower diversity values reported for Hawaiʻi and Oʻahu. Our mitchondrial haplotype network results, which included only individuals from the Hawaiian Islands, were similar to the mitochondrial clade networks described in phylogenetic studies of Russell et al. (2015) and Baird et al. (2017) when individuals with point mutations suggestive of heteroplasmy were removed. Our increased sampling and geographic coverage presented 11 additional haplotypes. One clade occurred across four islands, although it was most abundant on Hawaiʻi. The second clade, most abundant on Maui, also was present in low numbers on Oʻahu and Hawaiʻi. Peak mitochondrial haplotype diversity on Maui reported by Russell et al. (2015) and Baird et al. (2017) was supported by our results. Sampling to include Lanai and Molokaʻi Islands may uncover additional diversity increasing our understanding of the unique haplotype dynamics occurring on Maui.

The presence of heteroplasmy in five individuals from Oʻahu is unique among Hawaiian bats, although heteroplasmy occurs at high rates in vespertilionid bats elsewhere (Petri, Von Haeseler & Pääbo, 1996; Wilkinson et al., 1997). Heteroplasmy may result from paternal mtDNA leakage, or as a function of de novo somatic mutations, or from maternal transmission via heteroplasmic eggs. Pinzari et al. (2020) reported the highest number of non-synonymous mutations within Oʻahu bats, thus observations in our study could result from mutational processes or oxidative mutations linked to infection or stressful reproductive effort (note that four of five individuals in our study were males) (Jebb et al., 2018).

In comparison to northern hoary bat populations, ʻōpeʻapeʻa exhibited low heterozygosity with similar levels among islands. Also, fewer alleles per loci occur in ʻōpeʻapeʻa compared to the northern hoary bat (Korstian, Hale & Williams, 2014; Korstian, Hale & Williams, 2015; Keller et al., 2014). A caveat is that downward bias of genetic diversity estimates from microsatellite primers transferred across species may occur when a small number of loci are selected based on high polymorphism in the focal species because congeners may have fixed or relatively low polymorphisms (Allendorf, Luikart & Aitken, 2012). Alternatively, it is possible that fewer alleles in the ʻōpeʻapeʻa population is due to founder effects and subsequent genetic drift as the ancestral population expanded. If populations are approaching or have reached migration-drift equilibrium, genetic differentiation within islands is not expected to be greater than differentiation among islands. Drift is strongest in small populations and thus may explain the lower genetic diversity on small islands like Oʻahu and Kauaʻi. Isolated populations also may show reduced allelic richness and gene diversity (Broquet et al., 2010), which aligns with those measures for Oʻahu and Kauaʻi and the geographic distance of those islands from the others.

Population structure

Evidence of population structure using both mtDNA and microsatellite markers is absent in localized studies of L. cinereus (Korstian, Hale & Williams, 2015; Pylant et al., 2016), thus indicating substantial exchange of genetic material via dispersal and mating opportunities (Waples & Gaggiotti, 2006). In contrast, population genetic structure might be expected in island bat species where dispersal is limited by broad oceanic channels and/or resources are annually abundant within an island. Genetic isolation of island bat populations descending from continental populations has been demonstrated for molossids and vespertilionids, having large open water crossings as barriers to gene flow; examples include Myotis punicus in the Mediterranean (Biollaz et al., 2010), Nyctalus azoreum in the Azores (Salgueiro et al., 2004; Salgueiro et al., 2007; Salgueiro et al., 2008), Tadarida brasiliensis in the Bahamas (Speer et al., 2017), and two species of Miniopterus in Madagascar (Weyeneth et al., 2011). We found a weak signature of IBD across the four islands in this study and within Hawaií and Maui, indicating that geographic distance may restrict dispersal in ʻōpeʻapeʻa to some extent. Although heavily sampled, there was not strong support for IBD within Hawaiʻi from both methods used. Our number of sampling sites for Oʻahu and Kauaʻi were likely too few to enable detection of such patterns.

Nuclear microsatellite clustering in Structure analyses adds additional evidence to an overarching pattern of island structure for ʻōpeʻapeʻa reported with SNPs in Pinzari et al. (2020). In this study, three clusters have a clear signal using data from all available individuals. However, analysis conducted using similar numbers of individuals per island revealed that a greater number of genetic clusters could be present with sub-structure within Maui evident and that Kauaʻi was distinct from other islands. Increased sampling density could better test hypotheses of regional population structure within Oʻahu and Kauaʻi. Results from Bayesian clustering programs like Structure can overestimate genetic structure when IBD is also present across a continuously distributed population (Frantz et al., 2009; Perez et al., 2018); however, this may have limited effect on our results because ʻōpeʻapeʻa populations are strongly structured by island and the signal of IBD was not strong.

Estimates of genetic differentiation for nuclear microsatellites maintained significance in population structure between islands but were subtler than mitochondrial results. The large values of population structure estimated from mtDNA may be due to the smaller effective population size of mtDNA resulting in greater effects of genetic drift on population structure for mtDNA (Allendorf, Luikart & Aitken, 2012). The latter may occur when populations have been isolated over short periods of time (Zink & Barrowclough, 2008). In addition microsatellites may underestimate F_ST if high mutation rates cause size homoplasy even if migration rates are low (Balloux et al., 2000).

Migration history and gene flow

Model estimates of historical connectivity indicate that dispersal throughout Maui, Hawaiʻi, and Oʻahu may have always been limited; for example, if estimated movements were less than 1 individual per generation. However, these results present a challenge for comparisons due to the differing model assumptions, especially if the true scenario of island migration is not captured in the historical Migrate-n models. For example, during times of lower sea level when the regional ‘Maui Nui super-island complex’ was contiguous (Price & Elliot-Fisk, 2004), genetic separation would not have been maintained by spatial distance or large oceanic channel barriers. Our Migrate-n models selected a model that placed Maui as the source of Hawaiʻi and Oʻahu migrants. This model result and genetic diversity values, are additional evidence that Maui represents the founding population for Hawaiian ʻōpeʻapeʻa as hypothesized in Russell et al. (2015); Baird et al. (2017), and Pinzari et al. (2020). Additionally, historical patterns of gene flow in ʻōpeʻapeʻa contrast with the biogeographical stepping-stone model or “progression rule” of colonization modeled to occur across the Hawaiian archipelago. In this model, Kauaʻi is the oldest of islands that formed along the volcanic hot spot (Shaw & Gillespie, 2016; Wagner & Funk, 1995). Contemporary migration estimates from microsatellites using BayesAss did not support gene flow between islands over recent generations, nor did we detect recent migrants using alternative methods in Structure. Although the estimates of recent movement between islands we generated were not significant, even a few individuals moving between populations per generation could maintain levels of gene flow. It appears that dispersal events between islands are incredibly rare for ʻōpeʻapeʻa, but could potentially be triggered by extreme weather events, e.g. (such as hurricanes) (Fleming & Murray, 2009). It is also possible that the number of bats and the time intervals sampled were not broad enough to detect long term temporal gene flow patterns.

Diversity differences found between microsatellite loci and mtDNA may indicate sex-biased dispersal in bat species, in addition to inherent differences between mtDNA and nuclear DNA patterns of divergence (Moussy et al., 2013; Petit, Balloux & Goudet, 2001). Male-biased dispersal occurs in forest and crevice dwelling vespertilionid bat species including, Myotis bechsteinii (Kerth, Mayer & Petit, 2002), Myotis septentrionalis (Arnold, 2007), and Vespertilio murinus (Safi, König & Kerth, 2007). Female-biased dispersal has been documented in Nyctalus azoreum in the Azores Salgueiro et al. (2008). If female bats return to the same roost areas each breeding season and dispersal is male-biased, this will result in mtDNA showing greater population differentiation than microsatellites. We found greater population differentiation in mtDNA, potentially indicating a pattern of female philopatry to islands. However, sex-biased hypothesis testing revealed no significant sex-biased dispersal across or within Hawaiʻi, Oʻahu, and Kauaʻi for either sex. On Maui, our data did support male-biased sex dispersal. Sex-biased tests have limited performance when dispersal rates are low, and limited power in cases where the bias in dispersal is greater than 80:20 (Goudet, Perrin & Waser, 2002). When dispersal is less than 10%, the vAIc statistic can be informative, but normally F_ST is considered robust to changes in sampling and magnitude, whereas mAIc has been described as a measure between the two (Goudet, Perrin & Waser, 2002). Our vAIc statistic supported male-biased dispersal patterns within Hawaiʻi and Oʻahu but F_ST and mAIc metrics did not align and conversely indicated potential female-biased dispersal. The proportion of the population and number of individuals sampled strongly influences the power to detect sex-biased dispersal, thus our limited sampling on Oʻahu and Kauaʻi may not be sufficient to detect true patterns within those islands.

Although radio tracking studies on Hawaiʻi and Maui have demonstrated varying scales of nightly movement for both sexes of ʻōpeʻapeʻa (Bonaccorso et al., 2015; HT Harvey and Associates, 2020), no telemetry studies have demonstrated seasonal or long-range movements within or between islands. Seasonal movements have thus far have been inferred from acoustic studies across elevational gradients without a sex specific context for interpretation (Menard, 2001; Gorresen et al., 2013; Bonaccorso et al., 2016; Todd, Pinzari & Bonaccorso, 2016; West Inc, 2020; West Inc, 2021). The contrasting signals of possible female sex-biased dispersal in our results, are difficult to untangle from sex-biased dispersal because of testing limitations and panmixia effects within islands and warrant future investigation. Continued field-based capture studies focusing on female roost fidelity, the identification of seasonal roosts, swarming (mating) sites, and examination of genotypes at more distant collection sites within islands could contribute to a greater understanding of mating systems, gene flow, and sex-specific dispersal patterns in ʻōpeʻapeʻa.

Estimates of genetic effective population size

Extinction risk can be predicted if both population census size, genetic effective population size, and levels of gene flow are known (Luikart et al., 2010). Although current census population sizes of ʻōpeʻapeʻa on each Hawaiian island remain unknown and confidence intervals are wide, estimates of contemporary genetic effective population sizes from this study indicate that the Oʻahu and Maui populations may be lower in abundance than those of Hawaiʻi and Kauaʻi. Effective population size estimates aid in evaluating conservation programs for endangered species as they may flag inbreeding depression and potential loss of adapability to environmental change and disease (Waples, 2002; Funk et al., 2019; Hohenlohe, Funk & Rajora, 2021). Methods that infer effective population size assume simple population models; however, migration, spatial structure, reproductive success, and overlapping generations can confound results (Waples & England, 2011). Effective population estimates are not equivalent to the true “census” population size, and all estimates we report here should be interpreted with caution and subject to change with additional data. Thresholds for population viability, known as the 100/1000 rule (Franklin, 1980; Frankham, Bradshaw & Brook, 2014), state that a population should have an N_e of at least 100 individuals to avoid inbreeding depression, and an N_e of at least 1,000 to avoid erosion in evolutionary potential. Using the linkage-disequilibrium (LD) method, contemporary estimates of genetic effective population size (N_eC) vary by island, with the Hawai′i estimate exceeding 1,000 individuals. On Hawaiʻi and Kauaʻi estimates were imprecise having infinity as an upper bound in confidence, and indicating N_eC could be greater than 500, as the LD method is known to be more precise when effective population size is <500 (Gilbert & Whitlock, 2015). We estimated N_eC with bounded CI to be <500 on Maui and <40 on Oʻahu, thus management actions that reduce mortality from threats, such as collisions with wind turbines, are important to the long-term persistence of hoary bat populations on those islands (Frick et al., 2017; Friedenberg & Frick, 2021). Reductions in genetic diversity caused by population decline, disconnection, or increased inbreeding, to either the highly diverse Maui population or the differentiated Oʻahu and Kauaʻi populations could erode the total genetic diversity and adaptive potential of ʻōpeʻapeʻa across the Hawaiian Islands.

Statistical estimation techniques of contemporary effective population sizes using the LD method are prone to uncertainty, in species such as bats with overlapping generations (Waples, Antao & Luikart, 2014). However, Waples & Do (2010) suggest that estimates using LD on population samples that include multiple cohorts, as we have presented, may approximate N_eC, if the population sample included enough cohorts to encompass one generation. Although generation time for ʻōpeʻapeʻa is unknown, we assume generation time to be at least two years whereas our population samples spanned nearly a decade. An alternative model for species with overlapping generations is to use a single-sample method based on parentage (Wang, 2009), which requires the age and sex for each individual. New DNA methylation techniques show promise for estimating generation time and age estimation in bats (Wilkinson et al., 2021). Additionally, utilizing single nucleotide polymorphisms from a larger number of individuals may refine estimates of ʻōpeʻapeʻa effective population sizes, as was recently examined for L. cinerus by Cornman et al. (2021).

Long-term effective population estimates for all ʻōpeʻapeʻa sampled in this study produced a range of 248 to 62,162 bats. We reported both SMM and IAM models using high, average, and low mutation rates because evolutionary mutation rates can vary within and across loci. These averaged models produced similar estimates for each island and indicate that long term maintenance of genetically diverse populations may require several hundred to tens of thousands of ʻōpeʻapeʻa. Historical female genetic effective population estimates for ʻōpeʻapeʻa presented wider ranges of 1,384 to 286,658 bats; however, the precision of these estimates is unknown because confidence intervals cannot be calculated. Molecular population genetic estimates of genetic effective population size between historical female (N_Ef) and long-term N_e overlap when confidence intervals of the latter are taken into consideration. Historical genetic effective population sizes likely reflect the ancestral high genetic diversity of founding L. cinereus individuals from a large continental population.

Demographic changes

Population bottlenecks are distinguished by loss of rare alleles, which typically has little effect on heterozygosity. When the population size or genetic effective size is reduced, inbreeding leads to an excess of common microsatellite alleles compared to the number of rare alleles expected under equilibrium conditions (Cornuet & Luikart, 1996). The M-ratio test indicated genetic bottlenecks with different intensity may have occurred across Maui, Oʻahu, and Kauaʻi. These measures do not detect or reflect possible recent declines in the abundance of ʻōpeʻapeʻa because genetic diversity measures lack the power to detect changes over short time spans for species with large effective population sizes, high pre-bottleneck diversity, or high connectivity (Peery et al., 2012). Genetic bottlenecks also can retain the signature of decline over many generations and be a product of population disconnection (e.g., establishing a population on a new island) rather than exclusive of reduction in population size or effective size (Broquet et al., 2010). Population growth can be signaled by significant negative Tajima’s D values, high haplotype diversity, and low sequence diversity. Tajima’s D was significantly negative for Hawaiʻi, indicating past signatures of population expansion and is congruent with our migration history, gene flow, and haplotype results. Signatures of population expansions using coalescent modeling were detected in ʻōpeʻapeʻa mitochondrial clades at 800, 10,000 (Russell et al., 2015), and 20,000 years before present (Baird et al., 2017), consistent with population expansion on Hawaiʻi and supported by the high density of minor alleles in site frequency spectra from Pinzari et al. (2020). The power to detect bottlenecks would be greatly improved with demographic coalescent methods as our current results do not offer a context on when in the past a bottleneck occurred for each the island. Incorporating samples from ʻōpeʻapeʻa in museum collections prior to 2005 could enable more accurate evaluations of bottlenecks and/or declines in genetic diversity that may have occurred in the last two centuries, as well as refine estimates of historical population sizes.