Comparative Cranial Morphology of Squamates

Analyses were largely conducted in R v3.4.1 (75 ). First, a disparity-through-time plot was created for the individual partitions by using the dtt function in the GEIGER R package v2.0.6 (76 (link), 77 (link)). Second, we used BayesTraitsV3 (78 ) to estimate evolutionary rates for the overall skull and its partitions based on PCs accounting for >95% of total shape variation (SI Appendix, Table S3). All available evolutionary models were tested, including BM, OU, δ, κ, and λ models. Bayes Factor was used to determine the best-supported evolutionary model. Estimated rates based on the preferred model were used to generate a rate-through-time plot for each partition. Total disparity (Procrustes variance) and mean evolutionary rates under the BM model were calculated within each model and for each landmark and sliding semilandmark by using the morphol.disparity and compare.multi.evol.rates functions in the geomorph R package v.3.0.5 (79 , 80 ), respectively. We simulated shape evolution under BM model by using the sim.char function in GEIGER R package. Ancestral shapes were calculated by using the anc.recon function in the Rphylopars R package v0.2.9 (81 (link)) on extant-only and combined datasets.
The geomorph R package was used to assess allometry, phylogenetic signal, evolutionary allometry, and phylogenetically corrected correlations with ecological and reproductive variables in partition-specific and whole skull-shape data by using the functions procD.allometry, physignal, and procD.pgls (SI Appendix, Tables S6, S9, and S11). Because of the proportionately low contribution of (evolutionary) allometric signal, we did not correct for allometry in the datasets for other analyses. To identify cranial modules from partitions, we employed EMMLi (26 (link)), a maximum-likelihood approach, and CR analysis (27 ), which is implemented in geomorph. By using the results output by EMMLi, we further assessed the pattern of integration by comparing estimated within-module correlations (ρ_within) with the estimated between-module correlation (ρ_between) for each pair of partitions (SI Appendix, SI Text). We applied a subsampling approach to ascertain the robustness of these results by running 100 repetitions of EMMLi with 10% of the total landmark dataset while keeping a minimum of three landmarks + sliding semilandmarks per partition, which consistently returned the same pattern of cranial modules, as did analyses limited to landmarks only (SI Appendix, Fig. S9). Results from EMMLi and CR analyses with and without phylogenetic correction yielded largely congruent results (SI Appendix, Tables S7 and S8). Least-squares linear regression was conducted on these values to determine any relationships between disparity, mean evolutionary rates, and within-module correlation (SI Appendix, Table S9). Integration and macroevolutionary patterns were assessed separately for extant lizards, extant snakes, all extant squamates, and combined extant and extinct squamates.