DOI: 10.5061/dryad.8892g0m
title: Data from: Dispersal out of Wallacea spurs diversification of Pteropus flying foxes, the world’s largest bats (Mammalia: Chiroptera)

description: Aim: Islands provide opportunities for isolation and speciation.
 Many landmasses in the Indo-Australian Archipelago (IAA) are oceanic
 islands, and founder-event speciation is expected to be the predominant
 form of speciation of volant taxa on these islands. We studied the
 biogeographic history of flying foxes, a group with many endemic species
 and a predilection for islands, to test this hypothesis and infer the
 biogeographic origin of the group. Location: Australasia,
 Indo-Australian Archipelago, Madagascar, Pacific Islands Taxon: Pteropus
 (Pteropodidae) Methods: To infer the biogeographic history of
 Pteropus, we sequenced up to 6169 bp of genetic data from 10 markers and
 reconstructed a multilocus species tree of 34 currently recognized
 Pteropus species and subspecies with 3 Acerodon outgroups using BEAST and
 subsequently estimated ancestral areas using models implemented in
 BioGeoBEARS. Results: Species-level resolution was occasionally
 low because of slow rates of molecular evolution and/or recent
 divergences. Older divergences, however, were more strongly supported and
 allow the evolutionary history of the group to be inferred. The genus
 diverged in Wallacea from its common ancestor with Acerodon; founder-event
 speciation out of Wallacea was a common inference. Pteropus species in
 Micronesia and the western Indian Ocean were also inferred to result from
 founder-event speciation. Main conclusions: Dispersal between
 regions of the IAA and the islands found therein fostered diversification
 of Pteropus throughout the IAA and beyond. Dispersal in Pteropus is far
 higher than in most other volant taxa studied to date, highlighting the
 importance of inter-island movement in the biogeographic history of this
 large clade of large bats.
description: <b>Taxon and Gene Sampling</b><br> We obtained
 tissue samples from 168 <i>Pteropus</i> individuals including
 34 currently recognized species or putative species out of 65 valid,
 extant species (Appendix S1). These samples included 10 of the 12
 <i>Pteropus </i>species groups recognized by Almeida et al.
 (2014) and 1 species <i>incertae sedis</i>. The omitted
 species group include an African lineage of two species closely allied to
 the <i>vampyrus</i> species group (Almeida et al. 2014), and a
 group comprising a single New Caledonian species. We initially included
 multiple individuals of widespread <i>Pteropus alecto, P.
 hypomelanus,</i> and <i>P. vampyrus </i>and treated them
 as distinct OTUs to examine the hypothesis that different populations
 might actually be cryptic species complexes (Bickford et al., 2007). Some
 of these OTUs correspond to different subspecies or previously recognized
 full species (Andersen, 1912; Corbet & Hill, 1992). We initially
 included data from 19 outgroup samples comprising 11 species in 7
 pteropodid genera, but this diverse outgroup sample dramatically increased
 analysis time and reduced branch support. We therefore trimmed the
 outgroups to three species of <i>Pteropus’</i> sister genus
 <i>Acerodon</i> (Appendix S1; Almeida et al., 2011). We
 attempted to amplify 8 nuclear and 2 mitochondrial loci from each sample,
 totalling up to 6169 bp of sequence data per specimen. Gene names, PCR
 primers, and thermal cycling conditions are provided in Appendix S2.
 Successfully amplified PCR products were cleaned with ExoSAP or a vacuum
 manifold through a Millipore filter plate. Products were sequenced on an
 Applied Biosystems 3730xl DNA Analyzer. All markers were aligned using
 Geneious 5.4.3 implementing MAFFT 7.0
 (Katoh & Standley, 2013) using the G-INS-i option. We used the
 <i>Query</i> function in the Species Identifier module of
 TaxonDNA 1.0
 (Meier, Shiyang, Vaidya, & Ng, 2006) and the <i>View as
 pairwise distances</i> option in SequenceMatrix
 1.8 (Vaidya, Lohman, & Meier, 2011) to find identical or
 suspiciously similar heterospecific gene sequences, then re-amplified and
 re-sequenced sequenced them to confirm that sequences were not attributed
 to the wrong species through contamination or other error.
 <b>Phylogenetic Analyses</b><br> Haplotypes for
 each nuclear marker were estimated using PHASE 2.0
 (Gowri-Shankar & Rattray, 2007) using 10 runs of 20,000
 iterations. The program jModelTest
 2.1.10 (Darriba, Taboada, Doallo, & Posada, 2012) was used to
 select appropriate substitution models using the corrected Akaike
 Information Criterion (AICc; Appendix S2). Individual gene trees were
 reconstructed using MrBayes 3.2.3
 (Ronquist & Huelsenbeck, 2003) and the analysis was run for 10
 million generations with 25% burn-in and <i>Acerodon
 celebensis</i> designated as the sister outgroup (Almeida et al.
 2011). To minimize noise and shorten analysis time for species tree
 inference, a maximum likelihood (ML) consensus tree was first inferred
 using IQ-TREE 1.6.10
 (Nguyen, Schmidt, von Haeseler, & Minh, 2015) on a dataset from
 which the uninformative loci BDNF, PLCB4, RAG-1, and RAG-2 were removed. A
 species tree was then estimated with BEAST 2.5.2
 (Bouckaert et al., 2019) using a log-normal relaxed clock birth-death
 model and the consensus ML tree as a starting guide tree. A Yule prior was
 not a better fit to our data, thus we chose the birth-model, which is
 better suited as a null hypothesis species diversification. Given the
 paucity of bat fossils, we used median divergence times and confidence
 intervals of secondary calibrations from Almeida et al. (2014) to
 calibrate our tree, and utilized the same divergences where possible: the
 <i>Acerodon</i>-<i>Pteropus</i> genus split (8.01
 ± 1 Mya), the emergence of the <i>P. gilliardorum+P.
 woodfordi</i> clade (0.93 ± 0.1 Mya), and the emergence of the
 <i>P. pelewensis</i>-<i>P. yapensis </i>clade
 (0.03 ± 0.01 Mya). Markov chain performance for species trees was verified
 using Tracer to ensure ESS values above 200, appropriate burn-in, and
 chain convergence. This required using LogCombiner
 2.5.2 (Bouckaert et al., 2019) to combine multiple BEAST runs with
 the same parameters with 20% burn-in; the final sample of trees was run
 for 200 million generations. There were 14401 trees in the final posterior
 distribution from BEAST. The
 resulting species tree of 34 <i>Pteropus </i>species and 3
 outgroup taxa with one tip per taxon was used in subsequent biogeographic
 analyses. <b>Biogeographic
 Analyses</b><br> We inferred the biogeographic history of
 <i>Pteropus</i> using BioGeoBEARS
 1.1.1 (Matzke, 2014) and allowed the <i>d </i>(dispersal) and
 <i>e </i>(extinction) parameters to vary freely. Results from
 the DEC, DEC+<i>j</i>, DEC*, DEC*+<i>j</i> and
 four other models implemented in BioGeoBEARS
 were compared using AICc to determine the effect of adding the <i>j
 </i>parameter on model performance. DEC* refers to models where
 transitions into the null range in the transition rate matrix were
 prohibited, resulting in better inference of local extinction (Massana,
 Beaulieu, Matzke, & O’Meara 2015). While we acknowledge criticism
 of the <i>j </i>parameter by Ree and Sanmartín (2018),
 peripatric (“jump”) speciation is often the most biologically realistic
 scenario for colonization and subsequent divergence on remote, oceanic
 islands. For biogeographic analyses, we divided the
 range of <i>Pteropus</i> into six zoogeographic areas that
 were defined based on geologic history (Hall 2002) and known biogeographic
 barriers: Wallacea (A), New Guinea + Australia + South Pacific (B),
 Sundaland + South Asia (C), the Philippines (D), Micronesia (E), and
 western Indian Ocean (F) (Fig. 1). The broad geographic scope of this
 study required grouping many small islands into a single region for
 computational efficiency and to achieve our goal of assessing the
 biogeographic history of the genus as a whole. Improbable area
 combinations were excluded from analysis to model more realistically the
 spatial distribution of biogeographic areas (<i>e.g.</i>,
 direct dispersal between Wallacea and the western Indian Ocean without
 inhabiting areas in between was considered statistically unlikely). The
 maximum number of ancestral areas was set to four, as this was the maximum
 number of areas occupied by any OTU in our analysis. Species were
 considered “widespread” if their range encompasses three or more
 biogeographic areas. Area assignments for each species were based on known
 distribution data (Simmons, 2005). To account for
 poorly supported nodes in our phylogeny, we ran BioGeoBEARS on
 several possible trees to determine which biogeographic results were
 recovered repeatedly. To ensure that our taxon sample did not bias
 biogeographic inferences, the tree that we inferred was grafted to the
 pteropodid tree of Almeida et al. (2011) using the supertree Maximum
 Representation with Parsimony method (Sanderson, Purvis, & Henze,
 1998). This larger set of outgroups did not change ancestral area
 estimates for <i>Pteropus</i>, and the supertree was
 subsequently omitted from our figures for clarity.