Scincidae

Phylogenetic analyses The total concatenated dataset consisted of 4,758 base pairs (best fit model = GTR + F + R4). Best fit models for each partition are as follows: ND4 (755bp, TIM2+F+R4), 16S (480bp, TIM2+F+I+G4), COI (715bp, TPM2u+F+R3), RAG-1 (810bp, TPM3u+F+I+G4), β-Globin-2 (853bp, TVM+F+G4), C-mos (400bp, TN+F+G4), ALD (235bp, SYM+R2), MHC-2 (510bp, TN+F+G4), mitochondrial (1950 bp, GTR+F+R4), and nuclear (2,808bp, TVMe+R3). We find no well supported conflict between the partitioned and concatenated results (see Supplementary Figure 1), therefore results from the concatenated analysis are presented here. The concatenated Bayesian analysis (Fig. 1) strongly supports (1.0pp) the monophyly of all three genera of rainbow skinks, with Liburnascincus sister to Lygisaurus + Carlia (.97pp), as found in Dolman & Hugall (2008) but in contrast to Bragg et al. (2018). Within Lygisaurus our results support two distinct clades, one composed of Australian species, and the other with the Australian L. macfarlani, (referred to as L. novaeguineae Meyer, 1874, in Dolman & Hugall 2008; changed here based on locality) sister to all New Guinean samples. The placement of L. sesbrauna Ingram & Covacevich, 1988, is problematic and unresolved in the concatenated tree. Lygisaurus curtus, is monophyletic (.94pp) with long branch lengths and deep genetic distances across localities. Within Carlia, the fusca group is sister to the bicarinata group, with all remaining Australian taxa sister this clade (1.0pp). Support for the monophyly of the C. bicarinata clade is low (.77pp) and results in a basal tritomy, due to the placement of C. schmeltzii. Our results add to the support for the monophyly of the C. fusca group (1.0 pp) (Dolman & Hugall 2008; Donnellan et al. 2009; Stuart-Fox et al. 2002; Zug 2004) and divide it into four well supported clades, with the exception of C. caesius Zug and Allison, 2006 where placement is poorly supported (Fig. 2). Estimation of Divergence Times and biogeographic history We find the DEC+J model to be a significantly better fit to our phylogeny compared to the DEC model (AIC 268.4 v.s 333.7). There are problems with the use of the DEC+J model (Ree & Sanmartín 2018) with models returning zero dispersal and extinction probabilities. We find similar problems in our biogeographical analysis with d=0 (dispersal) and e=0 (extinction). Although this model seems appropriate for testing biogeographical hypotheses of Carlia considering their widespread distribution across the Islands of Oceania (Zug 2004; Zug & Allison 2010; Donnellan et al. 2009), we recommend caution in the interpretation of the results. We find Australian Carlia and the fusca + bicarinata group diverge ~16.9 mya (Fig. 3, for confidence intervals on dates see Supplementary Figure S2), followed by the rapid radiation of the fusca group beginning ~10.4 mya in Australia. Our data support four independent colonizations of New Guinea from Australia; single invasions by Lygisaurus and the C. bicarinata group, and two by the C. fusca group. New Guinean Lygisaurus diverged ~8.7 mya, dispersing to the Craton and then the EPCT, while the C. bicarinata group dispersed directly to the EPCT (~7.0mya). The first invasion of the C. fusca group to New Guinea occurred ~9.8 mya into either the Craton or the EPCT (the relationships between C. caesius from the Craton is poorly resolved in our phylogenetic analysis (.56, see Fig. 2) and the probability of colonization through the Craton is only slightly higher than the EPCT). The EPCT appears to have served as a “launching point” for dispersal to the fold belt and Accreted terranes as well as island dispersal and speciation, with Darnley, Fergusson, Woodlark, and Sudest populations all originating from the EPCT. The second fusca group colonization event involves the dispersal of the Island clade to the VogelKop peninsula (~4.7mya), with additional long-distance dispersal events to Halmahera, New Ireland, Palau, and New Britain. Morphometric analyses All morphological measurements are reported in supplementary data (Table S 3). MANOVA tests revealed significant sexual dimorphism in New Britain Carlia (p = 0.01), but not in Lakekamu (p=.174), as Welch’s T-test showed significant differences in HL (p=0.00065, t= 4.20), HW (p=0.0076, t = 3.014), and PlD (p=0.0039, t= 3.50). To remove the effects of sexual dimorphism in interpopulation comparisons, these characters (HL, HW, and PlD) were removed from subsequent analyses. The Broken stick method revealed only PC1 (pectoral width, pelvic width) and PC6 (pelvic depth, pelvic width) (Table 1) of our Principal Components Analysis to explain a substantial portion of variance, explaining 41.5% of the cumulative variance. The morphospace resulting from PCA reveals considerable overlap between regions (Fig. 4). Results of Tukey Honesty significant difference tests of ANOVA results are shown in Table 2, with significant differences marked with an asterisk. New Britain and Maimafu were found to contain no significant differences in characters analyzed. No more than 1–2 significant differences were found between other regions apart from New Britain and Lakekamu, with Lakekamu specimens possessing wider bodies and longer forelimbs. DNA Barcode COI haplotype networks were able to differentiate between recognized species and identify multiple clades that may represent new species (see Fig. 5; for simplicity only the network for Lygisaurus is shown, see Supplementary Figure S2 for the Carlia haplotype network).