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

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.

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.

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.

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.

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).

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.