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.

Taxon and Gene Sampling
We obtained tissue samples from 168 Pteropus individuals including 34 currently recognized species or putative species out of 65 valid, extant species (Appendix S1). These samples included 10 of the 12 Pteropus species groups recognized by Almeida et al. (2014) and 1 species incertae sedis. The omitted species group include an African lineage of two species closely allied to the vampyrus species group (Almeida et al. 2014), and a group comprising a single New Caledonian species. We initially included multiple individuals of widespread Pteropus alecto, P. hypomelanus, and P. vampyrus 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 Pteropus’ sister genus Acerodon (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 Query function in the Species Identifier module of TaxonDNA 1.0 (Meier, Shiyang, Vaidya, & Ng, 2006) and the View as pairwise distances 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.

Phylogenetic Analyses
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 Acerodon celebensis 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 Acerodon-Pteropus genus split (8.01 ± 1 Mya), the emergence of the P. gilliardorum+P. woodfordi clade (0.93 ± 0.1 Mya), and the emergence of the P. pelewensis-P. yapensis 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 Pteropus species and 3 outgroup taxa with one tip per taxon was used in subsequent biogeographic analyses.

Biogeographic Analyses
We inferred the biogeographic history of Pteropus using BioGeoBEARS 1.1.1 (Matzke, 2014) and allowed the d (dispersal) and e (extinction) parameters to vary freely. Results from the DEC, DEC+j, DEC*, DEC*+j and four other models implemented in BioGeoBEARS were compared using AICc to determine the effect of adding the j 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 j 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 Pteropus 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 (e.g., 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 Pteropus, and the supertree was subsequently omitted from our figures for clarity.