Genetic differentiation and demographic trajectory of the insular Formosan and Orii’s flying foxes

Insular flying foxes are keystone species in island ecosystems due to their critical roles in plant pollination and seed dispersal. These species are vulnerable to population decline because of their small populations and low reproductive rates. The Formosan flying fox (Pteropus dasymallus formosus) is one of the five subspecies of the Ryukyu flying fox. P. d. formosus has suffered from a severe decline and is currently recognized as a critically endangered population in Taiwan. On the contrary, the Orii’s flying fox (P. d. inopinatus) is a relatively stable population inhabiting Okinawa Island. Here, we applied a genomic approach called double digest restriction-site associated DNA sequencing to study these two subspecies for a total of seven individuals. We detected significant genetic structure between the two populations. Despite their contrasting contemporary population sizes, both populations harbor very low degrees of genetic diversity. We further inferred their demographic history based on the joint folded site frequency spectrum and revealed that both P. d. formosus and P. d. inopinatus had maintained small population sizes for a long period of time after their divergence. Recently, these populations experienced distinct trajectories of demographic changes. While P. d. formosus suffered from a drastic ~10-fold population decline not long ago, P. d. inopinatus underwent a ~4.5-fold population expansion. Our results suggest separate conservation management for the two populations—population recovery is urgently needed for P. d. formosus while long-term monitoring for adverse genetic effects should be considered for P. d. inopinatus.

Methods
Seven individuals of Pteropus dasymallus bats were collected. Among them, three are P. d. formosus from Taiwan (designated as 178, 617, and 690), while four are P. d. inopinatus from Okinawa Island (designated as R10, R36, R4, and R5). For the three P. d. formosus individuals, frozen muscle tissues were obtained from the Wildlife Cryobank of Taipei Zoo. For each of the four P. d. inopinatus individuals, skin samples were excised from the wing membranes using a biopsy punch (3 mm in diameter); eight punches were collected per individual. All biopsies were preserved in Allprotect Tissue Reagent (Qiagen) and stored at -20˚C. Genomic DNA was extracted from samples using the DNeasy Blood and Tissue kit (Qiagen) following the manufacturer’s instructions.

Double digest restriction-site associated DNA sequencing (ddRADSeq) method was applied to these samples. For a given individual genome, two restriction enzymes, EcoRI and MseI, were used to cleave the DNA. The cleaved DNA fragments were ligated with adaptors that contain barcodes and sequences complementary to the PCR primers. After these ligated fragments were amplified by PCR, fragments with lengths between 250 base pairs (bp) and 450 bp were retained and sent for paired-end sequencing. Sequencing was carried out using the Illumina HiSeqTM 2000 sequencing system in the core facility center of Genomics company in New Taipei City.

The demultiplexed paired-end sequence reads were processed with SOAPnuke2 to discard low-quality reads and unwanted sequences. First, reads with correct barcodes were retained exclusively and their barcodes were trimmed. Reads without correct restriction enzyme cleavage sites were subsequently removed from the dataset. Low-quality reads were also removed if the proportions of base pairs with Phred quality scores lower than 15 were >40% or the proportions of unidentified base pairs were >10%. Next, we removed the restriction enzyme cleavage sites from all remaining reads. The remaining forward reads were trimmed to 87 bp to accommodate the differences in sequence lengths caused by barcode trimmings and restriction enzyme cleavage site removals.

Finally, we assembled the clean reads using Stacks 2.2 with the following assembly parameters: m = 3, M = 1, n = 1. Secondary read incorporations were disabled and the deleveraging algorithm was enabled. Only sites that were sequenced in all individuals were assembled (r100) to prevent systematic biases on population genetics statistics including estimates of nucleotide diversity and Tajima's D. Since SNPs with excessive heterozygosity across all individuals can result from misaligning reads of paralogous or repetitive loci, we applied the Hardy-Weinberg equilibrium (HWE) test to exclude the assembled reads or SNPs that showed significant deviations from HWE (p < 0.05). Two types of heterozygosity filters were applied to the dataset to generate datasets for population genetics analyses that rely on accurate nucleotide diversity estimates or only genetically unlinked SNPs.