A Natural Colonisation of Asia: Phylogenomic and Biogeographic History of Coin Spiders (Araneae: Nephilidae: Herennia)

Abstract

Reconstructing biogeographic history is challenging when dispersal biology of studied species is poorly understood, and they have undergone a complex geological past. Here, we reconstruct the origin and subsequent dispersal of coin spiders (Nephilidae: Herennia Thorell), a clade of 14 species inhabiting tropical Asia and Australasia. Specifically, we test whether the all-Asian range of Herennia multipuncta is natural vs. anthropogenic. We combine Anchored Hybrid Enrichment phylogenomic and classical marker phylogenetic data to infer species and population phylogenies. Our biogeographical analyses follow two alternative dispersal models: ballooning vs. walking. Following these assumptions and considering measured distances between geographical areas through temporal intervals, these models infer ancestral areas based on varying dispersal probabilities through geological time. We recover a wide ancestral range of Herennia including Australia, mainland SE Asia and the Philippines. Both models agree that H. multipuncta internal splits are generally too old to be influenced by humans, thereby implying its natural colonisation of Asia, but suggest quite different colonisation routes of H. multipuncta populations. The results of the ballooning model are more parsimonious as they invoke fewer chance dispersals over large distances. We speculate that coin spiders’ ancestor may have lost the ability to balloon, but that H. multipuncta regained it, thereby colonising and maintaining larger areas.

1. Introduction

Biogeography is a scientific field that integrates evolutionary hypotheses, contemporary and fossil taxonomic distributions and time calibrated phylogenies. Nevertheless, modern biogeography has struggled to become an exact science for several reasons. First, time calibrated phylogenies often yield unreliable topologies and/or divergence times [1] or produce very wide margins of error [2]. Consequently, time estimates of divergence and speciation events remain vague, and hypothesis testing imprecise (but, see Magalhaes et al. [3]). Second, organism-specific biology is usually not accounted for within historical biogeographic reconstructions. When it is, the organism-specific resolution rarely goes beyond basic trait-state binning [4], e.g., winged versus pedestrian versus aquatic animals, etc. Third, although the Earth’s tectonic and climatic histories represent essential variables for the distribution of organisms, their precise reconstructions through insets of time have not been integrated into biogeographic algorithms.

We have recently discussed this gap in methodology in a biogeographic study of a globally distributed spider family [5]. We suggested and demonstrated a novel method of fine-tuning biogeographical analyses by combining a robust phylogeny and specific organismal biology with dispersal probability estimates, based on concrete measurements between geographical regions in the geological past. This approach proved suitable for analyses across large geographical areas, where geological reconstructions are sufficiently accurate and for organisms whose dispersal biology is well understood. When the geological past of the area of interest is not as clear and species biology is unknown, however, this methodology has to be modified. Here, we produce a comprehensive phylogeny of coin spiders (Nephilidae: Herennia) using phylogenomic data and test and discuss an alternative approach to biogeographical inference to the one proposed previously.

Among the golden orbweaver spiders of the family Nephilidae (catalogued as Nephilinae-Araneidae; we here follow the family classification proposed by [6]), coin spiders (genus Herennia Thorell, 1877 [7]) are the most species-rich genus with over a dozen species distributed in tropical Asia and Australasia. As all nephilids, they exhibit extreme sexual size dimorphism with males but a fraction of female size [6]. Unlike the remainder of genera though, all Herennia build arboricolous (“tree-hugging”) ladder webs on tree trunks [6,8]. A newly updated Herennia taxonomy (in preparation) recognises 14 species of coin spiders (three of which have not yet been formally described and thus feature manuscript names in quotation marks). With the exception of H. multipuncta Doleschall 1859 [9], they are distributed narrowly, mostly as island endemics. Herennia gagamba Kuntner 2005 [8] and H. tone Kuntner 2005 are found in the Philippines, H. “tsoi” Kuntner et al., (in preparation) in Taiwan, H. “maj” Kuntner (in preparation) in Vietnam, H. etruscilla Kuntner 2005 in Java, H. “eva” Kuntner (in preparation) in Sulawesi, H. deelemanae Kuntner 2005 in Borneo, H. jernej Kuntner 2005 in Sumatra, H. sonja Kuntner 2005 in Borneo and Sulawesi, H. papuana Thorell 1881 [10] in New Guinea and Australia, H. agnarssoni Kuntner 2005 in the Solomon Islands, H. milleri Kuntner 2005 in New Guinea and New Britain and H. oz Kuntner 2005 in northern Australia (Figure 1). In contrast, H. multipuncta is distributed throughout southern India, Indochina and the Philippine and Indonesian archipelagos (Figure 1). Unlike other species, which are obligatory arboricoles in pristine forests, H. multipuncta is synanthropic, frequently found in managed habitats, and lives in sympatry with other, narrower endemic species of the genus [8]. This fact has sparked speculation on the invasive origin of H. multipuncta super-range [8], but this hypothesis has remained untested.

Although coin spiders exhibit several intriguing biological features, many aspects of their biology remain unexplored. A prior revision of the genus discussed the taxonomy, biology and biogeography of the 11 then-known species [8]. It suggested purely Australasian speciation of coin spiders and proposed the “Herennia line”, west of Wallace’s and Huxley’s lines, which was only crossed by one then-known species, H. multipuncta. Since then, several new species have been recognised, some of them inhabiting areas west of the proposed line (Kuntner et al., in preparation). A recent nephilid biogeographic study [5] inferred historical biogeography of 10 species with available genetic data; however, that study was global and, thus, its biogeographical resolution was necessarily insufficient to resolve the Herennia biogeographic history. Nonetheless, Turk et al. [5] suggested an Indomalayan origin of the genus with recent colonisation of Australasia by the ancestor of H. milleri and H. multipuncta.

Here, we aimed to answer the following three main questions: (i) what is the sequence and chronology of coin spider dispersal from their origin to the present distribution, (ii) how would alternative dispersal biologies influence this pattern and (iii) is the unusually large range of H. multipuncta a result of human activity, and is the species invasive? We used phylogenomic data to construct a species-level phylogenetic scaffold, and then used classical phylogenetic markers to infer the most comprehensive population-level phylogeny and chronogram of coin spiders to date. We used this reference phylogeny to infer coin spiders’ historical biogeography by adapting the methodology proposed by Turk et al. [5]. We tested two alternative models, each assuming a different type of dispersal, while accounting for the complex geological past of Australasia (Figure 2). The first model (A) assumes active dispersal via ballooning [11]. Ballooning behaviour has been observed in other nephilid species [12], but not yet in coin spiders. This type of dispersal promotes island colonisation, but also facilitates gene flow maintenance across large distances, thus inhibiting island endemism. The second model (B) assumes short distance random walking dispersal during the search for vacant habitat as the main method of dispersal. It allows for passive dispersal over long distances of connected lands given enough time. Neither model completely excludes rare chance occurrences of long-distance dispersal with wind currents.

To address our main questions highlighted above, we inferred the most comprehensive species- and population-level phylogenies of coin spiders to date and calculated the most plausible biogeographic histories and dispersal trajectories using the two alternative models. If our data revealed H. multipuncta intraspecific splits are too recent to be resolved via phylogenetic time calibration, this would speak for a human-induced dispersal of H. multipuncta into other Herennia species ranges. In contrast, if intraspecific splits were resolved, and reconstructed as more ancient than human presence in this area (over 50,000 years [13]), then a natural colonisation would be implied. Generally, the estimated ages of nodes within H. multipuncta did not support the human-driven dispersal hypothesis. We speculate that ballooning ability and/or propensity was lost in coin spiders (except for H. multipuncta), resulting in a high occurrence of island endemism in the genus.

2. Materials and Methods

2.1. Species-Level Phylogeny

Currently, tissue material and genetic data are available for ten of the 14 known Herennia species, while the remaining four (H. agnarssoni, H. deelemanae, H. jernej and H. sonja) are only known from holotypes and could not be included here. Prior to our study, the best supported nephilid phylogeny used Anchored Hybrid Enrichment (AHE) data to resolve phylogenetic relationships among 22 species including five species of Herennia [6]. We here expand this taxon coverage to include eight Herennia species, of which H. multipuncta was represented by a specimen from Sri Lanka and one from Laos in order to test this widespread species monophyly (Supplementary Spreadsheet S1).

We employed the AHE targeted-sequencing approach for spiders (outlined in Hamilton et al. [14]) to target 585 single copy orthologous loci from across the genome. These loci have been shown to possess sufficient variation for resolving both shallow and deep-scale evolutionary relationships throughout the Araneae, e.g., [14,15]. These data have also been used to recover inter- and intrageneric relationships, as well as inter- and intraspecies relationships within a range of spider families [6,16,17,18,19,20,21]. Library preparation, enrichment, sequencing, assembly, alignment and phylogeny construction from AHE data followed the procedures described in Kuntner et al. [6]. The full AHE matrix is available as Supplementary File S1.

For two species, H. papuana and H. tone, we were unable to recover AHE data. To identify their phylogenetic placement, we ran maximum likelihood (RAxML) and Bayesian (MrBayes) phylogenetic inference analyses only on COI sequences, obtained for the ten species. The primers and PCR amplification protocols to obtain partial COI sequences followed Kuntner et al. [22] (matrix available as Supplementary File S2; COI GenBank accession codes OK017092 to OK017142). The placements of Herennia “maj” and H. milleri, the latter likely due to low COI sequence quality, were poorly supported and inconsistent with the AHE topology. To eliminate the influence of these “rogue taxa” on H. papuana (topological placement of H. tone was consistent regardless of their presence), these two terminals were removed from the COI only dataset and the analyses were repeated. Results produced by MrBayes were inspected with Tracer to ensure effective sample sizes were >200.

Divergence dating calibration was performed in BEAST2 [23], again using only COI sequences. We used a relaxed log normal clock and set the bModelTest [24] as the nucleotide substitution model and a birth-death tree prior. The ucldMean prior was set as normally distributed with a mean of 0.0199 and a standard deviation of 0.001 following Bidegaray-Batista and Arnedo [25]. The topology was constrained as described above, whereas clade ages and their confidence intervals were constrained in three nodes that correspond directly to those acquired by Kuntner et al. [6]. BEAST ran on four MCMC chains for 10 million generations. Results were checked with Tracer and summarised with TreeAnnotator with a 10% burnin.

2.2. Population-Level Phylogeny

In the population-level phylogeny, species were represented by a varying number of samples, ranging from one in H. tone, H. milleri and H. papuana to 26 in H. multipuncta. Again, we inferred the phylogeny in BEAST2, using COI and 28S sequences. For 28S sequences, primers and PCR amplification protocols again followed Kuntner et al. [22] (concatenated matrix available as Supplementary File S3; 28S GenBank accession codes OK017174 to OK017212). Where an individual lacked data for both genes, only one was used. We employed a relaxed log normal clock and set the bModelTest. The Coalescent Bayesian Skyline prior, allowing for stochastic changes in population sizes through time, was chosen as the tree prior following Ritchie et al. [26]. COI priors were set as before, while the ucldMean prior for 28S was set as normally distributed with a mean of 0.0011 and standard deviation of 0.0003 after Bidegaray-Batista and Arnedo [25]. The ages of all species-level splits were constrained to those recovered in the previous phylogeny. BEAST ran on four MCMC chains for 70 million generations to ensure large enough effective sample sizes. Again, a consensus tree was obtained with TreeAnnotator with a 10% burnin.

2.3. Inference of Biogeographic History

For biogeographic analyses, we pruned the reference population-level phylogeny so that each population (locality) was represented by a single specimen (hereafter referred to as “pruned phylogeny”). The individual with the most complete sequence was chosen as the representative of each population. This narrowed the population tree to 29 tips, with the number of representatives per species ranging from one in H. oz, H. milleri, H. papuana and H. tone to 13 in H. multipuncta. We treated the known areas of species occupancies within 10 biogeographic regions. These consisted of three continental landmasses—mainland South-East Asia, Australia and India—and seven islands or archipelagos—Sulawesi, Sumatra, New Guinea, Sri Lanka, Taiwan, the Philippines and the island pair Java and Lombok, for simplicity referred to as Lesser Sunda islands. As the tips in the phylogenetic tree are individual samples, necessarily only inhabiting a single biogeographic region, each tip was attributed one region.

Following the rationale developed in Turk et al. [5], dispersal probabilities were fine-tuned to reflect the varying geographical configuration of biogeographic regions during the area’s lively geological past. Physical distance among landmasses was used as a proxy for dispersal probabilities, scored separately in six time slices, each spanning 5 million years. As the geological history of the area, especially Indonesian islands, is extremely complex and thus difficult to reconstruct with precision, we binned dispersal probabilities into five categories. Following the argumentation in Turk et al. [5], we attributed a 95% dispersal probability to pairs of geographic regions in physical contact, where dispersal is likely but not necessary, and a 5% dispersal probability where distances exceed 4000 kilometres. For distances between regions of 1000 kilometres or less, a 75% dispersal probability was assigned, for distances between 1000 and 2500 kilometres, a 50% dispersal probability, and for distances between 2500 and 4000 kilometres, a 25% dispersal probability (Figure 2a). If a region had not yet emerged or was sunk during a time slice, it was disallowed in the “areas allowed” option in RASP (see below).

In contrast to the previous paper, however, we attributed these probabilities in two ways, differing in the definition of physical contact among regions (Figure 2b). In time slices where, for example, three regions were consecutively connected by land, model A attributed a 95% dispersal probability only between the middle and marginal areas, but not between the marginal areas themselves, thus accounting for the larger physical distance between them, despite physical connection via the middle area. This model, which we term “ballooning dispersal”, puts emphasis on long-distance dispersal via ballooning as the main method of dispersal in coin spiders. Model B, in contrast, attributed a 95% dispersal probability for all three pairs of regions in the previous example. Biologically, model B assumes coin spiders largely disperse over land with small-step, gradual expansion; we term this “walking dispersal”. Given enough time, this model assumes that spiders can reach all physically connected areas equally likely.

Geological data (Supplementary Spreadsheet S1) were compiled from a tectonic reconstruction model [27] via GPlates plate tectonics visualisation software [28] and geological literature [29,30,31,32,33]. We reconstructed the genus’ historical biogeography with RASP 4.0 beta [34], comparing all six included biogeographical models. The maximum number of ancestral areas occupied was set to three. We evaluated model fit through weighted AIC_c values (AIC_c_wt), expressing the model’s relative probability, corrected for small sample sizes.

3. Results

3.1. Phylogenies

The AHE phylogeny placed eight Herennia species unequivocally and with overwhelming support in all nodes (Figure 3a). The oldest split from the MRCA of all coin spiders was recovered in the Vietnamese H. “maj”. Herennia “eva” and H. gagamba then branched off as sister species, followed by H. “tsoi” and H. multipuncta, H. milleri and finally H. etruscilla and H. oz as sister species. As for the COI sequences, both RAxML and MrBayes consistently placed H. tone as sister to H. “eva” with a 99% support in MrBayes (Figure 3b). After the elimination of “rogue” taxa that interfered with topological stability, H. papuana was placed as sister to ((H. “tsoi”, H. multipuncta), (H. milleri, H. etruscilla, H. oz)) consistently using both methods, again with a high node support of 92% in MrBayes (Figure 3b).

In the population-level phylogeny (Supplementary Figure S1), samples of the same species always grouped together; however, samples from the same locality often did not (e.g., H. etruscilla populations from Java, and H. multipuncta populations from Laos, Vietnam, Malaysia, Yunnan and Hainan). In the pruned phylogeny containing one tip per population, used in biogeographical reconstruction (see Figure 4), the divergence dating revealed frequent within-species cladogenesis during the last few million years.

3.2. Biogeographical Reconstruction: Model A

RASP identified DIVALIKE+j as the best model for the data (Table 1, Figure 4). The node uniting all Herennia taxa received strong support for a wide ancestral distribution in Australia, mainland SE Asia and the Philippines (61%). Although H. “maj” persisted only in mainland SE Asia, the most recent common ancestor (MRCA) of all other species persisted in the Philippines and Australia (65%). The clade containing H. gagamba, H. tone and H. “eva” remained in the Philippines, with the latter species colonising Sulawesi sometime during the last four million years.

The MRCA of all other species remained in Australia only, with H. papuana persisting in the same area until the present. The MRCA of H. “tsoi” and H. multipuncta shifted from Australia to mainland SE Asia, from where H. “tsoi” colonised Taiwan, while H. multipuncta remained in SE Asia with consecutive colonisations of Sri Lanka, India, Sulawesi, Java and Sumatra during the last four million years. From the remaining Australian MRCA, H. milleri colonised New Guinea, while its sister Australian clade diversified into the Australian H. oz and H. etruscilla, which shifted to the Lesser Sunda islands in the last few million years.

3.3. Biogeographical Reconstruction: Model B

Again, DIVALIKE+j received the highest statistical support in RASP (Table 1, Figure 4). All nodes except the MRCA of H. “tsoi” and H. multipuncta and their intra-species nodes received nearly identical support as in model A. Here, the reconstruction for the MRCA of H. “tsoi” and H. multipuncta was ambiguous, with similar support for Sulawesi (35%), mainland SE Asia (23%) and India (21%). The basal H. multipuncta was attributed the same three regions in almost identical shares (36, 23 and 20%, respectively). The species was reconstructed to have colonised Sri Lanka and India from Sulawesi. Next, it colonised (and remained in) the Lesser Sunda islands (49%) or mainland SE Asia (51%), and then dispersed across the SE Asian mainland and to Sumatra.

Comparing the fit of DIVALIKE+j between the two models through the calculation of the Bayes factor from the log of the likelihood scores (B = 0.98; after Kass and Raftery [36]) revealed a non-significant difference.

3.4. Human-Induced Dispersal of H. multipuncta?

Neither of the two models supported the hypothesis of a recent, human-induced dispersal of H. multipuncta. Node ages, recovered for phylogenetic splits among the included populations in the pruned phylogeny, ranged from 3.72 million years ago (mya) in the MRCA of all included populations to 0.13 mya (132,200 years) in the split between the Yunnanese and Cambodian populations. In the latter, the confidence interval ranged from 516,800 to 800 years ago, a time frame potentially compatible with human-induced dispersal, albeit between two areas in relatively close proximity. The only other node whose estimated age fit within the timeframe of human presence in the area was the split between the Vietnamese and Lao populations. It was dated to 277,300 years ago, with a confidence interval ranging from 757,500 to 26,900 years ago. The two areas are adjacent with no obvious barriers between them, implying that natural dispersal is very likely.

4. Discussion

In this study, we provide a test case of what next generation biogeographic inference should optimally encompass: a robust phylogenetic/phylogenomic framework for a comprehensive population-level ancestral area reconstruction that at the same time accounts for geological dynamics and species biology. To this end, we used a modified version of our previously proposed methodology [5], comparing two methods of dispersal probability quantification. Both models suggested a wide ancestral range and relatively old splits (from 3.72 to 0.13 mya) between terminals of H. multipuncta, strongly implicating a natural, not anthropogenic colonisation of the areas that constitute an extremely wide contemporary range of this species.

4.1. Phylogenetic Placements

The species-level phylogeny recovered unexpected relationships among species with overwhelming support (Figure 3a,b). Within the available taxon sample, Herennia “maj” was the first species to split from the coin spider MRCA, a placement that was not recovered in a prior COI-based phylogeny in Turk et al. [5]. As surprising as this may be, we believe the relationship is not artefactual, given our understanding of the robustness of AHE phylogenomic topologies in nephilids (H. “maj” was not included in Kuntner et al. [6]). In the population-level phylogeny, species were represented by a varying number of specimens. In H. multipuncta, the largest sample, specimens from the same population failed to group together (Supplementary Figure S1). Perhaps genetic differences between them are not (yet) large enough, as they maintain modest gene flow among populations. This matches our emerging hypothesis of heightened dispersal propensity in this species, relative to others, as proposed below.

4.2. Biogeographic Inference with Two Models

Unlike our previous study [5], we only have limited knowledge about the dispersal biology of coin spiders. As in all motile animals, they can be assumed to, at least, undergo gradual, slow dispersal through suitable habitat during stochastic search for vacant space. On the other hand, we might expect coin spiders to practice active ballooning, as observed in the related Nephila pilipes [12], but no field or experimental evidence that would support this assertion currently exists. Even in the absence of ballooning, however, chance dispersal across geographic barriers, such as the sea, have to occur, otherwise they would not be able to simultaneously inhabit landmasses that have not been connected during the relevant timeframe, such as mainland SE Asia and Australia. The two proposed types of dispersal have different consequences: ballooning would make it easier to colonise new areas and maintain gene flow across the entire range, but the evolution of island endemism is less likely. Conversely, relying on “walking” dispersal with occasional, chance long-distance dispersal would facilitate genetic isolation and thus speciation and island endemism, but would restrict gene flow maintenance across large ranges. Regardless, we do not assume that all coin spider species necessarily exhibit identical dispersal behaviours and propensities.

Thus, in the absence of a better understood, and experimentally tested, organism-specific dispersal biology, we resorted to engaging in two biogeographical approaches (Figure 2). Model A assumed less likely dispersal between areas further apart, even if connected by land, because ballooning was taken as the default method of dispersal. In contrast, model B assigned equal dispersal probabilities to all physically connected areas, as it assumed spiders spread “on foot”, gradually, over millions of years. Although such passive, “walking” dispersal is much slower, range expansion of, for example, only 10 m per year (which we deem as a distance easily overcome by orbweaving spiders) would theoretically allow for a spread across 10,000 kilometres over the course of a single million years, which is approximately the extent of the entire extant genus range. While the comparison of the two best-fitting models using the Bayes factor did not show significant differences between them, they nonetheless provided somewhat discrepant results.

Deep nodes near the root of the phylogeny were reconstructed with nearly identical probability proportions for ancestral areas by both models (Figure 4). A wide ancestral distribution over mainland SE Asia, Australia and the Philippines fragmented, with the ancestral mainland population leading to the Vietnamese H. “maj”, and the Australian and Philippine populations giving rise to all other species. Interestingly, mainland SE Asia was not re-colonised until 24–29 million years later, depending on the model. In both models, H. tone and H. gagamba maintain the ancestral range in the Philippines, while H. “eva” disperses to Sulawesi sometime during the last four million years. This inferred colonisation is plausible assuming either model of dispersal, requiring an active or chance dispersal across approximately 600 km of sea [32]. Either type of dispersal on this route may have been further facilitated by a possible island chain connecting Sulawesi to the Philippines in the Neogene [33].

The most likely ancestral areas of the (H. “tsoi” + H. multipuncta + H. milleri + H. oz + H. etruscilla) and (H. milleri + H. oz + H. etruscilla) MRCAs are reconstructed to Australia by both models (Figure 4). The latter clade features another recent over-water colonisation of H. etruscilla from Australia to Lesser Sunda islands. A salient difference between the two models is the inferred dispersal of H. multipuncta (Figure 4). In model A, the MRCA of H. multipuncta and H. “tsoi” is already distributed in mainland SE Asia, from where H. “tsoi” colonises Taiwan during the last 6 million years. This does not exceed the understood age of the island (approximately 9 million years [29]), which has remained close to, and even connected with, the Asian mainland by a land bridge during Pleistocene glaciations [30]. In H. multipuncta, model A suggests gradual expansion from mainland SE Asia to Sri Lanka, India, across Indonesian islands and finally within mainland SE Asia itself. On the other hand, model B inconclusively places the MRCA of H. multipuncta and H. “tsoi” to Sulawesi, a challenging proposition. At 6 mya, the distance between Sulawesi and Taiwan, which H. “tsoi” supposedly crossed, was not considerably shorter than today (approximately 2300 kilometres), having been separated by the Philippine archipelago.

Deeper nodes within H. multipuncta are also ambiguous in model B, with similar probabilities for several ancestral areas, but it generally suggests gradual dispersal of the species from Sulawesi to Sri Lanka and India, colonisation of Java directly or via mainland SE Asia, and finally a spread within SE Asia and colonisation of Sumatra. Dispersal from Sulawesi to Sri Lanka during this time would require crossing a >4000-kilometre distance either by chance dispersal, which is rarely observed on such a scale in spiders in general (as discussed in Turk et al. [5]), or over land. As opposed to islands such as Sri Lanka, New Guinea and Taiwan, however, it is unclear whether Sulawesi formed land connections with the mainland or mainland-connected islands, such as Borneo, during the time of Plio-Pleistocene sea level changes [33]. Furthermore, a subsequent colonisation of India, requiring another >4000-kilometre dispersal event shortly after the first, seems highly unlikely. On the other hand, all sampled H. multipuncta populations except those from Sulawesi, India and Sri Lanka were connected by land during the Pleistocene on the landmass Sundaland [32]. Dispersing around this landmass could indeed be performed by “walking” dispersal, but, as seen above, even model B infers several oversea dispersals, either by chance or actively.

Curiously, species such as H. “maj” seem to maintain a narrow distribution despite living on the Asian mainland, where they could, across millions of years, passively disperse over a much wider area (assuming our model B with no other limiting factors). We speculate such species are either restricted to specific habitats that are not continuous, decreasing the chance of successful active or passive dispersal, or are confined to their range through ecological competition with other species and genera in adjacent areas.

4.3. A Natural “Coinquest” of H. multipuncta

The pruned population-level phylogeny showed relatively recent splitting between the sampled populations of H. multipuncta (median node ages ranging from 3.72 to 0.13 mya), however, not recent enough to infer human-related dispersal. If that were the case, nucleotide sequence divergencies would have been predictively low, having only accumulated changes over the past few thousand years (contra to what we see in our phylogeny). In the two cases where confidence intervals of node ages do overlap with human presence in the species’ range, the pairs of geographical areas are either adjacent or in close proximity, making the exclusion of natural colonisation difficult.

Why, then, is a super-wide present-day range of H. multipuncta unique among all Herennia? We hypothesise that although nephilid ancestors perform active ballooning, retained in Nephila and Trichonephila [9,34,35], coin spiders secondarily lost the ability to balloon. Such a loss of dispersal ability is a common phenomenon in island spider biology, severely limiting gene flow maintenance and often leading to single-island endemism [37,38]. Herennia multipuncta might have regained this ability, allowing it to disperse across suitable habitat and inhabit most of the genus range, sometimes overlapping with more ancient Herennia spp ranges. This would also allow it to sustain some degree of gene flow among populations, which would, in turn, explain the recovered phylogenetic picture of non-monophyly of sympatric specimens of certain populations (from Vietnam, Laos, Malaysia, Hainan and Yunnan). In many parts of its range, H. multipuncta is sympatric with other species (H. etruscilla in Java, H. “eva” in Sulawesi, H. “maj” in Vietnam, H. jernej in Sumatra, etc.) and, according to our dated population phylogeny, this has been the case for millions of years. This pattern suggests that H. multipuncta does not outcompete other sympatric species, perhaps due to subtle differences in ecological niches. In fact, adaptiveness to different habitats might be another trait, specific to H. multipuncta, that enabled its easier and more successful dispersal. These interpretations could be tested in an ecological framework in the near future, as could the ability of H. multipuncta, but not other Herennia, to disperse via ballooning.

4.4. Limitations in Methodology and Future Work

One source of statistical bias in the present study is the incomplete representation of coin spider species in the phylogeny and with it the lack of representation of the missing species’ geographical distributions. The precise phylogenetic placement of the Bornean H. deelemanae, Sumatran H. jernej, Bornean and Sulawesian H. sonja, as well as H. agnarssoni from the Solomon Islands has never been tested outside of a morphological, cladistic framework [8]; therefore, the influence of their distributions on biogeographical reconstruction remains unclear. Furthermore, at the population level, species are not equally represented in terms of specimen number and range cover, also potentially biasing evolutionary relationships and divergence times between populations. That said, specimen collections of coin spiders are scarce. The present study was performed with all genetic material currently available to us.

One of the topics addressed in the study was the dispersal behaviour of coin spiders. Ideally, the presence or absence of active ballooning ought to be tested experimentally. Considering that ballooning is difficult to observe in nature, future research could include subjecting juvenile coin spiders to wind tunnel experiments [12] in a laboratory environment. If performed on multiple species, such an experiment might serve as a test of our hypothesised regained ballooning behaviour in H. multipuncta, but not in other species.

5. Conclusions

We have demonstrated the importance of the understanding of organismal biology in biogeographic reconstructions. In organisms where dispersal is not well understood, testing alternative modes of dispersal through parallel statistical models might prove helpful in uncovering the most likely dispersal biology without direct field observation. By modifying our previously proposed pipeline to account not only for the specifics of the geological history of the area, but also dispersal specifics of the studied organisms, we are further contributing to the development of biogeographic methodology.

Author Contributions

Conceptualisation, E.T. and M.K.; methodology, E.T., K.Č., C.A.H. and M.K.; formal analysis, E.T., K.Č., C.A.H., J.E.B., R.-C.C. and M.K.; investigation, all authors; writing—original draft preparation, E.T. and M.K.; writing—review and editing, all authors; visualisation, K.Č.; funding acquisition, M.K. and S.K.-F. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Slovenian Research Agency, grants P1-0236, P1-0255, J1-9163 and 1000-17-0618.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The new genetic sequences are available on GenBank (accession codes: OK017092 to OK017142 (COI), OK017174 to OK017212(28S)) and also as a supplement to this paper. AHE and geological data are available as a supplement to this paper.

Acknowledgments

We thank A. V. Abramov, I. Agnarsson, G. J. Anderson, M. Bedjanič, D. Court, A. Harmer, M. S. Harvey, S. Huber, P. Jaeger, D. Li, W. Maddison, W. Schawaller and I.-M. Tso for their kind assistance in field work and/or specimen acquisition. We thank I. Magalhaes and two anonymous referees for constructive criticism.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Sampling locations of Herennia individuals used in this study. Encircled numbers denote the number of individuals from the same sampling location. Geographic areas, used for historical biogeography inference (A–J), are colour-coded.

Figure 2. Methods of dispersal probability attribution. Each pink circle represents a geographical area. (a) For non-adjacent pairs of areas, both models attributed dispersal probabilities based on physical distances between the areas. Probabilities were inversely proportionate to the distance and binned into four categories. (b) When areas were in physical contact, model A attributed the maximal, 95% dispersal probability only between directly neighbouring areas, while indirectly connected areas were scored as non-adjacent pairs of areas. In contrast, model B treated all physically connected areas as a single area, attributing a 95% dispersal probability to all pairs of areas.

Figure 3. Species-level phylogenies of coin spiders. (a) AHE-only phylogeny produced with ASTRAL [35], resolving the relationships between eight species of Herennia with Nephilengys as the outgroup as per Kuntner et al. [6]; (b) Species-level phylogeny of the available ten species of Herennia, calculated from COI data. Highlighted are two species lacking AHE data, H. papuana and H. tone. The lock symbols denote age-constrained nodes. Supports for nodes not present in the AHE-only phylogeny (marked with an asterisk) were recovered by MrBayes.

Figure 4. RASP ancestral area reconstruction of two alternative biogeographic models on a COI + 28S population-level phylogeny of Herennia. Model A assumes long-distance dispersal via ballooning as the main method of dispersal in coin spiders, while model B predicts that coin spiders mainly disperse over land with small-step, gradual expansion. Encircled letters signify the likeliest ancestral area in that node, with a combination of several letters indicating an inferred distribution in all those areas. Green area marks the major conflicts in ancestral area reconstruction between the two models.

Table 1

RASP model scores for models A and B. Best supported RASP model is shown in bold. An asterisk (*) in the last column indicates that the model variant allowing jump dispersal (+j) was a significantly better fit for the data than the variant without it.

Model A	LnL	AICc	AICc_wt	Δ −j/+j
DEC	−62.28	129	2.5 × 10−7	*
DEC+j	−47.08	101.1	0.29
DIVALIKE	−55.86	116.2	0.0002	*
DIVALIKE+j	−46.19	99.33	0.7
BAYAREALIKE	−72.01	148.5	1.5 × 10−11	*
BAYAREALIKE+j	−50.27	107.5	0.012
Model B
DEC	−62.86	130.2	3.1 × 10−7	*
DEC+j	−48	103	0.25
DIVALIKE	−56.68	117.8	0.0001	*
DIVALIKE+j	−46.92	100.8	0.74
BAYAREALIKE	−73.4	151.3	8.1 × 10−12	*
BAYAREALIKE+j	−50.89	108.7	0.014

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/d13110515/s1, File S1: AHE matrix (Fasta) with partition file, File S2: COI species-level matrix (Nexus), File S3: concatenated COI and 28S population-level matrix (Nexus), Figure S1: full population-level phylogeny, Spreadsheet S1: geology and geography data.

References

1. Jin, Y.; Brown, R.P. Partition number, rate priors and unreliable divergence times in Bayesian phylogenetic dating. Cladistics; 2018; 34, pp. 568-573. [DOI: https://dx.doi.org/10.1111/cla.12223]

2. Brown, J.W.; Smith, S.A. The Past Sure is Tense: On Interpreting Phylogenetic Divergence Time Estimates. Syst. Biol.; 2018; 67, pp. 340-353. [DOI: https://dx.doi.org/10.1093/sysbio/syx074]

3. Magalhaes, I.L.F.; Santos, A.J.; Ramírez, M.J. Incorporating topological and age uncertainty into event-based biogeography supports paleo-islands in Galapagos and ancient connections among Neotropical dry forests. Diversity; 2021; 13, 418.In press [DOI: https://dx.doi.org/10.3390/d13090418]

4. Crews, S.C.; Esposito, L.A. Towards a synthesis of the Caribbean biogeography of terrestrial arthropods. BMC Evol. Biol.; 2020; 20, pp. 1-27. [DOI: https://dx.doi.org/10.1186/s12862-019-1576-z]

5. Turk, E.; Čandek, K.; Kralj-Fišer, S.; Kuntner, M. Biogeographical history of golden orbweavers: Chronology of a global conquest. J. Biogeogr.; 2020; 47, pp. 1333-1344. [DOI: https://dx.doi.org/10.1111/jbi.13838]

6. Kuntner, M.; Hamilton, C.A.; Cheng, R.-C.; Gregorič, M.; Lupše, N.; Lokovšek, T.; Lemmon, E.M.; Lemmon, A.R.; Agnarsson, I.; Coddington, J.A. et al. Golden orbweavers ignore biological rules: Phylogenomic and comparative analyses unravel a complex evolution of sexual size dimorphism. Syst. Biol.; 2019; 68, pp. 555-572. [DOI: https://dx.doi.org/10.1093/sysbio/syy082] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30517732]

7. Thorell, T. Studi sui Ragni Malesi e Papuani. I. Ragni di Selebes raccolti nel 1874 dal Dott. O. Beccari. Ann. del Mus. Civ. Stor. Nat. Genova; 1877; 10, pp. 341-637.

8. Kuntner, M. A revision of Herennia (Araneae: Nephilidae: Nephilinae), the Australasian “coin spiders”. Invertebr. Syst.; 2005; 19, pp. 391-436. [DOI: https://dx.doi.org/10.1071/IS05024]

9. Doleschall, L. Tweede Bijdrage tot de kennis der Arachniden van den Indischen Archipel. Acta Soc. Sci. Indica-Neerl.; 1859; 5, pp. 1-60.

10. Thorell, T. Studi sui Ragni Malesi e Papuani. III. Ragni dell’Austro Malesia e del Capo York, conservati nel Museo civico di storia naturale di Genova. Ann. Mus. Civ. Stor. Nat. Genova; 1881; 17, pp. 7-27.

11. Bell, J.R.; Bohan, D.A.; Shaw, E.M.; Weyman, G.S. Ballooning dispersal using silk: World fauna, phylogenies, genetics and models. Bull. Entomol. Res.; 2005; 95, pp. 69-114. [DOI: https://dx.doi.org/10.1079/BER2004350] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15877859]

12. Lee, V.M.J.; Kuntner, M.; Li, D. Ballooning behavior in the golden orbweb spider Nephila pilipes (Araneae: Nephilidae). Front. Ecol. Evol.; 2015; 3, pp. 1-5. [DOI: https://dx.doi.org/10.3389/fevo.2015.00002]

13. O’Connell, J.F.; Allen, J.; Williams, M.A.J.; Williams, A.N.; Turney, C.S.M.; Spooner, N.A.; Kamminga, J.; Brown, G.; Cooper, A. When did Homo sapiens first reach Southeast Asia and Sahul?. Proc. Natl. Acad. Sci. USA; 2018; 115, pp. 8482-8490. [DOI: https://dx.doi.org/10.1073/pnas.1808385115] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30082377]

14. Hamilton, C.A.; Lemmon, A.R.; Lemmon, E.M.; Bond, J.E. Expanding Anchored Hybrid Enrichment to resolve both deep and shallow relationships within the spider tree of life. BMC Evol. Biol.; 2016; 16, 212. [DOI: https://dx.doi.org/10.1186/s12862-016-0769-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27733110]

15. Opatova, V.; Hamilton, C.A.; Hedin, M.; De Oca, L.M.; Král, J.; Bond, J.E.; Wiegmann, B. Phylogenetic Systematics and Evolution of the Spider Infraorder Mygalomorphae Using Genomic Scale Data. Syst. Biol.; 2020; 69, pp. 671-707. [DOI: https://dx.doi.org/10.1093/sysbio/syz064] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31841157]

16. Hamilton, C.A.; Hendrixson, B.E.; Bond, J.E. Taxonomic revision of the tarantula genus Aphonopelma Pocock, 1901 (Araneae, Mygalomorphae, Theraphosidae) within the United States. Zookeys; 2016; 2016, pp. 1-340. [DOI: https://dx.doi.org/10.3897/zookeys.560.6264] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27006611]

17. Chamberland, L.; McHugh, A.; Kechejian, S.; Binford, G.J.; Bond, J.E.; Coddington, J.A.; Dolman, G.; Hamilton, C.A.; Harvey, M.S.; Kuntner, M. et al. From Gondwana to GAARlandia: Evolutionary history and biogeography of ogre-faced spiders (Deinopis). J. Biogeogr.; 2018; 45, pp. 2442-2457. [DOI: https://dx.doi.org/10.1111/jbi.13431]

18. Maddison, W.P.; Evans, S.C.; Hamilton, C.A.; Bond, J.E.; Lemmon, A.R.; Lemmon, E.M. A genome-wide phylogeny of jumping spiders (Araneae, Salticidae), using Anchored Hybrid Enrichment. Zookeys; 2017; 695, pp. 89-101. [DOI: https://dx.doi.org/10.3897/zookeys.695.13852]

19. Godwin, R.L.; Opatova, V.; Garrison, N.L.; Hamilton, C.A.; Bond, J.E. Phylogeny of a cosmopolitan family of morphologically conserved trapdoor spiders (Mygalomorphae, Ctenizidae) using Anchored Hybrid Enrichment, with a description of the family, Halonoproctidae Pocock 1901. Mol. Phylogenet. Evol.; 2018; 126, pp. 303-313. [DOI: https://dx.doi.org/10.1016/j.ympev.2018.04.008]

20. Bond, J.E.; Hamilton, C.A.; Godwin, R.L.; Ledford, J.M.; Starrett, J. Phylogeny, evolution, and biogeography of the North American trapdoor spider family Euctenizidae (Araneae: Mygalomorphae) and the discovery of a new “endangered living fossil” along California’s Central Coast. Insect Syst. Divers.; 2020; 4, pp. 1-14. [DOI: https://dx.doi.org/10.1093/isd/ixaa010]

21. Hebets, E.A.; Bern, M.; McGinley, R.H.; Roberts, A.; Kershenbaum, A.; Starrett, J.; Bond, J.E. Sister species diverge in modality-specific courtship signal form and function. Ecol. Evol.; 2021; 11, pp. 852-871. [DOI: https://dx.doi.org/10.1002/ece3.7089] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33520171]

22. Kuntner, M.; Arnedo, M.A.; Trontelj, P.; Lokovšek, T.; Agnarsson, I. A molecular phylogeny of nephilid spiders: Evolutionary history of a model lineage. Mol. Phylogenet. Evol.; 2013; 69, pp. 961-979. [DOI: https://dx.doi.org/10.1016/j.ympev.2013.06.008]

23. Bouckaert, R.R.; Heled, J.; Kühnert, D.; Vaughan, T.; Wu, C.H.; Xie, D.; Suchard, M.A.; Rambaut, A.; Drummond, A.J. BEAST 2: A software platform for Bayesian evolutionary analysis. PLoS Comput. Biol.; 2014; 10, e1003537. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1003537] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24722319]

24. Bouckaert, R.R.; Drummond, A.J. bModelTest: Bayesian phylogenetic site model averaging and model comparison. BMC Evol. Biol.; 2017; 17, 42. [DOI: https://dx.doi.org/10.1186/s12862-017-0890-6]

25. Bidegaray-Batista, L.; Arnedo, M.A. Gone with the plate: The opening of the Western Mediterranean basin drove the diversification of ground-dweller spiders. BMC Evol. Biol.; 2011; 11, 317. [DOI: https://dx.doi.org/10.1186/1471-2148-11-317]

26. Ritchie, A.M.; Lo, N.; Ho, S.Y.W. The impact of the tree prior on molecular dating of data sets containing a mixture of inter- and intraspecies sampling. Syst. Biol.; 2017; 66, pp. 413-425. [DOI: https://dx.doi.org/10.1093/sysbio/syw095] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27798404]

27. Müller, R.D.; Zahirovic, S.; Williams, S.E.; Cannon, J.; Seton, M.; Bower, D.J.; Tetley, M.; Heine, C.; Le Breton, E.; Liu, S. et al. A global plate model including lithospheric deformation along major rifts and orogens since the Triassic. Tectonics; 2019; 38, pp. 1884-1907. [DOI: https://dx.doi.org/10.1029/2018TC005462]

28. Müller, R.D.; Cannon, J.; Qin, X.; Watson, R.J.; Gurnis, M.; Williams, S.; Pfaffelmoser, T.; Seton, M.; Russell, S.H.J.; Zahirovic, S. GPlates: Building a virtual Earth through deep time. Geochem. Geophys. Geosys.; 2018; 19, pp. 2243-2261. [DOI: https://dx.doi.org/10.1029/2018GC007584]

29. Sibuet, J.C.; Hsu, S.K. How was Taiwan created?. Tectonophysics; 2004; 379, pp. 159-181. [DOI: https://dx.doi.org/10.1016/j.tecto.2003.10.022]

30. Voris, H.K. Maps of Pleistocene sea levels in Southeast Asia: Shorelines, river systems and time durations. J. Biogeogr.; 2000; 27, pp. 1153-1167. [DOI: https://dx.doi.org/10.1046/j.1365-2699.2000.00489.x]

31. Park, Y.; Maffre, P.; Godderis, Y.; MacDonald, F.A.; Anttila, E.S.C.; Swanson-Hysell, N.L. Emergence of the Southeast Asian islands as a driver for Neogene cooling. Proc. Natl. Acad. Sci. USA; 2020; 117, pp. 25319-25326. [DOI: https://dx.doi.org/10.1073/pnas.2011033117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32973090]

32. Hall, R. The palaeogeography of Sundaland and Wallacea since the Late Jurassic. J. Limnol.; 2013; 72, pp. 1-17. [DOI: https://dx.doi.org/10.4081/jlimnol.2013.s2.e1]

33. Moss, S.J.; Wilson, M.E.J. Biogeographic implications of the Tertiary palaeogeographic evolution of Sulawesi and Borneo. Biogeography and Geological Evolution of SE Asia; Hall, R.; Holloway, J.D. Backhuys Publishers: Leiden, The Netherlands, 1998; pp. 133-163.

34. Yu, Y.; Harris, A.J.; Blair, C.; He, X. RASP (Reconstruct Ancestral State in Phylogenies): A tool for historical biogeography. Mol. Phylogenet. Evol.; 2015; 87, pp. 46-49. [DOI: https://dx.doi.org/10.1016/j.ympev.2015.03.008] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25819445]

35. Mirarab, S.; Warnow, T. ASTRAL-II: Coalescent-based species tree estimation with many hundreds of taxa and thousands of genes. Bioinformatics; 2015; 31, pp. i44-i52. [DOI: https://dx.doi.org/10.1093/bioinformatics/btv234] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26072508]

36. Kass, R.E.; Raftery, A.E. Bayes factors. J. Am. Stat. Assoc.; 1995; 90, pp. 773-795. [DOI: https://dx.doi.org/10.1080/01621459.1995.10476572]

37. Čandek, K.; Agnarsson, I.; Binford, G.J.; Kuntner, M. Biogeography of the Caribbean Cyrtognatha spiders. Sci. Rep.; 2019; 9, pp. 1-14. [DOI: https://dx.doi.org/10.1038/s41598-018-36590-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30674906]

38. Gillespie, R.G.; Baldwin, B.G.; Waters, J.M.; Fraser, C.I.; Nikula, R.; Roderick, G.K. Long-distance dispersal: A framework for hypothesis testing. Trends Ecol. Evol.; 2012; 27, pp. 47-56. [DOI: https://dx.doi.org/10.1016/j.tree.2011.08.009]