DOI: 10.5061/dryad.msbcc2fxh
title: Genetic differentiation and demographic trajectory of the insular Formosan and Orii’s flying foxes

description: 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.
description: Seven individuals of <em>Pteropus dasymallus</em>
 bats were collected. Among them, three are <i>P. d.
 formosus</i> from Taiwan (designated as 178, 617, and 690), while
 four are <i>P. d.</i> <i>inopinatus</i> from
 Okinawa Island (designated as R10, R36, R4, and R5). For the three
 <i>P. d. formosus</i> individuals,<i> </i>frozen
 muscle tissues were obtained from the Wildlife Cryobank of Taipei Zoo. For
 each of the four <i>P. d. inopinatus</i> 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
 HiSeq<sup>TM</sup> 2000 sequencing system in the core facility
 center of Genomics company in New Taipei City.
 The demultiplexed paired-end sequence reads were processed with
 <i>SOAPnuke2</i> 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 &gt;40% or the proportions of
 unidentified base pairs were &gt;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
 <i>Stacks</i> 2.2 with the following assembly parameters:
 <em>m</em> = 3, <em>M</em> = 1,
 <em>n</em> = 1. Secondary read incorporations were disabled
 and the deleveraging algorithm was enabled. Only sites that were sequenced
 in all individuals were assembled (<em>r</em>100) 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 (<i>p</i> &lt; 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.
description: <strong>RAD loci from all individuals
 (Pd-pooled_RADloci.fa)</strong><br> The assembled sequences of
 all RAD loci (191,263) in each individual are pooled in this FASTA file.
 These RAD loci have not been filtered by the HWE test but have passed
 through other quality control steps. Unique RAD loci are characterized by
 the serial number in the header line (CLocus_XX). Individual ID can also
 be found in the brackets at the end of each header. If the RAD locus is
 identified as homozygous with no variants, there will only be one sequence
 per individual. Otherwise, two alleles (labeled with Allele_0 and Allele_1
 in headers) can be found with the same RAD loci and individual labels.
 Unrecognized nucleotides (N) are included in these sequences.<br>
   <strong>SNPs after heterozygosity filters
 (Pd-SNPs-hf*.map and Pd-SNPs-hf*.ped)</strong><br> The SNPs
 used in the population genetics analyses in our study are provided in
 PLINK .map and .ped files. They were filtered with two types of HWE
 filters. The file suffix -hf implies that SNPs showing significant
 deviation from HWE were removed from this dataset, while the suffix
 -hf_loci implies that whole RAD loci were removed if any SNPs with
 deviation from HWE were present on them. The former dataset was used in
 analyses that require only SNPs (e.g. STRUCTURE, PCA), while the latter
 was used in analyses that need the estimates of nucleotide diversities
 (e.g. demographic inference). The numbers of SNPs provided in the file are
 33,005 (Pd-SNPs-hf) and 32,935 (Pd-SNPs-hf_loci).<br>  
 <strong>The lengths of RAD loci
 (Pd-loci_lengths.tsv)</strong><br> This tsv file records the
 lengths of each unique RAD loci, with or without unidentified nucleotides
 (N). This file was used, along with the above SNP files, in estimating the
 nucleotide diversity in <em>P. .d formosus</em> and
 <em>P. d. inopinatus</em>. All loci have passed through the
 quality controls including the hf_loci filter. The lengths for a total of
 188,588 loci are provided in this file.