Modularity.
The network modules of all species analyzed are shown in Supplementary Figs. 1–58. Despite the morphological diversity of lizard skulls, lizards generally possess separate left and right preorbital (purple and red), postorbital (blue and orange), and mandibular modules (light and dark gray) (Fig. 1). Nevertheless, in some taxa, the snout (light purple), including the premaxilla, nasal and frontal (e.g., Basiliscus vitattus, Draco volans, Tupinambis teguixin), or the braincase elements (yellow) (e.g., Anolis cristatellus, Elgaria panamintina) form a single module. Furthermore, in the gekkotans, the frontals, parietals, and postorbitals form a skull roof module (pink), while in other taxa, the nasals are included in or the parietals are excluded from the skull roof module. Only in chamaeleonids, the parietals are integrated into the preorbital module, while the iguanians with ornamentation similar to chamaeleonids (Phrynosom asio, Basiliscus vitattus) have their parietals integrated into the postorbital module. Rhineura floridana exhibits a unique pattern in which all the cranium bones are integrated into a single module on each side. In two species of geckos, Coleonyx variegatus and Oedura tryoni, the pterygoid forms a separate module with the epipterygoid (see Supplementary Figs. 10 and 32).
The boundaries between modules in the dorsal and ventral regions are not constant, and the modules, including the parietal and pterygoid, differ from species to species. Except for Heloderma horridum, Brookesia brygooi, Rhampholeon brevicaudatus, and Oplurus cyclurus, the lateral boundaries are almost constant in the jugal-postorbital (postorbitofrontal).
The skull of the tuatara differs significantly from lizards and shows the preorbital module containing the temporal bones, a braincase module, and left and right mandibular modules. Interestingly, the jugal of the tuatara is highly integrated with the postorbital (in the dendrogram, the jugal and postorbital are adjacent to each other (Supplementary Fig. 1)), while in lizards, they are in separate modules (Fig. 1; Supplementary Figs. 2–58).
Multivariate analyses of network parameters.
The PC 1 and PC 2 of the network parameters together account for more than 65% of the total variation (see Supplementary data 4 file). The PC 1 explains most of the parameters except for H. Negative PC 1 values relate to greater N, K, L, Q-modules, S-modules, and Qmax, and positive values relate to greater D and C. Amphisbaenia, a fossorial taxon with greater D and less N, K, and L, exhibits larger PC 1 scores. However, within Amphisbaenia, Rhineura floridana (C = 0.5418546) and Bipes biporus (C = 0.4955357) with a greater integration differ from Amphisbaena alba (C = 0.3447368) and Trogonophis wiegmanni (C = 0.3429654), resulting in larger PC 1 score. The basal Lacertoidea, Lacertidae and Teiidae, are intermediate and separated from the derived Lacertoidea, Amphisbaenia, along with PC 1. The PC 1 score for the tuatara was 0.5331 and intermediate.
Negative PC 2 values relate to the greater C and H and positive values relate to greater K. Most gekkotans with specialized skulls without postorbital bars and upper temporal bars are plotted on the negative side of PC 2 due to the small value of K. The Mann-Whitney U test strongly supports that Gekkota (n = 14) and other lizards (n = 44) differ from each other in present multivariate analyses (Table 1). In other words, only N is lower in the skull of Gekkota than in that of other lizards, and the reduction in connectivity due to the absence of postorbital bars and upper temporal bars in Gekkota does not seem to affect other parameters. Notably, the PC 2 score of the tuatara is the greatest (2.6805), which is due to relatively low C (0.3544974) and the lowest H (0.2763419). This result indicates that lizards evolved skulls that were highly integrated and had greater anisomerism than the tuatara.
Table 1. Comparison of network parameters and principal components scores using the Mann-Whitney U test. Values with significant differences are shown in bold.
|
fossorial (n = 6) vs.
non-fossorial (n = 52)
|
Gekkota (n = 14) vs.
non-Gekkota (n = 44)
|
|
z-value
|
p-value
|
z-value
|
p-value
|
N
|
3.30537
|
0.00095
|
1.903872
|
0.05693
|
K
|
2.80996
|
0.00496
|
3.18151
|
0.00147
|
D
|
3.37071
|
0.00075
|
0.890493
|
0.3732
|
C
|
0.178725
|
0.8582
|
1.055153
|
0.2914
|
L
|
2.80854
|
0.00498
|
1.237532
|
0.2159
|
H
|
0.663846
|
0.5068
|
0.981239
|
0.3265
|
S-modules
|
1.707661
|
0.0877
|
0.427791
|
0.6688
|
Q-modules
|
2.84809
|
0.0044
|
1.912208
|
0.05585
|
Qmax
|
3.24258
|
0.00119
|
0.620519
|
0.5349
|
PC 1
|
3.24258
|
0.00119
|
0.964282
|
0.3349
|
PC 2
|
0.204257
|
0.8382
|
2.84763
|
0.00441
|
The pPC 1 and pPC 2 of the network parameters together account for about 60% of the total variation (see Supplementary data 4 file). The distributions of pPC1 and pPC2 are essentially unchanged compared to the PCA results, which indicates that there is not much phylogenetic signal in network parameters (Fig. 2). However, the plot for most phylogenetically basal and fossorial species, Dibamus novaeguineae, apparently shifts its placement compared to the PCA results, becoming more similar to phylogenetically distant and alike fossorial Amphisbaenia (Fig. 2a, g).
In each ecological category, the network parameters did not differ by diet, but they by habitats and locomotion. Analyses on habitats and locomotion (Fig. 2c, d) result in greater PC 1 scoring in the fossorial and digger lizards due to their lower N, K, L, Q-Modules, Qmax, and higher D than those of other species. Thus, the skulls of fossorial (digger) species are morphologically more complex and have evolved higher functional efficiency and morphological complexity than those of other species. The Mann-Whitney U test supports that fossorial and digger lizards (n = 6) and other species (n = 52) differ from each other (Table 1).
In the morphological categories, the presence or absence of the upper temporal bars does not appear to be explained in parameters. In the PCA plots, the groups with the upper temporal bars cluster, while those without the upper temporal bars are scattered (Fig. 2e, k). The FDA on upper temporal bars shows a misclassification error rate of 34.48%, which indicates that the presence of upper temporal bars has no effect on the parameters (Fig. 3a). On the other hand, the presence or absence of postorbital bars does not have a strong association with the differences in the parameters, where the distribution of each group overlaps in the PCA plots (Fig. 2f, l). Nonetheless, the group without postorbital bars tends to score greater PC 1 values. On the other hand, the group with postorbital bars scores lower PC 1 and greater PC 2 values, while the group with an incomplete postorbital bar tends to scores lower PC 1 and PC 2 values. The FDA on postorbital bars indicate a misclassification error rate of 20.69%, suggesting that the presence of postorbital bars has a weak effect on the parameters (Fig. 3b).