Northward geographic diversification of a kleptoparasitic spider Argyrodes lanyuensis (Araneae, Theridiidae) from the Philippine Archipelago to Orchid Island

Abstract Oceanic islands are unique geographic systems that promote local adaptations and allopatric speciation in many of their highly endemic taxa. This is a common case in the Philippine Archipelago, where numerous unrelated taxa on islands have been inferred to have diversified in isolation. However, few cases have been reported in invertebrates especially among parasitic organisms. Here, we tested for biogeographical structure in novel populations of the “generalist" kleptoparasitic spider, Argyrodes lanyuensis Yoshida, Tso & Severinghaus, 1998 in the Philippines. Results showed that, in addition to Orchid/Lanyu Island, this species has a wide geographic distribution in the Philippine Archipelago. The estimated divergence time of this lineage using the mitochondrial cytochrome oxidase 1 (mt‐CO1) suggests that this species diverged ca 3.12 MYA, during the Pliocene. Two reciprocal monophyletic clades were elucidated in A. lanyuensis, but with limited differentiation across Pleistocene Aggregate Island Complex (PAIC) boundaries and modern‐day islands. However, in our analyses of morphological variation, we identified two phenotypically differentiated units in males (Orchid Island, Taiwan + Luzon, Philippine PAIC populations vs. Palawan + West Visayan + Mindanao PAIC populations). We infer that this species diverged in the southern portion of the Philippine Archipelago and only recently colonized Orchid Island. Our study provides new information on the extensive distribution of A. lanyuensis outside Orchid Island, Taiwan, but we documented a very limited geographically associated genetic variation. Our study points to behavioral phenomena such as foraging behavior as essential contributor to the evolutionary process of species diversification, in contrast to the traditionally invoked geographic drivers of divergence.


| INTRODUC TI ON
Oceanic island chains usually host high levels of endemic terrestrial biodiversity because of strong geographic isolation, which promotes the partitioning of their fauna and flora (Gillespie, 2007;Lomolino et al., 2010). Dispersal plays an important mechanism in the process of diversification of taxa in an oceanic island (De Queiroz, 2005;Gillespie et al., 2012). Once a historical oceanic island emerged above the surface of the ocean, it is then available for colonization of taxa from distant land areas (Cowie & Holland, 2006). However, effective colonization to oceanic islands particularly by terrestrial biota depends on several factors, such as climatic conditions, wind speed variation, local adaptation, sizes of islands, distance to source biota, and geographic boundary fluctuations. (Leihy & Chown, 2020;Lomolino et al., 2010;Whittaker et al., 2007). All of which may contribute to the historical subdivision of populations on oceanic islands.
The Philippine Archipelago is known as one such highly partitioned case: a dynamic, highly fragmented geographical template.
It consists of more than 7,000 oceanic islands situated at a unique location-spanning portions of the Australasian and Asian faunal regions (Brown & Diesmos, 2009;Brown et al., 2013;Lohman et al., 2011). It hosts substantial genetic structure, both within species and among highly differentiated lineages Hosner et al., 2014;Siler, Oaks, et al., 2012;Wood et al., 2020). The subdivision of populations, species, and even higher taxa have been hypothesized to be the result of dynamics current and historical geographic processes of the archipelago (Hall, 1998(Hall, , 2002Yumul et al., 2009). With the relatively clear understanding of the geographic boundaries, and dynamic nature of their corresponding geological history, reassessment of species diversity and mechanisms of diversification has been explored comprehensively in multiple clades (Hosner et al., 2014;Linkem et al., 2010;Siler, Jones, et al., 2012;Weinell & Brown, 2017). This resulted in the identification of localized evolutionary trends and many instances of allopatric speciation following bouts of dispersal (Barley et al., 2020;Brown et al., 2016;Oaks et al., 2019;Siler et al., 2010). It is intuitive to consider that pronounced subdivision of the Philippine Islands might cause or be related to diversification, presumably resulting in the formation of new endemic species once their ancestors invaded relatively isolated islands (Heaney, 1985;Inger, 1954). However, whether such species continued to expand their range via recent dispersals among islands has rarely been reported (but see Brown et al., 2010;Siler et al., 2014).
The Philippines is located approximately 390 km south of Taiwan, but Orchid/Lanyu Island (Taiwan) and the Batanes and Babuyan island groups (Philippines) span the intervening seas with a series of small island chains (Figure 1a). Initially documented on Orchid/Lanyu Island, the kleptoparasitic spider, Argyrodes lanyuensis ( Figure 1b), has been considered endemic to this small island since it was described in 1998 (Yoshida et al., 1998). However, our recent sampling of argyrodinae spiders in the Philippines has revealed the occurrence of A. lanyuensis in at least six of the archipelago's islands ( Figure 1a). We used this species to assess whether a strongly subdivided geographic system (the oceanic portion of the archipelago) would be effective in generating pronounced geographical structure in genetic variation among populations of this kleptoparasitic spider.
The foraging behavior of the subfamily Argyrodinae is remarkable in that they rely on either araneophagy or kleptoparasitism-or sometimes both-as their main feeding strategy (Cobbold & Su, 2010;Vollrath, 1979;Whitehouse, 2011). Argyrodes lanyuensis is closely related to the Philippine endemic A. tripunctatus Simon 1877, and two Australasian species, A. nasutus Pickard-Cambridge 1880, and A. rainbow Roewer 1942(Su & Smith, 2014. Based on the first reports of Yoshida et al. (1998), A. lanyuensis forages prey items and silk from the webs of a wide range of orb-weaving spider hosts, that is, Nephila, Gasteracantha, and Cyrtophora. Aside from orb-weaving hosts, it was also observed to hunt prey items and consume silk on Achaearanea (Theridiidae) host. Thus, an ecological "generalist" kleptoparasites like A. lanyuensis tend to have high tolerance on a wide array of spider hosts than ecological "specialist" kleptoparasites which specifically utilize one species/genus of orb-weaving hosts (e.g., A. fissifrons and A. miniaceus kleptoparasites). Since most of the geographic variability, biogeography, and individual species distributions of Argyrodinae on oceanic islands have not been fully characterized, we focused on A. lanyuensis as a fitting representative of ecological "generalist" kleptoparasitic spiders to distinguish it from "specialist" kleptoparasites.
The Pleistocene Aggregate Island Complex (PAIC) model of speciation (Inger, 1954;Heaney, 1985Heaney, , 1986review: Brown & Diesmos, 2002 has been used as an operational hypothesis to generate testable predictions related to the analysis of diversification patterns among Philippine biota (Evans et al., 2003;Sánchez-González & Moyle, 2011;. The Pleistocene glacial cycles (between 2.5 MYA to 18 KYA) resulted in the repeated rising and lowering of sea levels (100-140 m). In the Philippines, this led to the repetitive isolation and formation of land bridges between neighboring islands separated by shallow seas (Figure 1a). With the tracing of bathymetric contours (100-140 m) within this period, Pleistocene islands can be estimated with the maximum extent of land bridges. This resulted in six major larger island-amalgamations known as PAICs: Luzon, Mindanao, Western Visayas, Mindoro, Sulu, and Palawan (Brown & Diesmos, 2002;Heaney, 1985Heaney, , 1986. These paleoisland connections among islands in the Philippines served as a basis for predicting patterns of species diversity and distribution. To date, several vertebrate taxa like mammals, lizards, frogs, and birdsshowed nearly complete concordance to PAIC boundaries (Evans et al., 2003;Heaney, 1985Heaney, , 1986McGuire & Alcala, 2000;Sánchez-González & Moyle, 2011). However, applying the PAIC speciation model to highly dispersive arthropod species is sparse in literature, except for one pilot study . Even though it has not been utilized more often to terrestrial invertebrate species due to the characteristic of flight and ballooning, it is also worth noting that this speciation model has been used to explain the diversification patterns of widely distributed volant mammals and birds (Heaney et al., 2005; Sánchez-González & Moyle, 2011).
The predictions derived from a strict interpretation of the PAIC Paradigm would include (1) a homogenized (or nearly so) gene pool of island populations within PAICs and (2) limited gene flow, leading to pronounced geographical structure, among and between PAICs. It follows, then, that if a particular taxon colonized the archipelago before or during the Pleistocene, the distribution of its species (or populations) would likely be found today in concordance with the PAIC model's six major faunal regions. The Philippines Archipelago has a dynamic geologic history (Hall, 2002;Yumul et al., 2003Yumul et al., , 2009, which likely influenced the diversification of its fauna and flora (Brown & Diesmos, 2009;Brown et al., 2013). Therefore, we assumed that heterogeneous, interrupted, and partitioned geographic template of land area throughout the archipelago might have led to distinct populations of A. lanyuensis across oceanic islands including the island banks stretching north toward Taiwan and Orchid Island (Figure 1c). However, if we consider the dispersal ability of spiders through long-distance ballooning (Bell et al., 2005;Bishop & Riechert, 1990), then we would expect to see little to no differentiation of A. lanyuensis populations, as it would greatly affect the gene flow of this species. Additionally, the behavior of this species, which is a generalist kleptoparasite, would also explain a little to no differentiation of taxa because generalists do not need to specifically adapt to a particular host (Su et al., 2018).
To ascertain how A. lanyuensis may have dispersed and colonized in the Philippine Archipelago and Orchid Island, we first update its geographical distribution and used time-calibrated phylogenetic analyses. We infer the ancestral area range evolution using biogeographical reconstruction models. Initially, we hypothesized that this species diverged from Sundaland and colonized Philippine islands via Eastern and Western arcs ( Figure 1a; Route 1, 2), through the northern-most islands (Babuyan and Batanes island groups; Route 3), and eventually colonized Orchid Island, as suggested by the results of Su and Smith (2014). Alternatively, if the current distribution of the species came about by recent southward colonization (<1MYA), the species may have originated on Orchid Island (Route 4; Figure 1a) and subsequently colonized the Philippines via the Taiwan-Batanes-Babuyan island chain (Dickerson, 1928;Esselstyn & Oliveros, 2010;Oliveros et al., 2011). Thus, we undertook the current study to test the north-to-south versus south-to-north predictions derived from PAIC and analyze dispersal or vicariance events. F I G U R E 1 Paleogeological features of Taiwan and the Philippine Archipelago (Hall, 1998). (a) Arrangement of Taiwan-Philippine archipelago during the Pleistocene epoch, <2 Million Years ago (MYA) with the −120 m contours (gray) of the Pleistocene aggregate island complexes (PAICs). Taiwan was connected with mainland China, while Philippines formed different PAICs as indicated by each color. We hypothesized south-to-north colonization via ① eastern island arc hypothesis, ② western island arc hypothesis, and ③ Batanes-Babuyan islands route; and north-to-south colonization via ④ Taiwan route. (b) Photograph of male Argyrodes lanyuensis Yoshida et al., 1998 2 | MATERIAL S AND ME THODS

| Taxon sampling
We collected A. lanyuensis samples from the main islands in the Philippines and Taiwan between 2005 and 2007, and from July to August 2019. We found A. lanyuensis at only 13 collection sites on six islands: Orchid Island, Luzon, Palawan, Negros, Panay, and Mindanao ( Figure 1a). Samples were collected from the webs of orb-spinning spiders of families Araneidae, Tetragnathidae, and Uloboridae. We preserved the specimens in 95% ethanol and stored at −30°C, for subsequent morphological examination and DNA extraction. All specimens were deposited in Evolution and Ecological Genomics (EEG) Laboratory, Kaohsiung Medical University, Kaohsiung, Taiwan.
Specific collection information and sample accession number for each specimen are reported in Appendix 1.

| Morphological variation
To assess the geographic variation of A. lanyuensis populations from the Philippines and Orchid Island, we examined adult specimens for variation in continuous morphometric measurements. Male (n = 37) and female (n = 38) samples were observed under a Leica stereomicroscope. We embedded each specimen in a gel-loaded calibration slide (1 division = 0.1 mm; 1 division = 0.01) and used tethered Nikon camera D5600 to capture high-resolution images (Appendix 5). We utilized measure3 software (Tsai, 2021; https://github.com/yucen wan/Spide r-measure) to generate calibrated measurements of specific body characteristics from captured high-resolution images. We normalized body morphometrics using carapace length, following character definitions of Yoshida et al. (1998). The measured body characteristics include total length (TL), carapace length (CL), carapace width (CW), total length of each leg (L1TOT; L2TOT; L3TOT, and L4TOT)), and the length of each leg (I-IV) segment: femur (L1F-L4F); tibia + patella (L1PT-L4PT); metatarsus (L1M-L4M), and tarsus (L1T-L4T). We additionally measured palp morphometrics from the male specimens, which include total palp length (PL), bulb length (BL), median apophysis (MA), accessorial apophysis (AP), and embolus length (EL). Bulb length was used to normalize all the palp morphometrics. All body and palp measurements used in this study were displayed in Appendix 5. Variation in morphometric dimensions (separately for males vs. females) was summarized in Principal Component Analysis (PCA) using the "prcomp" function in R 3.6.1 (R Core Team, 2017). Data visualization was carried out using the R package ggfortify (Horikoshi & Tang, 2018). We used nonlinear iterative partial least squares (NIPALS, followed Wold, 1973) in which the algorithm conducts local regressions using the latent components to predict and impute missing values caused by poor preservation conditions (Female, n = 9; 1.03% of the data matrix; Male, n = 50; 4.83% of the data matrix). To avoid multicollinearity problems among the measurements of our morphological data, we followed Vignon (2011) to adopt the Partial Least Square-Discriminant Analysis (PLS-DA), assessing if individuals clustered into geographical distributions based on morphology. We used the "plsda" function within the R package mixOmics (Rohart et al., 2017), where all measurements were included as response variables. Permutational test with 9,999 repetitions was performed based on cross-model validation procedures, where estimation of the classification error rate (CER) was used as the test statistics. Additionally, the function "pairwise.
MVA.test" in the same R package was implemented for pairwise comparisons of clusters.

| DNA extraction, marker choice, and PCR amplification
We extracted the genomic DNA from legs and prosomal tissues of preadult and adult specimens following the Maxwell® RSC Blood DNA Kit AS1400 protocol. Tissues were homogenized in 300 μl Lysis Buffer and 30 μl Proteinase K (PK) Solution and incubated at 56°C for 2 hr. We purified the genomic DNA through the Maxwell ® RSC Instrument following the manufacturer's instructions. The extracted genomic DNA was stored at −30°C condition until used for polymerase chain reaction (PCR) amplification.
We sequenced the mitochondrial cytochrome oxidase I (CO1) partial gene region, which is an effective genetic marker in species identification and taxonomic delimitation (Hebert et al., 2004), especially for invertebrates (Cao et al., 2016;Carew et al., 2007;Gutiérrez et al., 2014). The CO1 fragment was targeted and amplified using primer pairs, CO1-F and CO1-r designed by Su and Smith (2014). PCR amplification was performed in a TurboCycler 2 thermal cycler (TCST-9622, Taiwan) with a total volume of 25 μl with 12 μl of premix, 10 μl of nuclease-free water, and 0.5 μl to each of the primers. PCR products were visualized through 1.5% agarose gel electrophoresis to check amplified DNA fragments of the expected size and sequenced at the genetic sequencing facility of Genomics Co. Ltd., Taiwan.

| Sequence alignment and molecular data analysis
We filtered all the sequences according to the quality control reports and obtained a total of 95 CO1 sequences. Some samples used in morphological analyses have poor quality and thus were not included in the population genetic analyses. Contigs were generated from merged forward, and reverse, sequences and their consensus sequences were aligned using Genious Prime 2020.2. Alignment was refined manually to generate a complete alignment of 840 base pairs. We reconstructed a time-calibrated phylogenetic tree using BEAST v1.10.4 . We incorporated seven species (nine sequences in total) from GenBank as an outgroup (Appendix 1b). Species included in the outgroup are the closest relatives of A. lanyuensis according to the phylogenetic tree inferences of Su and Smith (2014). The program jModelTest2 v. 2.1.10 was used to calculate the best-fit nucleotide substitution model for the CO1 gene using the Akaike Information Criterion (AIC) (Posada, 2008).
Additionally, nucleotide and haplotype diversity of the in-group sequences were calculated based on the PAIC boundaries and current island boundaries using DnaSp v.6.12.03 (Rozas et al., 2017).
Haplotype networks were also created in TCS v.1.21 (Clement et al., 2000) and displayed as a final network using tcsBU v.1.0 (Múrias dos Santos et al., 2016). We conducted an isolation by distance (IBD) test among PAIC islands through Mantel's test of correlation between Edward's distances and Euclidian geographic distances.
IBD test was implemented in R package "adegenet" using the mantel.randtest function (Jombart, 2008). Cline and distant patches of points were checked using the 2-dimensional kernel density estimation (kde2d) in R package "MASS." Gene flow among current islands was further assessed by calculating pairwise Fixation indices (F ST ) using the R package "StAMPP" (Pembleton & Pembleton, 2013).

| Biogeographical analyses
The ancestral geographic ranges were reconstructed by two programs: R package "BioGeoBEARS" (Matzke, 2014), and Reconstruct Ancestral State in Phylogenies (RASP) (Yu et al., 2015). The best-fit historical biogeographical model selection was conducted among six available models in "BioGeoBEARS": DEC, DEC+j, DIVALIKE, DIVALIKE+j, BAYAREALIKE+j (Matzke, 2014). We applied the bestfit historical model (BAYAREALIKE+j) with the highest corrected Akaike information criterion (AICc) weights to the time-calibrated BEAST trees dataset and consensus tree dataset. Additionally, we applied the Bayesian Binary MCMC (BBM) and Statistical Dispersal-Vicariance models in RASP as alternative biogeographical reconstruction analyses.
We designated the geographical distributions of A. lanyuensis according to PAIC islands, while the known geographical distribution of the outgroup was based on the descriptions from World Spider Catalog (2020) and other published literature. There were five current distinct geographical areas included for the in-group: Orchid and Mindanao PAIC (E). Five geographical areas were also included for the outgroup, namely Sundaland (F), Papua New Guinea (G), Japan (H), China (I), and Australia (J). The geographical range allowed at each node was set up to four geographical areas since no extant species occupied more than four geographical areas. Additionally, the density of evolutionary events such as dispersal and vicariance events was calculated and visualized along a time-calibrated tree.

| RE SULTS
Argyrodes lanyuensis samples were collected from 13 sampling sites distributed across Orchid Island, Taiwan, the main (northern) component of Luzon Island, its southern Bicol Peninsula, Palawan, Negros, and Panay islands; plus, the northern, eastern, and southwestern Alternatively, we used the PLS-DA, which emphasized a dimension reduction technique for handling multicollinearity data (Vignon, 2011), to detect the morphological clustering among samples. Because individual samples were assigned according to PAICs a priori, the PLS-DA score plot was able to discriminate PAIC clusters in both males and females (Figure 2c Table 1a). The overall discrimination method based on PLS-DA among the male samples was found to be significant (CER = 0.512, p-value < .05; Table 1a). Thus, we inferred two morphologically discrete clusters for male data as Mindanao+Palawan+West Visayan populations and Luzon+Orchid Island populations were undifferentiated ( Figure 2c). In contrast, we did not find obvious differentiation in PLS-DA plot with the female data; however, the Orchid Island and Mindanao PAIC clusters were significantly different from each other (Orchid Island vs. Mindanao PAIC: CER = 0.277; p-value < .05; Table 1b). Nonetheless, the overall discrimination method based on PLS-DA among the female samples was found to be nonsignificant (CER = 0.567, p-value > .05; Table 1b), which is consistent with our initial PCA results.
To visualize and explore the correlations among variables, we used the latent components in the PLS-DA to display a loading vector plot. The loading vector plot demonstrates the importance of each variable and its contribution to the overall variance in males and females. Figure 3 shows the results of the male and female load- The same lack of pattern is also apparent in our BEAST maximum clade credibility (MCC) tree (Figure 4b), which shows two, strongly supported (Posterior Probability, or PP = 1.00) major clades, each of which exhibits no differentiation among PAIC or current island boundaries. Furthermore, all nodes within these two major clades have low posterior probability support (PP < 0.5), which is surprising given that CO1 is a rapidly evolving mitochondrial gene region. The  Figure 5b shows the probability density of evolutionary events along the time-calibrated tree. We observed a consistent higher probability density of dispersal events than vicariance events that started from node 198, specifically at ~3 MYA (Miocene-Pliocene epochs) when A. lanyuensis diverged from the outgroup. Dispersal events continued toward later nodes wherein more dispersal events have occurred (Figure 5b). Therefore, based on our phylogenetic analyses and biogeographical reconstruction analyses, we reject the strict PAIC biogeographical patterns/ predictions and the recent southward colonization (north-to-south prediction) and thus accept the south-to-north colonization as our best interpretation, but with little to no differentiation due to recent dispersal events and in response to a wide array of host species during range expansion.

| D ISCUSS I ON
Our study demonstrated an updated geographic distribution of Argyrodes lanyuensis that covers almost the entire Philippine Archipelago, aside from Orchid Island, Taiwan, on which this species previously was thought to be endemic (Figure 1a). This species exhibits two phenotypically differentiated units in male morphology (Orchid Island Taiwan+Luzon, Philippines populations vs. Palawan+West Visayan+Mindanao populations; Figure 2c). Our can be seen through their strong habitat affinities (Gillespie, 2016).
For example, genetic structure was observed in excellent dispersalist, Nephila pilipes (Kuntner & Agnarsson, 2011;Su et al., 2007), and Argiope bruennichi (Krehenwinkel et al., 2016). However, we could not observe local adaptations in the case of A. lanyuensis. The specific behavioral phenotype of this species, which is a "generalist" kleptoparasite, could explain the limited differentiation exhibited in this species and implies higher tolerance on different host webs (a case of ecological adaptation) without specialized functions in host-specific feeding strategies. Other spider kleptoparasites (e.g.,

A. fissifrons and A. miniaceus) utilize webs of specific host spiders to
forage prey items (Tso & Severinghaus, 2000) and are thus called ecological "specialist" kleptoparasites (Su et al., 2018). These specialists demonstrate a strong association of these kleptoparasites to their specific host species which in turn may have caused geneticstructured populations across different islands in the Australasian region (Su & Smith, 2014). Thus, we assume that specialized kleptoparasitism could interrupt gene flow between different groups or  . Each subpopulation of P. bifoliate appears to specialize on a species-specific host plant, per PAIC island . In contrast, the results presented here could be a special case for the Philippines archipelago in that we estimate a deeper, pre-Pleistocene temporal divergence time, and yet we did not detect any clear differentiation among PAICs or modern, current-day islands.
The phenotypic clustering evident in males from Orchid Island Orchid Island, respectively. The phenotypic variations observed in these two populations, specifically in the leg I and palp bulb, could be attributed to sexual selection in males. Male A. lanyuensis typically have longer Leg I than females (Yoshida et al., 1998), in which similar observations were recorded in this study (Appendix 2). Leg 1 was usually used by both male and female A. lanyuensis, for moving around the web to locate the host's silk and prey items for food consumption (Yoshida et al., 1998). For most of the Argyrodes spiders, leg I is very important in the male-male competition for female copulation, which results from highly modified intrasexual selection in males (Whitehouse, 1991(Whitehouse, , 2016. With these phenotypic variations in the leg I and palp bulb lengths between the ancestral (Mindanao) and recent populations (Orchid Island), we hypothesize that different mating strategies may evolve in recent populations given the selective pressures in the new environment. In argyrodinae spiders, species-specific differences and intersexual differences in foraging strategies have been noted (Cangialosi, 1990;Kerr, 2005;Tso & Severinghaus, 2000). Female A. lanyuensis may be able to utilize the same foraging strategies across different spider hosts in which the functional genes for a specialized foraging behavior are not well expressed, even though a unique form of foraging strategy (host silk consumption) has been noted on this species in Orchid Island (Yoshida et al., 1998). Hence, our results on the population structure of females could be related to their foraging behavior. On the other hand, mating strategies in males could lead to the morphological differentiation of this spider kleptoparasite generalist (Whitehouse, 2016). However, these results should be further validated due to the limited sample size and genetic markers.
Our study demonstrates the possible exchange of taxa between two geographical entities. In this case, faunal transfers (dispersal) were possible between Taiwan and the Philippines through the Luzon-Taiwan strait, in which dispersal events originated from the Philippines. This study added to the cases of Philippine fauna that have been recorded to disperse from the Philippines to Orchid Island. These include Eutropis cumingi (skink; Ota & Huang, 2000), Polypedates leucomystax (frog; Kuraishi et al., 2009;Ota, 2004); five species of geckos (Ota, 1987;Siler et al., 2014;Wang, 1962), two species of butterflies, Macroglossum ungues cheni (Yen et al., 2003), and Catopyrops ancyra almora (Lu & Hsu, 2002). In contrast, other taxa (e.g., plants, shrews) have been recorded to disperse in the Philippines from Orchid Island (south-to-north colonization; Dickerson, 1928;Esselstyn & Oliveros, 2010;Oliveros & Moyle, 2010). Thus, careful analyses should be done for the diversification of taxa along the Philippine-Lanyu oceanic island chain.

| CON CLUS ION
In conclusion, our results revealed the presence and widespread distribution of A. lanyuensis in the Philippines, far beyond its originally assumed microendemic distribution in Orchid Island, Taiwan.
Our study also emphasized northward colonization of A. lanyuensis from the Philippines toward adjacent Orchid Island, Taiwan, through recent dispersal events. The molecular data highlight the importance of behavioral phenotype such as foraging behavior, rather than isolation by distance, sea-level vicariance, and climatic oscillations (e.g., PAICs Paradigm biogeographic isolation) as drivers of diversification of kleptoparasitic spiders. However, it is also important to note the structured variations we observed in males for the northern populations, which possibly attributed to mating strategies. Future work on this study system may be best informed by a higher coverage of genomic data, to get a more robust and finely differentiated characterization of its population structure across the Philippines, northwards directing to Orchid Island.

ACK N OWLED G M ENTS
We thank the members of the Ecology and Evolutionary Genomics

CO N FLI C T O F I NTE R E S T S
The authors declare no competing interests.

A PPE N D I X 6
PLS-DA score plots of A. lanyuensis males (a-b) and females (c-d) based on 28 and 265 morphometrics, respectively. Individuals are plotted against components 1 and 2, grouped according to current island boundary (a and c) and locality (b and d)