Arrival and diversification of mabuyine skinks (Squamata: Scincidae) in the Neotropics based on a fossil-calibrated timetree

Introduction

The subfamily Mabuyinae comprises 24 genera and 197 species of lizards and belongs to the highly diverse worldwide family, Scincidae (Lepidosauria; Squamata) (Uetz, Hošek & Hallermann, 2017). The Mabuyinae, as well as the Scincidae, are distributed on nearly all continents. Approximately one-third of the species in this subfamily are from the Neotropical region, and they are the only representatives of scincids in South America.

The entire diversity of mabuyine species was traditionally assigned to the single genus, Mabuya, but a previous analysis proposed four new monophyletic genera with well-defined geographical distributions (Mausfeld et al., 2002): (1) Trachylepis (previously Euprepis), comprising African and Madagascan species with one South American representative (T. atlantica); (2) Eutropis, containing Asian species; (3) Chioninia, from the Cape Verde islands; and (4) Mabuya, containing New World species, recently rearranged into 16 new genera by Hedges & Conn (2012). Hedges & Conn (2012) treated South American Mabuyinae as a clade of family Mabuyidae, within the superfamily Lygosomoidea. The new Mabuyidae consisted of four new subfamilies: Mabuyinae, Chioniniinae, Dasiinae and Trachylepidinae. Hedges and Conn’s arrangement disregarded many genera previously related to Mabuyinae, such as Eutropis, Lankaskincus, and Ristella (Pyron, Burbrink & Wiens, 2013), and recent papers have questioned this classification (e.g., Pyron, Burbrink & Wiens, 2013; Pinto-Sanchez et al., 2015; Karin et al., 2016). Recently, alternative classifications were suggested. For instance, Karin et al. (2016) placed the Middle-Eastern Trachylepis (T. aurata, T. vittatus, and T. septemtaeniatus) into the genus Heremites and Eutropis novemcarinata into Toenayar novemcarinata. Moreover, Pinto-Sanchez et al. (2015) assigned species from the genera Maracaiba and Alinea back to genus Mabuya and Metallinou et al. (2016) reclassified Trachylepis ivensii as Lubuya ivensii. In this study, we followed the Reptile Database (Uetz, Hošek & Hallermann, 2017) taxonomy as of January 2017, which classified this clade as subfamily Mabuyinae of the family Scincidae.

Therefore, at least two phylogenetically distinct lineages of Mabuyinae are distributed in the New World, namely, T. atlantica and the Continental American Mabuyinae (CAM) clade. These lineages are distinguished by both morphological features—presacral vertebrae counts, keeled dorsal scales, coloration, and oviparity (Greer, 1970)—and molecular evidence (e.g., Mausfeld & Vrcibradic, 2002; Carranza & Arnold, 2003; Whiting et al., 2006). It is customary to assume that the history of Mabuyinae in South American continent consisted of two independent transoceanic dispersal events from the Old World (Mausfeld & Vrcibradic, 2002; Carranza & Arnold, 2003; Whiting et al., 2006). The spatial distribution of the single representative of the genus Trachylepis in an island closer to South America than to Africa, T. atlantica, is of particular interest. This issue is so intriguing that it has reached the pages of nontechnical literature (De Queiroz, 2014). T. atlantica is found in Fernando de Noronha, a small volcanic archipelago in the Atlantic Ocean that lies 375 km off the northeastern coast of Brazil and that was geologically formed from the Miocene (12.3 Ma) to the upper Pliocene, from 3.3 to 1.7 Ma (Almeida, 2002). Although other Trachylepis species have spread to several islands near the African continent, the presence of T. atlantica in Fernando de Noronha likely represents the farthest dispersal registered for the genus.

The CAM lineage, on the other hand, contains approximately 60 species of Mabuyinae. Previous studies have suggested that the split between this clade and the African Mabuyinae (genus Trachylepis) occurred from 28 to 34 Ma (Hedges & Conn, 2012; Karin et al., 2016). Additionally, the age of the last common ancestor (LCA) of the CAM was dated at 7–9 Ma (Carranza & Arnold, 2003; Pinto-Sanchez et al., 2015) and 11–14 Ma (Miralles & Carranza, 2010; Hedges & Conn, 2012; Karin et al., 2016). Considering the present-day Atlantic Ocean currents, Mausfeld et al. (2002) suggested that T. atlantica could have dispersed from the coast of Southwest Africa to South America, but no work to date has comprehensively evaluated this hypothesis. Older ages could be consistent with the supposed faunal transoceanic connection between Africa and South America during the Eocene/Oligocene (e.g., hystricognath rodents, Loss-Oliveira, Aguiar & Schrago, 2012; Voloch et al., 2013; anthropoid primates, (Schrago et al., 2012; Schrago, Mello & Soares, 2013); amphisbaenians, Vidal et al., 2008; emballonurid bats, Teeling, 2005; Leigh, O’Dea & Vermeij, 2013; testudinid turtles, Le et al., 2006) through a single or a series of islands constituting an Atlantic Ocean Ridge (the Atlantogea paleo province) (Simpson, 1950; Poux et al., 2006; De Oliveira, Molina & Marroig, 2009; Ezcurra & Agnolin, 2012; De Queiroz, 2014). This transatlantic island corridor, associated with a drop in the sea level in the Oligocene, could explain faunal exchanges in this period (De Queiroz, 2014). Regarding T. atlantica, there is a lack of estimates of the age of the separation between this species and its African sister lineage.

Chronological information is indispensable for a full understanding of the scenario underlying the current geographic distribution and evolutionary history of extant lineages (Sanmartín, Van der Mark & Ronquist, 2008; Loss-Oliveira, Aguiar & Schrago, 2012). However, estimates of divergence times within the Mabuyinae have been hampered by the lack of fossils for Mabuyinae, which made previous researchers rely on evolutionary rates borrowed from the literature (e.g., Carranza & Arnold, 2003; Miralles & Carranza, 2010; Lima et al., 2013; Barker et al., 2015; Karin et al., 2016) or to employ biogeographic events as calibrations (Hedges & Conn, 2012; Pinto-Sanchez et al., 2015) to derive the timescale of this lineage. The fossil record, however, has been demonstrated to be much more informative than biogeographic events as a calibration tool (Heads, 2011). The use of biogeographic events to time-calibrate phylogenies requires the assumption of vicariant scenarios of diversification, which entails that the speciation is synchronous with the breaking apart of landmasses, or that island colonization occurs immediately following geological formation (Heads, 2011; Mello & Schrago, 2012).

Revealing the origin and diversification of the CAM clade and T. atlantica also requires a robust phylogenetic hypothesis. Although it is generally accepted that T. atlantica is a member of the genus Trachylepis and is therefore excluded from the main diversification of the CAM, its evolutionary affinity remains controversial (Mausfeld et al., 2002; Carranza & Arnold, 2003; Whiting et al., 2006). Early proposed phylogenies of Mabuyinae used a maximum of 35 species, a number that significantly underrepresents the diversity of this subfamily (Mausfeld et al., 2002; Carranza & Arnold, 2003; Whiting et al., 2006). Hedges & Conn (2012) studied the CAM clade exclusively and included 40 species, whereas the large-scale Squamata phylogeny of Pyron, Burbrink & Wiens (2013), which contains 4,161 species, sampled 71 mabuyine skink species. Recently, Pinto-Sanchez et al. (2015) used 250 specimens to infer the phylogeny and species diversity of neotropical mabuyinaes, focusing on Colombian populations. Karin et al. (2016) analyzed the higher-order relationships of Mabuyinae using 22 species (24 specimens). Considering that improvements in phylogenetic inference and divergence time estimation can be obtained by increasing taxon sampling (Linder, Hardy & Rutschmann, 2005; Albert et al., 2009; Soares & Schrago, 2012), this matter requires further investigation.

Moreover, the occupation of the South American mainland by the CAM motivates an analysis of differential rates of diversification and rates of phenotypic traits evolution in this clade. It has been broadly reported that the ecological opportunity of a new environment can induce acceleration of macroevolutionary rates (Yoder et al., 2010; Liedtke et al., 2016). This phenomenon was documented for insular vertebrates (Losos & Mahler, 2010; Jønsson et al., 2012); but, as proposed by Simpson (1953), the occupation of a new continent could also trigger evolutionary radiations (Yoder et al., 2010; Pires, Silvestro & Quental, 2015).

In this study, as a means of exploring the continental biogeographic and evolutionary patterns associated with the occupation of South America by mabuyine skinks, we estimated the molecular phylogeny and inferred the divergence times of its members employing fossil information. The inferred time-tree was used to perform, for the first time, a formal statistical analysis of the historical biogeography and macroevolutionary diversification rates of Mabuyinae. To this end, by combining previously published data, we assembled the largest dataset of species and gene sampling composed to date, with the aim of uncovering the evolution of the Mabuyinae.

Materials and Methods

Results

Mabuyinae was recovered as a monophyletic group (BS = 76). The first split in this subfamily isolated the genus Dasia from the remaining Mabuyinae and was inferred to have occurred at 30 Ma, with an HPD interval ranging from 20 to 48 Ma (File S4). The biogeographical model with the highest likelihood was the DEC + J (lnL: −94.34). Therefore, the addition of the J parameter for founder events significantly increased the likelihood of the DEC model (p = 1.2e⁻⁹, Table 3). Our results supported a model in which the genera Lankascincus and Ristella split from Mabuyinae (BS = 20) in the Oriental region.

After dispersal from the Oriental region to Africa, the ancestor of the CAM clade split from Saharo-Arabian Trachylepis approximately 25 Ma (±16–41 Ma) (BS = 60). Subsequently, the stock that gave rise to New World Mabuyinae dispersed from the Saharo-Arabian region to the Neotropical Brazilian Subregion (Fig. 1). The New World Mabuyinae were recovered as monophyletic (BS = 97), and the age of their LCA was estimated as 19 Ma (12–31 Ma). Our results supported the hypothesis that this lineage reached the Neotropics via the Brazilian Subregion. From this region, the ancestor of Spondylurus and Copeoglossum dispersed to islands of the West Indies. The continental clade later dispersed along the Brazilian Subregion. Additional dispersal events to the West Indies were inferred to have occurred at least twice: in the ancestor of Mabuya and in that of Alinea and Marisora. Marisora was not recovered as monophyletic, with the Central American clade as sister to Alinea and the South American clade as sister to Aspronema.

The origin of Trachylepis atlantica consisted of a different history. Firstly, the monophyly of genus Trachylepis was not recovered. This genus was split into two lineages: one sister lineage of the CAM clade and another major clade to which T. atlantica belongs. Our analysis supported the phylogenetic position of T. atlantica as a sister lineage of the clade including the Europa Island T. infralineata, the Madagascan T. comorensis, and the continental African T. maculilabris (BS = 29). We found that the T. atlantica ancestor diverged from the remaining Trachylepis in the Miocene, approximately 17 Ma (between 10 and 27 Ma) in tropical Africa. The younger diversification of T. atlantica indicates that the crossing of this species to Fernando de Noronha occurred more recently than did the occupation of South America by the CAM clade. Our analyses also suggested that once the ancestors of this major Trachylepis clade reached Africa, dispersal events from the African continent occurred to Madagascar and nearby islands at least three times, in addition to the dispersal to the Neotropics (T. atlantica). Our results showed that a monophyletic Madagascan Trachylepis (BS = 94) diverged from its sister group approximately 19 Ma (12–32 Ma) (BS = 23).

The ANOVA was used to test whether the full four-parameter GeoSSE model of the $NW∕\overline{NW}$ comparison (model 1) was statistically favored over simpler alternative models (models 2–4). When we constrained speciation rates to be equal (model 2) and to present equal speciation and extinction rates (model 4), the full model was supported over these alternative models (p = 0.00782 and p = 0.02900, respectively). The comparison with the model that constrained extinction rates (model = 3), however, was not significantly different (p = 0.15575).

GeoSSE analysis of the difference between the posterior distributions of speciation rates in $\overline{NW}∕\overline{NW}$ lineages indicated a positive credible interval (CI) (0.005–0.113), suggesting that the speciation rates of South American clades (NW) were indeed higher than those of the remaining lineages ($\overline{NW}$). In BAMM, however, we failed to find a significant difference between the posterior distributions of speciation rates (CI: −0.030–0.031). We also found no evidence for diversification shifts during mabuyine diversification. However, the macroevolutionary cohort matrix suggested heterogeneous macroevolutionary rate regimes in the South American Mabuyinae and the remaining clades (Fig. 2A). In both GeoSSE and BAMM analyses, we failed to find differences between the extinction rates (CI_geoSSE: −0.036–0.074, CI_BAMM: −0.027–0.027) and between the net diversification rates (CI_geoSSE: −0.010–0.090, CI_BAMM: −0.027–0.027) (Fig. 2B).

Concerning the influence of viviparity on diversification rates, no submodels were supported over the unconstrained model (p > 0.05). STRAPP analysis failed to find any correlation between diversification rates and mode of reproduction (p > 0.05). We also found no differences between the means obtained from the MCMC of the BiSSE (CI: −0.018–0.125, −0.038–0.055, −0.051–0.049, respectively for speciation, extinction, and transition rates).

The ancestor of CAM and their Trachylepis sister clade was recovered as viviparous by BayesTraits (P_BTML = 0.751; P_BTBA = 0.787) and as oviparous by BiSSE (P_BiSSE = 1), while the ancestor of CAM alone was recovered as viviparous by all approaches (P_BiSSE = 0.953; P_BTML = 0.992; P_BTBA = 0.983), suggesting that the ancestral species that colonised South America was likely viviparous (Fig. 2A). The ancestor of the subfamily Mabuyinae was inferred as oviparous (P_BiSSE = 1; P_BTML = 0.983; P_BTBA = 0.980). Viviparity seems to have emerged several times in the genus Trachylepis, as well as in the ancestor of Chioninia (P_BiSSE = 0.952; P_BTML = 0.984; P_BTBA = 0.971).

According to PIC results, SVL measurement and body mass were correlated (p_LRP < 0.01), however, we had found differences in their evolutionary rates. The mean rates of trait evolution of SVL and body mass were smaller in the New World clade than in the $\overline{NW}$ lineages. However, these differences were not significant. In BAMM, we found an increase in the rate of SVL trait evolution in genus Chioninia. Only one clade presented an increase in the rate of body mass trait evolution: the Seychelles Trachylepis (T. wrightii and T. brevicollis). STRAPP analysis failed to find a significant correlation between the rate of species diversification and the rate of SVL or body mass trait evolution (p > 0.05).

Discussion

Our analyses based on a fossil-calibrated timetree support the hypothesis that the occupation of South America by the ancestors of the Continental America Mabuyinae (CAM) clade occurred between the Eocene and the Oligocene from the Saharo-Arabian region, corroborating recent estimates (e.g., Pinto-Sanchez et al., 2015; Karin et al., 2016). This colonization occurred independently of that of Trachylepis atlantica, which split from its tropical African sister clade between the Oligocene and the Miocene. Therefore, our results corroborate the hypothesis that the ancestor of T. atlantica dispersed from tropical Africa as suggested by Mausfeld et al. (2002). The divergence between T. atlantica and its African ancestors was older than the estimated age of the formation of the Fernando de Noronha archipelago (12.3–1.7 Ma; Almeida, 2002). Thus, it is unlikely that the ancestors of T. atlantica reached South America in a single transoceanic dispersal event. It is more probable to hypothesize a complex dispersal scenario, with a series of dispersals through a transatlantic island corridor leading to Fernando de Noronha, consistent with the Atlantogea model (De Oliveira, Molina & Marroig, 2009; Ezcurra & Agnolin, 2012). This stepping-stone mode of dispersal could explain the discrepancy between the ages of genetic divergence and the formation of the archipelago because it would dissociate the geological formation of the Fernando de Noronha archipelago from the genetic isolation of the ancestor of T. atlantica. This result must be interpreted cautiously because choosing between alternative biogeographical scenarios is influenced by the inferred timescale. The lack of extant and extinct mabuyine species may have impacted divergence time estimates.

The age of the crown node of living CAM species was inferred to be much older than the estimates reported by Carranza & Arnold (2003), 7–9 Ma. The credibility interval of our estimate (12–31 Ma), however, contains the ages inferred by Miralles & Carranza (2010), Hedges & Conn (2012) and Karin et al. (2016), which dated this node at approximately 14 Ma. Our results were still inconclusive to establish the precise sequence of events that gave rise to the current geographical distribution of both CAM and T. atlantica. Nevertheless, the evolution of these two American lineages of Mabuyinae corroborates a faunal connection between Africa and the Neotropics, and suggests that dispersals may have occurred through island hopping, which is more likely than a single sweepstake event across the Atlantic Ocean. In this sense, it is interesting to mention that a species has been described in the mid-Atlantic Ascension Island by Gray in 1838, which was subsequently assigned to genus Mabuya; however, no recent systematic re-evaluation of this specimen is available (Mausfeld & Vrcibradic, 2002).

The arrival of Mabuyinae to South America seems to have changed the diversification rate regime of this lineage. The evolutionary rates of both body mass and SVL traits were not significantly altered after the arrival of the ancestors of New World clades, but the increases in rates of both traits were found in island groups, such in Cape Verde’s Chioninia and Seychelles’ Trachylepis. The viviparous CAM clade presents peculiar reproductive traits, such as a specialized chorioallantoic placenta that provides fetal nutrition, similarly to eutherian mammals (Vrcibradic & Rocha, 2011). This character was suggested to be synapomorphous by Mausfeld et al. (2002). According to our analysis, the CAM ancestor was already viviparous, in agreement with the large-scale analysis of Squamates (Pyron & Burbrink, 2013). However, Pyron & Burbrink (2013) suggested an early origin of viviparity in Mabuyinae and multiple transitions to oviparity, including the clade consisting of Chioninia and CAM + Heremites. On the other hand, possibly because of topological differences, we estimated that the ancestor of Mabuyinae was oviparous and multiple reversions to viviparity occurred later along the evolution of the lineage.

Regarding the biogeographic history of the subfamily Mabuyinae, Greer (1970), based on the primitive characteristics of Asian Mabuyinae, suggested that the ancestors of this subfamily first dispersed from Asia. Other works reasserted this hypothesis, although no formal statistical analyses were conducted (Honda et al., 1999; Honda et al., 2003; Karin et al., 2016). Our results suggest that the early members of Mabuyinae dispersed to Africa and then reached South America and several oceanic islands. The inferred time-dated phylogeny and biogeographic reconstruction similarly suggest an Oriental ancestral distribution, as corroborated by the distribution of the genera Lankascincus and Ristella, which were recovered as sister clades of Mabuyinae and are geographically distributed in the Oriental region (Sri Lanka and the Indian, respectively). Karin et al. (2016) inferred that Eutropis, instead of Dasia, was the first Mabuyinae offshoot. Nevertheless, as both genera are distributed in the Oriental region, this topological rearrangement would not significantly alter the ancestral area reconstruction.

Finally, this study provides an updated timescale and estimates of macroevolutionary regimes of the diversification of the subfamily Mabuyinae. Our focus was the occupation of the South American continent by the subfamily Mabuyinae through an Oligocene/Eocene transoceanic connection (Atlantogea), which could be responsible for approximately 29% of the South American mammal diversity, suggesting the importance of this event to the extant vertebrate diversity in South America (Marshall et al., 1982). CAM and T. atlantica are the only representatives of the family Scincidae in South America accounting for 4% of the approximately 1,560 squamate species on this continent (Uetz, 2000). Although there are species of other subfamilies of Scincidae in Central America, these species seem to have never crossed the Isthmus of Panama. The converse also seems correct because the only mabuyine genus (Marisora) that partially occupied Central America based on our results had arrived via an intermediate dispersion through the West Indies. If this is the case, no representative of Scincidae would be found in present day South American fauna in the absence of this Atlantogea connection. Our results give an overall picture of the timing, biogeography and macroevolutionary dynamics associated with the arrivals of the ancestors of this exceptional case of transoceanic dispersal in two closely related lineages.

Supplemental Information