Phylogenomic analyses reveal a Gondwanan origin and repeated out of India colonizations into Asia by tarantulas (Araneae: Theraphosidae)

Abstract

The study of biogeography seeks taxa that share a key set of characteristics, such as timescale of diversification, dispersal ability, and ecological lability. Tarantulas are ideal organisms for studying evolution over continental-scale biogeography given their time period of diversification, their mostly long-lived sedentary lives, low dispersal rate, and their nevertheless wide circumtropical distribution. In tandem with a time-calibrated transcriptome-based phylogeny generated by PhyloBayes, we estimate the ancestral ranges of ancient tarantulas using two methods, DEC+j and BBM, in the context of their evolution. We recover two ecologically distinct tarantula lineages that evolved on the Indian Plate before it collided with Asia, emphasizing the evolutionary significance of the region, and show that both lineages diversified across Asia at different times. The most ancestral tarantulas emerge on the Americas and Africa 120 Ma–105.5 Ma. We provide support for a dual colonization of Asia by two different tarantula lineages that occur at least 20 million years apart, as well as a Gondwanan origin for the group. We determine that their current distributions are attributable to a combination of Gondwanan vicariance, continental rafting, and geographic radiation. We also discuss emergent patterns in tarantula habitat preferences through time.

Introduction

Due to the historically close geographies of several continental land masses prior to the breakup of the supercontinents Laurasia and Gondwana, many related taxa are presently found on opposite sides of the world. For instance, it has been shown that rodents and primates (De Oliveira, Molina & Marroig, 2009) in the Neotropics likely arrived there from Africa. Several mass faunal exchanges between continents have occurred when two land masses eventually meet, such as the “Great American Interchange” (Marshall et al., 1982) after the formation of the Isthmus of Panama (O’Dea et al., 2016), the dispersal of several taxa “out of India” and into Asia/Indomalaya (Karanth, 2006), and exchanges between Africa and Eurasia (Hedges, 2001; Springer et al., 2011). If we are to understand the evolutionary patterns involved in the diversification of taxa and the generation of biodiversity, it becomes important to consider the roles that continental drift may play, particularly in cases where evolutionary novelties have emerged in the diversification process.

Studies that use model taxa to understand biogeographic patterns are limited by the fact that the ideal model organism should have both a low dispersal ability and a wide (if not global) distribution. Given their low motility and wide array of dispersal strategies, plants have often been the focus of biogeographic studies (e.g., Zhou et al., 2018; McGlone et al., 2018; Fariña et al., 2018; Echeverría-Londoño et al., 2018). Mygalomorph spiders from the family Theraphosidae, also known as tarantulas, are also useful models for understanding biogeographic patterns. They possess several intriguing adaptations that have received significant attention in recent years, such as their stridulatory capacity (Galleti-Lima & Guadanucci, 2018; Galleti-Lima & Guadanucci, 2019), urticating bristle evolution (Bertani, 2003; Kaderka et al., 2019; Foley et al., 2019), colouration (Foley, Saranathan & Piel, 2020), and venom compositions (Santana et al., 2017; Lüddecke, Vilcinskas & Lemke, 2019). Robust, subfamily-level phylogenies are late in coming (Lüddecke et al., 2018; Foley et al., 2019), but these phylogenies have verified the monophylies of all studied subfamilies, except Ischnocolinae (with Selenogyrinae being the only subfamily not represented by these studies). It has only recently become possible to study tarantula adaptations in an explicitly phylogenetic context.

This diverse suite of intriguing adaptations undoubtedly contributed to the success of these creatures. They are quite widespread and are found throughout the subtropical regions of every continent (Gallon, 2000; World Spider Catalog, 2020). However, while many mature male tarantulas are known to wander in search of females (Prentice, 1992; Pérez-Miles et al., 2005; Pérez-Miles et al., 2007), most juvenile and female tarantulas rarely venture far from their retreats (Yáñez & Floater, 2000; Schultz & Schultz, 2009). Spiderlings seldom disperse far from the burrows of their mothers, leading to dense aggregations of burrows in relatively small areas (Reichling, 2000; Shillington & McEwen, 2006). There are no reports of ballooning tarantulas like in other mygalomorph families (Coyle, 1983; Coyle et al., 1985; Buzatto, Haeusler & Tamang, 2021), but some mygalomorphs are known to be highly divergent across small areas (Starrett & Hedin, 2007). These behaviours do not portend that tarantulas would be successful dispersers, yet they have spread across the globe and have colonized strikingly different ecological niches.

Previous studies estimate that tarantulas emerged between 150 Ma–71 Ma (Bond et al., 2014) or ∼107 Ma (Opatova et al., 2020), which is compatible with a Gondwanan origin. Indeed, some tarantulas (Selenocosmiinae) are suggested to be North Gondwanan taxa (West, Nunn & Hogg, 2012). Given their habitat specificity and low dispersal ability, yet wide geographic distribution, tarantulas are good model organisms for addressing biogeographic questions. In particular, two distinct lineages of Asian tarantulas were recovered in the most recent subfamily-level tarantula phylogeny (Foley et al., 2019). Estimates surrounding the timing of the collision between the Indian Subcontinent and Asia vary, but most studies agree on a late Paleocene/early Eocene timescale around 55 Ma (Hu et al., 2016; Zheng & Wu, 2018; Westerweel et al., 2019). Migrations from the Indian Subcontinent to Asia are known to have played a significant role in generating the high faunal diversity observed in Asia today (Karanth, 2006; Morley et al., 2016; Morley, 2018; Garg & Biju, 2019), even amongst spiders (Li, Shao & Li, 2020) and other invertebrates (Joshi, Karanth & Edgecombe, 2020). It has also been shown that other spiders of a comparable age originated on the Gondwanan supercontinent, i.e., South America and Africa, before eventually colonizing Asia and Oceania (Chamberland et al., 2018). We expect to recover a similar pattern for tarantulas.

Here, we test the Gondwanan origin hypothesis for tarantulas using a broad sampling of molecular data and taxa. We also hypothesize that Asian tarantulas may have originated on the Indian Subcontinent and subsequently dispersed throughout Asia. We test this by inferring the historical patterns that contributed to present-day theraphosid distributions. We use a subfamily-level chronogram of Theraphosidae and estimate their ancestral ranges under two different models, allowing us to interpret the node ages in the context of geographical events. We then discuss a timeline of events that may have led to the emergence of theraphosid subfamilies, with emphasis on those found in Asia, and present some hypotheses on how those events impacted their current distributions.

Methods Data acquisition

Although theraphosids are the primary focus of this study, the poor fossil record in this group required that we expand our sampling to include other Mygalomorphae with richer fossil records so that these could assist in molecular clock calibration. Only two tarantula fossils are known, and neither can be firmly placed onto the tarantula phylogeny; one from Burmese amber dating to ∼100 Mya (Wunderlich & Müller, 2020), and one from Chiapas amber dating to somewhere between the late Oligocene and mid Miocene (possibly as young as 16 Mya, Dunlop, Harms & Penney, 2008). Inferring calibration points from prior studies (i.e., secondary calibration) has also been discouraged (Schenk, 2016).

Transcriptome data for 29 theraphosids and 18 other mygalomorphs was obtained from the publicly available sequence read archive (SRA) database (http://www.ncbi.nlm.nih.gov/sra). Additionally, to improve the rooting of the family, we included a new transcriptome from a field-collected barychelid, Rhianodes atratus (collected under permit number NP/RP18-046, issued by the National Parks Board of Singapore). RNA was extracted using a TRIzol total RNA extraction protocol (Simms, Cizdziel & Chomczynski, 1993). The RNA was subsequently sent to a commercial company (NovogeneAIT) for paired-end sequencing on an Illumina HiSeq 4000. The resulting transcriptome is available on NCBI’s Sequence Read Archive (SRA) under bioproject number PRJNA635363 (accession GIXJ00000000). Accession numbers for a total of 48 transcriptomes are included in Table S1. Assembly data for the Rhianodes atratus (Barychelidae) transcriptome is available in Table S2.

A core ortholog variant to obtain a DNA supermatrix

Transcriptomes were assembled using Trinity v2.6.6 (Grabherr et al., 2011), and protein coding regions were predicted using TransDecoder v5.5.0 (Haas et al., 2013). The core-ortholog (OG) pipeline as per Garrison et al. (2016) was implemented for spider-specific ortholog selection across our dataset. The initial set of core orthologs was generated from a selection of publicly available transcriptomes, as per Cheng & Piel (2018), and a set of 4,446 profile hidden Markov models (pHMMs) was generated following their methods. HaMStR v13.1 (Ebersberger, Strauss & Von Haeseler, 2009) inferred orthology between these pHMMs and our dataset. Groups of orthologous genes were then pooled together and subjected to several filtering steps.

Software used during the filtering process include the following: MAFFT (Katoh et al., 2005) for alignment during filtering; ALISCORE (Misof & Misof, 2009; Kück et al., 2010) for alignment trimming; ALICUT (Kück, 2009) for excising ambiguous regions; Infoalign (Rice, Longden & Bleasby, 2000) for generating consensus sequences from OGs. Generally, the filtering process and criteria was the same as in Foley et al. (2019) with one key exception—instead of removing any OG that was sampled for half (or fewer) of the taxa included, we chose to only retain OGs that were present in all taxa, discarding the rest. Despite this stringency, an adequate set of 743 OG alignments remained, each of which included all 48 taxa. To obtain a more neutral, granular dataset, these protein alignments were converted back to DNA alignments using localized tBLASTn searches for each taxon (Camacho et al., 2009) against their corresponding Trinity outputs, in tandem with a set of custom Perl scripts for parsing. Gene trees for each of the 743 alignments were estimated using RAxML (Stamatakis, 2014). These departures from previous filtering criteria were informed by the downstream needs of our analysis pipeline.

To exclude OGs with conflicting signals, we identified divergent gene trees using Robinson-Foulds (RF) metrics (Robinson & Foulds, 1981) in the R package “phytools” (Revell, 2012). Given that the theraphosid tree from Foley et al. (2019) received 100% bootstrap support at almost all nodes, we used this as the reference tree to compare with each gene tree. In order to calculate RF values, each of our 743 gene trees was pruned to match the set of taxa in the reference tree. Any OG DNA alignment whose corresponding gene tree received an RF score of ≥ 6 was discarded, leaving 125 OGs.

Next, TranslatorX (Abascal, Zardoya & Telford, 2010) was used to realign the DNA sequences for the 125 OGs according to our protein alignments. Visualization of OGs in DAMBE v7 (Xia & Xie, 2001) showed that phylogenetically informative content was being distorted due to sequence saturation (Xia et al., 2003; Philippe et al., 2011), so every third base was subsequently removed to mitigate this effect. Lastly, FASconCAT (Kück & Meusemann, 2010) was used to concatenate all 125 OGs into a supermatrix, and GBLOCKS (Castresana, 2002) was subsequently used to trim and excise ambiguous regions. The result was a 96% complete DNA alignment of 48 taxa and 89,302 characters.

Time calibration with PhyloBayes

BEAST v1.10.4 (Drummond & Rambaut, 2007) was used to generate an initial starting phylogeny for time calibration, running for 100 million generations with 20% of the trees discarded as burnin. Divergence times were estimated in PhyloBayes v4.1 (Lartillot, Lepage & Blanquart, 2009) using this tree, the DNA supermatrix, and four fossils to infer calibration points. We followed Opatova et al. (2020) in treating the oldest mygalomorph fossil (Rosamygale grauvogely, Selden & Gall, 1992) as an Avicularioidea crown group, and use a subset of their fossil calibrations.

The first calibrated node corresponds to the Avicularioidea-Atypoidea split. The Rosamygale grauvogely fossil dated to 242 Million years ago (Ma), which represents the minimum bound for this node, while 323 Ma serves as the maximum bound and corresponds to the oldest age of the Bashkirian stage from which the oldest known spider fossil (Arthrolycosa sp.) hails (Selden et al., 2014; Garwood et al., 2016). Secondly, the split between Nemesioidina and Bemmeridae/Theraphosidae was assigned a minimum bound of 125 Ma based on the Cretamygale chasei fossil (Selden, 2002), and a maximum bound of 242 Ma based on Rosamygale grauvogely. Lastly, the split between Antrodiaetidae and Atypidae was assigned a minimum bound of 100 Ma based on a fossil of Ambiortiphagus ponomarenkoi (Eskov & Zonstein, 1990), with Rosamygale grauvogely once again providing the maximum bound of 242 Ma for this node.

Four independent PhyloBayes chains were each run for c. 24,000 cycles under an autocorrelated lognormal model with a GTR matrix, and otherwise default parameters. Several advantages of autocorrelated lognormal models over uncorrelated rates models have been demonstrated for a variety of datasets (Lepage et al., 2007; Paradis, 2013), and they have been shown to yield more plausible results in similar studies to ours (Sharma et al., 2018). Burn-in was adjusted to 10% to maximize the number of estimated sample size values that exceeded 200, yielding a total of 86,933 cycles across all chains. Divergence times were estimated using the “readdiv” command in PhyloBayes for each chain, and the results were summarized with TreeAnnotator. The resulting time-calibrated mygalomorph phylogeny is available in Fig. S1, and we proceeded with the tarantula subsection of this phylogeny for further analyses.

Biogeographic analysis

The tarantulas included in this study were assigned to one or more of the following five biogeographic ranges as determined by the World Spider Catalog (2020): (A) Indian Subcontinent; (B) Asia (non-Indian Indomalaya); (C) East of the Wallace Line; (D) Americas; and (E) Africa. Taxon range codings are included in Table S1. We also accounted for proximity between areas by excluding certain combinations of ancestral area (e.g., Asia plus Americas. Full list in Table S3). These areas were chosen based on the biogeographic realms first proposed by Udvardy (1975), but we treat India separately to the remainder of Indomalaya given that we expect the age of theraphosids to exceed the collision of India and Asia.

The tarantula subsection (with outgroups) of our tree from PhyloBayes served as input for RASP (Yu et al., 2015), which estimated ancestral range distributions under two different models that were chosen via model testing (Table 1) in BioGeoBEARS (Matzke, 2013). Despite BAYAREALIKE+j being predicted as the most favorable model, its assumption of a static geological history (Landis et al., 2013) rendered it unsuitable. Hence, we ran the analysis using our second highest scoring model—the dispersal-extinction cladogenesis with jump dispersal (DEC+j) model (Ree & Smith, 2008). However, such “+j” models have been divisive, with some criticizing the “+j” parameter (Ree & Sanmartín, 2018) and others advocating for its inclusion (McDonald-Spicer et al., 2019; Klaus & Matzke, 2019). Therefore, Bayesian Binary MCMC (BBM) was added as a second method, and was ran under the JC69 model for 1 million generations across ten chains with sampling every 100 generations (Ali et al., 2011). As the model testing found no statistically significant difference between DEC and DEC+j (Table 1), DEC was added as a third method. RASP was also used to calculate an event matrix for node 21.

Table 1: Model testing. Model appropriateness as determined by model testing in BioGeoBEARS. While BAYAREALIKE+j emerges as the most appropriate, its assumption of a static geologic history led us to choose DEC+j instead (as indicated by bold font). Asterisks indicate a significant difference between the appropriateness of a model with and without the +j. Hence, we also include DEC in supplementary Fig. S3.

Model	LnL	AICc	AICc_wt	Δ−j/+j
DEC	−34.83	74.13	0.22
DEC+J	−33.2573.51	0.3
DIVALIKE	−38.51	81.51	0.0055	*
DIVALIKE+J	−34.18	75.35	0.12
BAYAREALIKE	−48.13	100.7	3.70E−07	*
BAYAREALIKE+J	−33.1	73.21	0.35

DOI: 10.7717/peerj.11162/table-1

Results Ancestral ranges and ages of theraphosid clades

Our DEC+j (Fig. 1) and BBM (Fig. S3) results are largely in agreement, disagreeing at only at two nodes (per Fig. 1)—node 20 (which DEC+j estimates as present on both Africa and the Americas, whereas BBM estimates the Americas only), and node 28 (which DEC+j estimates as present on Asia only, whereas BBM estimates it as present on Asia and East of the Wallace line). Full range probabilities for DEC+j and BBM are given in Figs. S2 and S3 respectively. The DEC reconstruction is available in Fig. S4, and although there are some differences between this method and the others, the larger concepts discussed are not impacted.

Figure 1: Ancestral ranges as estimated under the DEC+j model. Subfamilies are given for each species, along with a timescale in millions of years. Node numbers are provided, and each node is color coded to correspond with the areas highlighted on the map. Letters on the map correspond to range codings given in the Methods section. The stars on the tree represent independent colonizations of Asia by theraphosids, having arrived there from the Indian Subcontinent. Map credit: https://mapchart.net/, 2021. Licensed under CC BY 4.0 SA. DOI: 10.7717/peerj.11162/fig-1

Table 2 contains divergence time estimates as 95% HPDs for the theraphosid subsection of the mygalomorph phylogeny, while Table S4 contains those estimates for all nodes. Our earliest theraphosid node (20) is recovered 120 Ma–115.5 Ma in the Americas (Fig. 1), and is found sister to Rhianodes atratus, indicating that Barychelidae are sister to Theraphosidae. When considered alongside the divergence times, these estimations are consistent with a Gondwanan origin for Theraphosidae.

Table 2: Divergence time estimates. Divergence times with error margins per 95% HPDs as estimated by PhyloBayes. Times are given in millions of years and are rounded to the closest 0.5 Ma. Node numbers as per Fig. 1.

Node number	Estimated age	95% HPD Max	95% HPD Min
20	119	120	115.5
21	111	112	107.5
22	57	58	55.5
23	45	45.5	43.5
24	109	111	106.5
25	107	108	103.5
26	98	99	95
27	57	57.5	55
28	47	47.5	45.5
29	99	100	96
30	72	73	69.5
31	28	28.5	27
32	16	16.5	15
33	71.5	72	69
34	34.5	35	33
35	23.5	24	22
36	100	101	97
37	95	96	92.5
38	89	90	86.5
39	77	78	74.5
40	42	42.5	40.5
41	39	39.5	37.5
42	18.5	18.5	17.5
43	41.5	42	39.5
44	55.5	56	53.5
45	50.5	51	48.5
46	29.5	30	28.5

DOI: 10.7717/peerj.11162/table-2

The Indian Plate is first recovered as an ancestral area at node 25 108 Ma–103.5 Ma. All ancestral range estimations recover two ecologically distinct Indian Plate lineages descending from node 25. The first such node (26) emerges 99 Ma–95 Ma, yielding several ground-dwelling taxa. It diverges into the Indian Thrigmopoeinae lineage, and the ancestor to Selenocosmiinae (node 27), which emerges 57.4 Ma–55 Ma. Node 27 represents the first appearance of Asia as an ancestral area in both DEC+j and BBM, albeit alongside the Indian Subcontinent (Fig. 1). This node then diverges into the Haplocosmia lineage, and the ancestor to Selenocosmia and Phlogiellus, which also emerges as present east of the Wallace line around 47.5 Ma-45.5 Ma (Fig. 1).

The second descendant of the first Indian Subcontinent node (25) corresponds to the last common ancestor of Harpactirinae and Ornithoctoninae (node 29), emerging 100 Ma–96 Ma. This node subsequently diverges into the ancestor of Harpactirinae and Stromatopelminae (node 30), and yet another Indian Subcontinent node corresponding to the ancestors of the primarily arboreal Ornithoctoninae, and arboreal Poecilotheriinae (node 33). After the divergence of Ornithoctoninae and Poecilotheriinae, we see Asia represented as an ancestral area for the second time in both DEC+j and BBM as the ancestor to Ornithoctoninae (node 34) emerges 35 Ma–33 Ma. DEC, however, places the ancestor to Poecilotheriinae and Ornithoctoninae on both the Indian Plate and Asia (Fig. S4).

The ancestor to the primarily African Eumenophorinae, and the Catumiri lineage (node 21) emerges in the Americas 112 Ma–107.5 Ma in all three analyses, with DEC also placing it in Africa. The event matrices of all analyses for this node are given in Table 3.

Table 3: Event matrix at node 21. The ancestor to both Eumenophorinae and the Catumiri lineage, node 21, can provide key insights into tarantula biogeography. DEC+j and BBM favor a dispersal event influencing this node, which is more consistent with members of Eumenophorinae crossing from South America to Africa. However, DEC favors vicariance, which would attribute the present-day distribution to continental drift.

Event	DEC+j	BBM	DEC
Dispersal	2	2	0
Vicariance	1	1	1
Extinction	0	0	0

DOI: 10.7717/peerj.11162/table-3

Discussion The importance of India and the dual colonizations of Asia

The “out of India” hypothesis (McKenna, 1973) has been accepted as a potential distribution pattern for the radiations of a diverse array of taxa (Karanth, 2006; Conti et al., 2002; Joshi, Karanth & Edgecombe, 2020). Our results are consistent with two independent “out of India” dispersals of tarantulas into Asia (Figs. 1 and 2). This “dual colonization of Asia via India” pattern has also been demonstrated in centipedes (Joshi, Karanth & Edgecombe, 2020). Both of our Asian tarantula lineages were estimated to arrive there during timescales consistent with Li, Shao & Li (2020), who determined that Ochyroceratid spiders likely arrived in Asia via India 55 Ma–38 Ma, and also with the India/Asia collision timescales (Hu et al., 2016; Westerweel et al., 2019).

Figure 2: The dual colonization of Asia. Two independent tarantula lineages emerge out of India. (A) The first theraphosids on the Indian Plate diverge into the ancestor of Selenocosmiinae and Thrigmopoeinae (98 Ma) and into the ancestor of Ornithoctoninae and Poecilotheriinae (99 Ma); (B) present-day distributions of Selenocosmiinae. They entered Asia via India 57 Ma and crossed the Wallace line sometime after 47 Ma. Some Selenocosmiinae (e.g., Haplocosmia nepalensis) may also be found on the Indian Subcontinent today, where Thrigmopoeinae remain; (C) present-day distributions of Ornithoctoniinae, which entered Asia via India 34.5 Ma but did not manage to cross the Wallace line. Meanwhile, Poecilotheriinae remain on the Indian Subcontinent. Map credit: https://dinosaurpictures.org/, 2021. Licensed under CC BY 4.0 SA. DOI: 10.7717/peerj.11162/fig-2

All three analyses place the ancestor to Selenocosmiinae (node 27) as present on the Indian Subcontinent and Asia 57.5 Ma–55 Ma (Fig. 1)—the first entry of theraphosids into Asia per DEC+j and BBM. Furthermore, both of these analyses recover the ancestor to Selenocosmia and Phlogiellus (node 28) in both Asia and east of the Wallace line, with DEC recovering the node as Asia-only. Interestingly, despite our analyses suggesting a Gondwanan origin for Theraphosidae, this pattern suggests that tarantulas were not always present in Oceania, and instead is consistent with Selenocosmiinae having diversified across Asia, eventually crossing the Wallace line sometime after the India/Asia collision (possibly as early as 47 Ma, Fig. S3), while the terrestrial Thrigmopoeinae remained in India (Fig. 2). The second colonization event of Asia by theraphosids occurred with the common ancestors to Ornithoctoninae (node 34, Fig. 1), which emerged in Asia at least 20 million years later (35 Ma–33 Ma), while Poecilotheriinae remained in India (Fig. 2). DEC suggests that this lineage colonized Asia even earlier, placing the ancestors to Poecilotheriinae (node 33) on both Asia and the Indian Plate 72 Ma–69 Ma, but this seems too early to be plausible as (i) it predates the collision of the Indian Plate into the rest of Asia, and (ii) no extant members of Poecilotheriinae exist in the rest of Asia, and DEC did not predict an extinction event at this node.

Our results indicate that both of these Asian lineages diverged while the Indian Plate was still rafting towards Asia, as the collision would likely not have happened before node 25 diverges to yield nodes 26 and 29 (Table 2), both of which also emerge on the Indian Plate in our DEC+j and BBM analyses (Fig. 1 and Fig. S3), while DEC also places node 29 on the African continent in addition to the Indian Plate. Interestingly, the two lineages also appear to be ecologically divergent. Selenocosmiinae is a subfamily that almost exclusively contains terrestrial members, whereas Ornithoctoninae is a largely arboreal subfamily with terrestrial members.

A theraphosid fossil was recently uncovered in Burmese amber, dating to ∼100 Ma (Wunderlich & Müller, 2020). While the identity of the fossil cannot adequately be established (though the authors did speculate that it may represent an ancestor to Selenocosmiinae), Burmese amber is thought to be of Gondwanan origin (Oliveira et al., 2016; Poinar Jr, 2018). The location and age of this fossil could plausibly correspond to one of the ancestral tarantula groups at either node 21 or 24, and confers some support towards the validity of our ancestral area estimations and dating.

The Catumiri lineage as a biogeographic relic, and a journey to Africa for Eumenophorinae

Our estimate for the age of our most ancestral theraphosid node (20) was 120 Ma–105.5 Ma (Table 2), which is consistent with previous estimates based on fewer ingroups that dated the divergence of tarantulas to occur 150 Ma–71 Ma (Bond et al., 2014) or ∼107 Ma (Opatova et al., 2020). All three analyses recover this node in the Americas, but DEC also places it on Africa. South America finally split from Africa 105 Ma–100 Ma (Heine, Zoethout & Müller, 2013; Matthews et al., 2016; Olyphant, Johnson & Hughes, 2017). While oceanic crossings were still required, it has been suggested that dispersal between both continents may have been facilitated via the Rio Grande Rise and Walvis Ridge up until ∼40 Ma, as these might have shortened the terrestrial distance between both continents (De Oliveira, Molina & Marroig, 2009).

The ancestor to the primarily African Eumenophorinae and the Catumiri lineage (node 21, Fig. 1) is recovered by all analyses in the Americas 112 Ma–107.5 Ma (with DEC also placing it in Africa)—a period when South America and Africa were likely still connected. As a member of the taxonomically troubled Ischnocolinae, Catumiri and its relatives have been considered to be among the most basal theraphosids (Raven, 1985; Schmidt, 2003; Guadanucci, 2014). The findings from Opatova et al. (2020) and Foley et al. (2019) seem to agree with this with respect to Catumiri, which is recovered as an early-branching theraphosid in both works. Catumiri was also the only neotropical taxon from Foley et al. (2019) that was not recovered in an otherwise monophyletic group of neotropical theraphosids. With a potential maximum age of just under 107.5 Ma (Table 2), we suggest that the Catumiri lineage represents an ancient biogeographical relic that diversified and remained in South America while the continents drifted apart—potentially even before South America split from Africa. This could also plausibly explain why Ischnocolinae (the subfamily in which Catumiri is included) was the only paraphyletic subfamily recovered by Foley et al. (2019).

Meanwhile, the ancestor to Eumenophorinae emerges solely in Africa in all three analyses 58 Ma–55.5 Ma. This result leaves two potential biogeographic origins for Eumenophorinae: (i) a vicariant event caused by the split between South America and Africa, which resulted in them being separated from the Catumiri lineage; or (ii) a dispersal event, where some of the ancestors of the Catumiri lineage crossed from South America to Africa and subsequently diversified. Per the event matrix in Table 3, the first scenario is supported by DEC (Fig. S4), which recovers the ancestors to both groups on both South America and Africa, whereas DEC+j and BBM favor a dispersal event, which is more consistent with a scenario where Eumenophorinae crossed back to Africa from South America, which may still have been possible as per the results from De Oliveira, Molina, and Marroig (2009).

A second emergence of tarantulas on Africa occurs at the ancestor to Harpactirinae and Stromatopelminae 73 Ma–69.5 Ma (node 30, Fig. 1), having apparently crossed back via the Indian Plate. The upward bound of this arrival can be attributed to the minimum age of its parent node (29), which is estimated at 96 Ma (Table 2), and a minimum of 69.5 Ma (minimum age of the ancestor to Harpactirinae and Stromatopelminae). Crossing between Africa and the Indian Plate via Madagascar could have been possible up until around 90 Ma to 85 Ma (Ali & Aitchison, 2008), and our estimates hence afford a short window for this node to cross to Africa. There are a handful of extant Malagasy theraphosids (World Spider Catalog, 2020) that could provide further support for this link, though sample material from these taxa was not available. This “back to Africa” event would represent an intriguing departure from typical patterns associated with Gondwanan vicariance.

Tarantula radiations and lifestyle switches

There are three nodes in particular (20, 35, and 37, Fig. 1) where the two descendent lineages demonstrate high fidelity to habitat choices. For node 20, one descendant lineage corresponds to the fossorial Selenocosmiinae and Thrigmopoeinae, while the other corresponds to the arboreal Poecilotheriinae and Ornithoctoninae (except Haplopelma, which is fossorial). Node 35 splits into Harpactirinae (fossorial) and Stromatopelminae (arboreal). Node 37 splits into Theraphosinae (fossorial), and Aviculariinae, Psalmopoeinae (both arboreal) and Schismatothelinae (fossorial). These habitat choices are maintained throughout radiations across large areas, indicating that, based on our phylogeny and sampling, tarantula lineages may be demonstrating niche conservatism in showing high fidelity to lifestyle and habitat choices.

This idea is also relevant to our two independent lineages of tarantulas that colonized Asia. After the terrestrial Selenocosmiinae diversified across Asia, there may have been an opportunity for Ornithoctoninae to capitalise on an open arboreal niche, and subsequently colonize Asia and further diversify while largely retaining this newly acquired arboreal preference. The two lineages that did not colonize Asia/remained on the Indian Plate are also ecologically divergent, with Poecilotheriinae being an arboreal subfamily and Thrigmopoeinae occupying terrestrial niches, suggesting that their lifestyle differences may have facilitated their co-diversification throughout the Indian Subcontinent. The divergence times imply that the arboreal Poecilotheriinae could be as old as 69 Ma, and the terrestrial Thrigmopoeinae could be as old as 95 Ma (Table 2), which implies that this co-diversification might have occurred before the Indian Plate fully collided with Asia.

Perhaps these radiations can be attributed to an ancient switch in lifestyle that each ecologically distinct subfamily to become successful by exploiting different ecological niches. Although the subfamilies in the phylogeny used here are not represented by all members, we expect that the addition of further taxa would verify this idea, given that the niches of most higher theraphosid taxa appear to be conserved over large timescales, and that lifestyle changes are observed at the base of many theraphosid radiations. With the advent of further data from more representative taxa, we encourage future studies to test whether these lifestyle switches directly precede evolutionary rate shifts.

Conclusions

The Indian Plate played a pivotal role in the history of Theraphosidae and facilitated the colonization of Asia by two ecologically distinct tarantula lineages at least 20 million years apart. Ancient tarantulas appear to have undergone several diversifications on India while it was still rafting, affirming the evolutionary significance of the subcontinent. We show that tarantulas demonstrate pronounced niche conservation, which could explain the ecological divergences observed in the two independent Indian subcontinent lineages that colonized Asia. We propose a Gondwanan origin for Theraphosidae, although their present-day distributions suggest that continental rafting and geographic radiations facilitated their colonizations of subsequent landmasses.

Additional Information and Declarations

Competing Interests

The authors declare there are no competing interests.

Author Contributions

Saoirse Foley conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.

Henrik Krehenwinkel, Dong-Qiang Cheng and William H. Piel conceived and designed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft.

Field Study Permissions

The following information was supplied relating to field study approvals (i.e., approving body and any reference numbers):

The National Parks Board of Singapore granted permission to collect Rhianodes atratus (Permit number NP/RP18-046).

DNA Deposition

The following information was supplied regarding the deposition of DNA sequences:

The resulting transcriptome is available on NCBI’s Sequence Read Archive (SRA): PRJNA635363, GIXJ00000000.

Accession numbers for a total of 48 transcriptomes and assembly data for the Rhianodes atratus (Barychelidae) transcriptome is available in the Supplemental Files.

Data Availability

The following information was supplied regarding data availability:

Data are available in the Supplemental Files.

Funding

The computational results reported in this study were supported by a Yale-NUS College Shared Instrumentation Grant (IG15-SI101) awarded to William H. Piel. William H. Piel was also supported by the South East Asian Biodiversity Genomics (SEABIG) centre, the Singapore Ministry of Education (grant number MOE2016-T2-2-137), and by Yale-NUS College (grant numbers IG14-SI002 and R-607-265-052-121). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

Abascal F, Zardoya R, Telford MJ. 2010. TranslatorX: multiple alignment of nucleotide sequences guided by amino acid translations. Nucleic Acids Research 38:W7-W13

Ali JR, Aitchison JC. 2008. Gondwana to Asia: plate tectonics, paleogeography and the biological connectivity of the Indian sub-continent from the Middle Jurassic through latest Eocene (166–35 Ma) Earth-Science Reviews 88(3–4):145-166

Ali SS, Yu Y, Pfosser M, Wetschnig W. 2011. Inferences of biogeographical histories within subfamily Hyacinthoideae using S-DIVA and Bayesian binary MCMC analysis implemented in RASP (Reconstruct Ancestral State in Phylogenies) Annals of Botany 109(1):95-107

Bertani R. 2003. Release of urticating hairs by Avicularia versicolor (Walckenaer, 1837)(Araneae, Theraphosidae) Bulletin de la Société entomologique de France 44:29-92

Bond JE, Garrison NL, Hamilton CA, Godwin RL, Hedin M, Agnarsson I. 2014. Phylogenomics resolves a spider backbone phylogeny and rejects a prevailing paradigm for orb web evolution. Current Biology 24(15):1765-1771

Buzatto BA, Haeusler L, Tamang N. 2021. Trapped indoors? Long-distance dispersal in mygalomorph spiders and its effect on species ranges. Journal of Comparative Physiology A 207:1-14

Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. 2009. BLAST+: architecture and applications. BMC Bioinformatics 10(1):421

Castresana J. 2002. GBLOCKS: selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Copyrighted by J. Castresana, EMBL.

Catalog WS. 2020. The World Spider Catalog, Natural History Museum Bern, online at http://wsc.nmbe.ch, version 20.0. (accessed26 November 2019 ) software

Chamberland L, McHugh A, Kechejian S, Binford GJ, Bond JE, Coddington J, Dolman G, Hamilton CA, Harvey MS, Kuntner M. 2018. From Gondwana to GAAR landia: evolutionary history and biogeography of ogre-faced spiders (Deinopis) Journal of Biogeography 45(11):2442-2457

Cheng D-Q, Piel WH. 2018. The origins of the Psechridae: web-building lycosoid spiders. Molecular Phylogenetics and Evolution 125:213-219

Conti E, Eriksson T, Schönenberger J, Sytsma KJ, Baum DA. 2002. Early Tertiary out-of-India dispersal of Crypteroniaceae: evidence from phylogeny and molecular dating. Evolution 56(10):1931-1942

Coyle FA. 1983. Aerial dispersal by mygalomorph spiderlings (Araneae, Mygalomorphae) Journal of Arachnology 11:283-286

Coyle FA, Greenstone MH, Hultsch A-L, Morgan CE. 1985. Ballooning mygalomorphs: estimates of the masses of Sphodros and Ummidia ballooners(Araneae: Atypidae, Ctenizidae) Journal of Arachnology 13(3):291-296

de Oliveira FB, Molina EC, Marroig G. 2009. Paleo-geography of the South Atlantic: a route for primates and rodents into the New World? In: Garber PA, Estrada A, Strier KB, eds. South American primates:comparative perspectives in the study of behavior, ecology, and conservation. New York, NY: Springer. 55-68

Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7(1):214

Dunlop JA, Harms D, Penney D. 2008. A fossil tarantula (Araneae: Theraphosidae) from Miocene Chiapas amber, Mexico. Revista Ibérica de Aracnología 15:9-17

Ebersberger I, Strauss S, Von Haeseler A. 2009. HaMStR: profile hidden markov model based search for orthologs in ESTs. BMC Evolutionary Biology 9(1):157

Echeverría-Londoño S, Enquist BJ, Neves DM, Violle C, Boyle B, Kraft NJ, Maitner BS, McGill B, Peet RK, Sandel B. 2018. Plant functional diversity and the biogeography of biomes in North and South America. Frontiers in Ecology and Evolution 6:219

Eskov KY, Zonstein S. 1990. First Mesozoic mygalomorph spiders from the Lower Cretaceous of Siberia and Mongolia, with notes on the system and evolution of the infraorder Mygalomorphae (Chelicerata: Araneae) Neues Jahrbuch für Mineralogie, Geologie und Paläontologie, Abhandlungen 178:325-368

Fariña JM, He Q, Silliman BR, Bertness MD. 2018. Biogeography of salt marsh plant zonation on the Pacific coast of South America. Journal of Biogeography 45(1):238-247

Foley S, Lüddecke T, Cheng D-Q, Krehenwinkel H, Künzel S, Longhorn SJ, Wendt I, Von Wirth V, Tänzler R, Vences M, Piel WH+1 more. 2019. Tarantula phylogenomics: a robust phylogeny of deep theraphosid clades inferred from transcriptome data sheds light on the prickly issue of urticating setae evolution. Molecular Phylogenetics and Evolution 140:106573

Foley S, Saranathan V, Piel WH. 2020. The evolution of colouration and opsins in tarantulas. Proceedings of the Royal Society B 287(1935):20201688

Galleti-Lima A, Guadanucci JPL. 2018. Morphology of setae on the coxae and trochanters of theraphosine spiders (Mygalomorphae: Theraphosidae) The Journal of Arachnology 46(2):214-226

Galleti-Lima A, Guadanucci JPL. 2019. Comparative morphology of stridulating setae of Theraphosinae (Araneae: Theraphosidae) Zoologischer Anzeiger

Gallon RC. 2000. The natural history of tarantula spiders. Natural History 8:18

Garg S, Biju S. 2019. New microhylid frog genus from Peninsular India with Southeast Asian affinity suggests multiple Cenozoic biotic exchanges between India and Eurasia. Scientific Reports 9(1):1906

Garrison NL, Rodriguez J, Agnarsson I, Coddington JA, Griswold CE, Hamilton CA, Hedin M, Kocot KM, Ledford JM, Bond JE. 2016. Spider phylogenomics: untangling the spider tree of life. PeerJ 4:e1719

Garwood RJ, Dunlop JA, Selden PA, Spencer AR, Atwood RC, Vo NT, Drakopoulos M. 2016. Almost a spider: a 305-million-year-old fossil arachnid and spider origins. Proceedings of the Royal Society B 283(1827):20160125

Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, Di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A+11 more. 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology 29(7):644-652

Guadanucci JPL. 2014. Theraphosidae phylogeny: relationships of the ‘Ischnocolinae’ genera (Araneae, Mygalomorphae) Zoologica Scripta 43(5):508-518

Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M. 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nature protocols 8(8):1494-1512

Hedges SB. 2001. Afrotheria: Plate tectonics meets genomics. Proceedings of the National Academy of Sciences 98(1):1-2

Heine C, Zoethout J, Müller RD. 2013. Kinematics of the South Atlantic rift. preprint

Hu X, Garzanti E, Wang J, Huang W, An W, Webb A. 2016. The timing of India-Asia collision onset—facts, theories, controversies. Earth-Science Reviews 160:264-299

Joshi J, Karanth PK, Edgecombe GD. 2020. The out-of-India hypothesis: evidence from an ancient centipede genus, Rhysida (Chilopoda: Scolopendromorpha) from the Oriental Region, and systematics of Indian species. Zoological Journal of the Linnean Society 189(3):828-861

Kaderka R, Bulantová J, Heneberg P, Řezáč M. 2019. Urticating setae of tarantulas (Araneae: Theraphosidae): morphology, revision of typology and terminology and implications for taxonomy. PLOS ONE 14(11):e0224384–e0224384

Karanth PK. 2006. Out-of-India Gondwanan origin of some tropical Asian biota. Current Science 90(6):789-792

Katoh K, Kuma K-I, Toh H, Miyata T. 2005. MAFFT version 5: improvement in accuracy of multiple sequence alignment. Nucleic Acids Research 33(2):511-518

Klaus KV, Matzke NJ. 2019. Statistical comparison of trait-dependent biogeographical models indicates that podocarpaceae dispersal is influenced by both seed cone traits and geographical distance. Systematic Biology 69(1):61-75

Kück P. 2009. ALICUT: a Perlscript which cuts ALISCORE identified RSS. Version 2. Bonn (Germany): Department of Bioinformatics, Zoologisches Forschungsmuseum A. Koenig (ZMFK) [cited 22nd June 2018]

Kück P, Meusemann K. 2010. FASconCAT: convenient handling of data matrices. Molecular Phylogenetics and Evolution 56(3):1115-1118

Kück P, Meusemann K, Dambach J, Thormann B, Von Reumont BM, Wägele JW, Misof B. 2010. Parametric and non-parametric masking of randomness in sequence alignments can be improved and leads to better resolved trees. Frontiers in Zoology 7(1):10

Landis MJ, Matzke NJ, Moore BR, Huelsenbeck J. 2013. Bayesian analysis of biogeography when the number of areas is large. Systematic Biology 62(6):789-804

Lartillot N, Lepage T, Blanquart S. 2009. PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating. Bioinformatics 25(17):2286-2288

Lepage T, Bryant D, Philippe H, Lartillot N. 2007. A general comparison of relaxed molecular clock models. Molecular Biology and Evolution 24(12):2669-2680

Li F, Shao L, Li S. 2020. Tropical niche conservatism explains the Eocene migration from India to Southeast Asia in ochyroceratid spiders. Systematic Biology 69(5):987-998

Lüddecke T, Krehenwinkel H, Canning G, Glaw F, Longhorn SJ, Tänzler R, Wendt I, Vences M. 2018. Discovering the silk road: Nuclear and mitochondrial sequence data resolve the phylogenetic relationships among theraphosid spider subfamilies. Molecular Phylogenetics and Evolution 119:63-70

Lüddecke T, Vilcinskas A, Lemke S. 2019. Phylogeny-guided selection of priority groups for venom bioprospecting: harvesting toxin sequences in tarantulas as a case study. Toxins 11(9):488

Marshall LG, Webb SD, Sepkoski JJ, Raup DM. 1982. Mammalian evolution and the great American interchange. Science 215(4538):1351-1357

Matthews KJ, Maloney KT, Zahirovic S, Williams SE, Seton M, Mueller RD. 2016. Global plate boundary evolution and kinematics since the late Paleozoic. Global and Planetary Change 146:226-250

Matzke NJ. 2013. BioGeoBEARS: BioGeography with Bayesian (and likelihood) evolutionary analysis in R Scripts. R package, version 0.2, 1 software

McDonald-Spicer C, Knerr NJ, Encinas-Viso F, Schmidt-Lebuhn AN. 2019. Big data for a large clade: bioregionalization and ancestral range estimation in the daisy family (Asteraceae) Journal of Biogeography 46(2):255-267

McGlone MS, Heenan P, Millar T, Cieraad E. 2018. The biogeography, origin and characteristics of the vascular plant flora and vegetation of the New Zealand mountains. In: Hoorn C, Perrigo A, Antonelli A, eds. Mountains, Climate and Biodiversity. Hoboken, New Jersey: Wiley-Blackwell. 375-389

McKenna MC. 1973. Sweepstakes, filters, corridors, Noah’s arks, and beached Viking funeral ships in palaeogeography. In: Tarling DH, Runcors SK, eds. Implications of continental drift to the earth sciences. New York: Academic Press. 21-46

Misof B, Misof K. 2009. A Monte Carlo approach successfully identifies randomness in multiple sequence alignments: a more objective means of data exclusion. Systematic Biology 58(1):21-34

Morley R. 2018. Assembly and division of the South and South-East Asian flora in relation to tectonics and climate change. Journal of Tropical Ecology 34(4):209-234

Morley R, Plath M, Zhang Y, Li JT. 2016. Biotic interchange between the Indian subcontinent and mainland Asia through time.

O’Dea A, Lessios HA, Coates AG, Eytan RI, Restrepo-Moreno SA, Cione AL, Collins LS, Queiroz ADe, Farris DW, Norris RD. 2016. Formation of the Isthmus of Panama. Science Advances 2(8):e1600883

Oliveira IDS, Bai M, Jahn H, Gross V, Martin C, Hammel JU, Zhang W, Mayer G. 2016. Earliest onychophoran in amber reveals Gondwanan migration patterns. Current Biology 26(19):2594-2601

Olyphant JR, Johnson RA, Hughes AN. 2017. Evolution of the Southern Guinea Plateau: implications on Guinea-Demerara Plateau formation using insights from seismic, subsidence, and gravity data. Tectonophysics 717:358-371

Opatova V, Hamilton CA, Hedin M, De Oca LM, Kral J, Bond JE. 2020. Phylogenetic systematics and evolution of the spider infraorder Mygalomorphae using genomic scale data. Systematic Biology 69(4):671-707

Paradis E. 2013. Molecular dating of phylogenies by likelihood methods: a comparison of models and a new information criterion. Molecular Phylogenetics and Evolution 67(2):436-444

Pérez-Miles F, Costa FG, Toscano-Gadea C, Mignone A. 2005. Ecology and behaviour of the ‘road tarantulas’ Eupalaestrus weijenberghi and Acanthoscurria suina (Araneae, Theraphosidae) from Uruguay. Journal of Natural History 39(6):483-498

Pérez-Miles F, Postiglioni R, Montes-de Oca L, Baruffaldi L, Costa FG. 2007. Mating system in the tarantula spider Eupalaestrus weijenberghi (Thorell, 1894): evidences of monandry and polygyny. Zoology 110(4):253-260

Philippe H, Brinkmann H, Lavrov DV, Littlewood DTJ, Manuel M, Wörheide G, Baurain D. 2011. Resolving difficult phylogenetic questions: why more sequences are not enough. PLOS Biology 9(3):e1000602

Poinar Jr G. 2018. Burmese amber: evidence of Gondwanan origin and Cretaceous dispersion. Historical Biology 31(10):1304-1309

Prentice TR. 1992. A new species of North American tarantula, Aphonopelma paloma (Araneae, Mygalomorphae, Theraphosidae) Journal of Arachnology 20:189-199

Raven RJ. 1985. The spider infraorder Mygalomorphae (Araneae): cladistics and systematics. Bulletin of the AMNH 182 (article 1)

Ree RH, Sanmartín I. 2018. Conceptual and statistical problems with the DEC+J model of founder-event speciation and its comparison with DEC via model selection. Journal of Biogeography 45(4):741-749

Ree RH, Smith SA. 2008. Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Systematic Biology 57(1):4-14

Reichling SB. 2000. Group dispersal in juvenile Brachypelma vagans (Araneae, Theraphosidae) The Journal of Arachnology 28(2):248-251

Revell LJ. 2012. phytools: an R package for phylogenetic comparative biology (and other things) Methods in Ecology and Evolution 3(2):217-223

Rice P, Longden I, Bleasby A. 2000. EMBOSS: the European molecular biology open software suite. Elsevier current trends. software

Robinson DF, Foulds LR. 1981. Comparison of phylogenetic trees. Mathematical Biosciences 53(1–2):131-147

Santana R, Perez D, Dobson J, Panagides N, Raven R, Nouwens A, Jones A, King G, Fry B. 2017. Venom profiling of a population of the theraphosid spider Phlogius crassipes reveals continuous ontogenetic changes from juveniles through adulthood. Toxins 9(4):116

Schenk JJ. 2016. Consequences of secondary calibrations on divergence time estimates. PLOS ONE 11(1):e0148228

Schmidt G. 2003. Die Vogelspinnen: Eine weltweite Übersicht. Westarp-Wissenschaften. Hohenwarsleben, Germany: Neue Brehm-Bücherei. Vol. 641

Schultz SA, Schultz MJ. 2009. The Tarantula keeper’s guide (3rd Edition). New York: Barron’s Educational Series.

Selden PA. 2002. First British Mesozoic spider, from Cretaceous amber of the Isle of Wight, southern England. Palaeontology 45(5):973-983

Selden PA, Gall J-C. 1992. A Triassic mygalomorph spider from the northern Vosges, France. Palaeontology 35:211-235

Selden PA, Shcherbakov DE, Dunlop JA, Eskov KY. 2014. Arachnids from the Carboniferous of Russia and Ukraine, and the Permian of Kazakhstan. Paläontologische Zeitschrift 88(3):297-307

Sharma PP, Baker CM, Cosgrove JG, Johnson JE, Oberski JT, Raven RJ, Harvey MS, Boyer SL, Giribet G. 2018. A revised dated phylogeny of scorpions: phylogenomic support for ancient divergence of the temperate Gondwanan family Bothriuridae. Molecular Phylogenetics and Evolution 122:37-45

Shillington C, McEwen B. 2006. Activity of juvenile tarantulas in and around the maternal burrow. The Journal of Arachnology 34(1):261-266

Simms D, Cizdziel PE, Chomczynski P. 1993. TRIzol: a new reagent for optimal single-step isolation of RNA. Focus 15(4):532-535

Springer MS, Meredith RW, Janecka JE, Murphy WJ. 2011. The historical biogeography of Mammalia. Philosophical Transactions of the Royal Society B 366(1577):2478-2502

Stamatakis A. 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30(9):1312-1313

Starrett J, Hedin M. 2007. Multilocus genealogies reveal multiple cryptic species and biogeographical complexity in the California turret spider Antrodiaetus riversi (Mygalomorphae, Antrodiaetidae) Molecular Ecology 16(3):583-604

Udvardy MD. 1975. A classification of the biogeographical provinces of the world (Vol. 8). Morges, Switzerland: International Union for Conservation of Nature and Natural Resources.

West RC, Nunn SC, Hogg S. 2012. A new tarantula genus, Psednocnemis, from West Malaysia (Araneae: Theraphosidae), with cladistic analyses and biogeography of Selenocosmiinae Simon 1889. Zootaxa 3299(1):1-43

Westerweel J, Roperch P, Licht A, Dupont-Nivet G, Win Z, Poblete F, Ruffet G, Swe HH, Thi MK, Aung DW. 2019. Burma Terrane part of the Trans-Tethyan arc during collision with India according to palaeomagnetic data. Nature Geoscience 12:863-868

Wunderlich J, Müller P. 2020. New and already described fossil spiders (Araneae) of 20 families in Mid Cretaceous Burmese amber with notes on spider phylogeny, evolution and classification. Beitr. Araneol. 13:22-164

Xia X, Xie Z. 2001. DAMBE: software package for data analysis in molecular biology and evolution. Journal of Heredity 92(4):371-373

Xia X, Xie Z, Salemi M, Chen L, Wang Y. 2003. An index of substitution saturation and its application. Molecular Phylogenetics and Evolution 26(1):1-7

Yáñez M, Floater G. 2000. Spatial distribution and habitat preference of the endangered tarantula, Brachypelma klaasi (Araneae: Theraphosidae) in Mexico. Biodiversity & Conservation 9(6):795-810

Yu Y, Harris AJ, Blair C, He X. 2015. RASP (Reconstruct Ancestral State in Phylogenies): a tool for historical biogeography. Molecular Phylogenetics and Evolution 87:46-49

Zheng Y, Wu F. 2018. The timing of continental collision between Indian and Asia. Science Bulletin 63(24):1649-1654

Zhou B, Tu T, Kong F, Wen J, Xu X. 2018. Revised phylogeny and historical biogeography of the cosmopolitan aquatic plant genus Typha (Typhaceae) Scientific Reports 8(1):8813