High-Density Geometric Morphometric Analysis of Intraspecific Cranial Integration in the Barred Grass Snake (Natrix helvetica) and Green Anole (Anolis carolinensis)

Synopsis How do phenotypic associations intrinsic to an organism, such as developmental and mechanical processes, direct morphological evolution? Comparisons of intraspecific and clade-wide patterns of phenotypic covariation could inform how population-level trends ultimately dictate macroevolutionary changes. However, most studies have focused on analyzing integration and modularity either at macroevolutionary or intraspecific levels, without a shared analytical framework unifying these temporal scales. In this study, we investigate the intraspecific patterns of cranial integration in two squamate species: Natrix helvetica and Anolis carolinensis. We analyze their cranial integration patterns using the same high-density three-dimensional geometric morphometric approach used in a prior squamate-wide evolutionary study. Our results indicate that Natrix and Anolis exhibit shared intraspecific cranial integration patterns, with some differences, including a more integrated rostrum in the latter. Notably, these differences in intraspecific patterns correspond to their respective interspecific patterns in snakes and lizards, with few exceptions. These results suggest that interspecific patterns of cranial integration reflect intraspecific patterns. Hence, our study suggests that the phenotypic associations that direct morphological variation within species extend across micro- and macroevolutionary levels, bridging these two scales.


Introduction
Physical and molecular interactions among structures dictate how anatomical structures develop and evolve. The notion of modularity and integration has provided an analytical framework for investigating how these interactions contribute to phenotypic changes during ontogeny and evolution ( Olson & Miller 1958 ;Wagner 1996 ;Klingenberg 2008Klingenberg , 2014Goswami et al. 2014 ;Evans et al. 2022 ). For example, a more modular structure has been historically thought to promote the evolution of novel and disparate morphologies, where its components are able to evolve semi-autonomously from one another ( Wagner & Altenberg 1996 ). More recent simulation studies have demonstrated that inte-gration could indeed promote the evolution of extreme phenotypes if selection acts along the preferred axes of variation ( Goswami et al. 2014 ;. Given these conceptual scenarios, understanding how patterns of modularity and integration within species translate into clade-wide patterns would help bridge micro-and macroevolutionary trends. With broad interspecific sampling, recent macroevolutionary studies have examined the links between phenotypic integration and morphological diversification (e.g., Porto et al. 2009 ;Goswami et al. 2014 ; ; Bardua et al. 2020 ;Fabre et al. 2020 ;Felice et al. 2020 ). In concept, these large-scale evolutionary patterns of integration are expected to reflect the genetic control (e.g., pleiotropy), developmental mechanisms, and functional coordination that occur within individuals. Because integration patterns may be altered throughout development ( Hallgrímsson et al. 2009 ), the link between developmental patterns of integration and macroevolutionary patterns could be difficult to evaluate. In comparison, the pattern of modularity and integration between populations of a species sampled at equivalent growth stages (static modularity ;Mitteroecker 2009 ;Klingenberg 2014 ;Evans et al. 2022 ), particularly somatically mature individuals, is likely to be more stable in signal due to completion of major anatomical changes, and thus more readily comparable to macroevolutionary trends.
However, most studies on modularity and integration to date have analyzed either evolutionary or intraspecific patterns, without a unified framework to bridge these two levels. Using the same morphometric scheme across intraspecific and interspecific levels, some studies have found generally congruent patterns between select species and clade-wide patterns in caecilians ( Marshall et al. 2019 ) and birds ( Mitchell et al. 2021 ). Within squamates, Uroše vi ć and colleagues (2019) used a common landmarking scheme to investigate patterns of modularity across static, ontogenetic, and evolutionary levels within lacertids and found overlapping patterns. In this study, we explore the variational integration within the squamate species Natrix helvetica (barred grass snakes) and Anolis carolinensis (green anoles). These two taxa were selected because they are collectively abundant in natural history collections and the genera are commonly used as genetic and phylogeographic models ( Tollis et al. 2012 ;Fritz & Schmidtler 2020 ). Intraspecific patterns of cranial integration have been explored in these two taxa previously ( Sanger et al. 2012 ;Andjelkovi ć et al. 2017 ), but without sampling taxa beyond the genus level to analytically bridge micro-and macroevolutionary scales. Here, the landmark scheme employed on N. helvetica and A. carolinensis is equivalent to the high-density geometric morphometric data collected for a squamate-wide study  ). This experimental design allows for proper comparisons between intraspecific and interspecific patterns of integration, while sampling the genera Natrix and Anolis permit comparisons with results of previous investigations on these two taxa ( Sanger et al. 2012 ;Andjelkovi ć et al. 2017 ). Through these comparisons, we assess (i) differences in intraspecific cranial integration patterns in N. helvetica and A. carolinensis ; and (ii) correspondence of intraspecific integration patterns to respective interspecific patterns in snakes and non-snake extant squamates (hereby "lizards").

Specimens and imaging
We sampled 29 N. helvetica specimens and 41 A. carolinensis specimens for this study ( Table 1 ), where one specimen from each taxon was used as template. Specimens exhibiting somatic maturity were selected, with visually intact head. Because the skull of snakes is highly kinetic, specimens that showed symmetric and closed or nearly closed gape were selected to reduce variation in cranial shape due to jaw position. The Anolis specimens were scanned using a GE Phoenix v|tome|xm CT Scanner at the University of Florida. N. helvetica specimens were imaged with Nikon Metrology HMX ST 225 μCT scanner at the Natural History Museum, London (NHMUK). The segmentation of the cranial elements for both species were performed using Avizo (Thermo Fisher Scientific) to digitally segment and reconstruct the skulls. Upon imaging, many N. helvetica specimens were found to unilaterally lack the maxillary bone, presumably removed for another study. However, the remaining side was intact. The resulting PLY mesh files were then imported into GeoMagic Wrap (formerly 3D Systems, now Artec 3D) to remove extraneous objects (e.g., vertebrae, mandible), and we virtually filled in foramina to prevent error in projection of surface semi-landmarks. If the left maxilla was missing from the specimen, then the skull models were mirrored (mirrored specimens specified in Table 1 ). For N. helvetica , we virtually removed teeth from the left side of the pterygoid and palatine bones to eliminate obstructions for the placement of surface semi-landmarks on these elements. Once the teeth were removed, the resulting cavities on the bones were fil le d using the Fill Single tool and its Tangent option in GeoMagic Wrap. The QuickSmooth function was used to apply global smoothing to the mesh surfaces. The skull models were then exported as PLY mesh files.

Coordinate data
The IDAV Landmark Editor ( Wiley et al. 2005 ) was used to place fixed (discrete) and curved (semi-)landmarks to delineate the following osteological elements : premaxilla, nasal, maxilla, jugal in A . carolinensis , frontal, parietal, squamosal in A. carolinensis , supratemporal in N. helvetica , jaw joint of the quadrate, supra-otoccipital, basioccipital, pterygoid, palatine, and occipital condyle. The jugal is excluded in N. helvetica and the snake-wide dataset due to the absence of this bone in snakes. After placing fixed landmarks on both the left and right sides of the skull and curved semi-landmarks on only the right side, the coordinate data were exported as a .pts file. Using R language Table 1. List of specimens sampled for this study. Mirrored specimens refer to those that were digitally mirrored prior to landmarking due to damaged or missing elements on the right side of the skull. Please refer to Watanabe et al. (2019) for information of specimens sampled for the interspecific datasets.   Bardua et al. 2019 ). The fixed and curve (semi-)landmark scheme is identical to the one used previously for a squamate-wide analysis of cranial integration . The surface (patch) semi-landmarks were placed onto the mesh surface using the placePatch function in the Morpho R package ( Schlager 2017 ) using N. helvetica (NHMUK 1917.6.22.1) and A. carolinensis (FMNH 242298) specimens as templates. With the exception of surface semi-landmarks on the pterygoid and palatine bones for N. helvetica , these templates contain the same number of surface semi-landmarks as the hemispherical template used previously for the squamate-wide data  ). These template specimens were not included in the subsequent analysis to maintain equivalency in landmarking procedure across analyzed specimens. In the previous study, surface semi-landmarks were not placed on the pterygoid and palatine bones in snakes due to the presence of teeth, which would lead to inconsistent placement of points. As mentioned above, we virtually removed the palatal teeth and closed the alveoli in GeoMagic Wrap to enable the placement of surface semi-landmarks on the pterygoid and palatine bones in N. helvetica .

Shape data
Prior to aligning the coordinate data, we used the mirrorfil l function in the paleomorph R package to mirror the curve and surface semi-landmarks placed only on the right side of the skull based on position of bilateral and midline fixed landmarks. This allows creation of artificially bilateral coordinate data to minimize artifacts of aligning one-sided data for bilaterally symmetric structures ( Cardini 2016( Cardini , 2017 ; Bardua et al. 2019 ).
Following the mirroring process, we performed generalized Procrustes analysis with sliding semi-landmarks, while minimizing total bending energy ( Gunz et al. 2005 ;Gunz & Mitteroecker 2013 ) on the bilateral coordinate data. Once the generalized Procrustes analysis was performed, semi-landmarks and fixed landmarks on the left-side were removed. Following the curve subsampling and patching process, intraspecific coordinate data of 28 N. helvetica specimens (i.e., with the template specimen removed) included 41 fixed, 535 curve, and 549 surface (semi-)landmarks, while those of 40 A. carolinensis comprised 47 fixed, 595 curve, and 580 surface (semi-)landmarks ( Fig. 1 ; Table 2 ). Prior to main analysis, we used the procD.lm function in the geomorph R package ( Adams & Otárola-Castillo 2013 ) to check for potential biases in our dataset due to surgical removal of one of the maxillae in some of the Natrix specimens.
The results indicated a lack of significant skull shape differences between Natrix with and without intact maxillae ( R 2 = 0.0196; P = 0.808).

Analyses
Analyses were performed using R v4.1.0 ( R Core Development Team 2022 ). To visualize the pattern of cranial shape variation, we created morphospaces based on the first two principal components (PCs). The strength of integration within and between cranial elements was calculated through the covariation-ratio (CR) method ( Adams 2016 ) and EMMLi, a maximum-likelihood approach ( Goswami & Finarelli 2016 ). Greater ρ (rho) value from EMMLi analysis and CR value indicate greater correlation and covariation. We assessed the relationship between cranial shape to log-transformed centroid size using procD.lm function in geomorph R package on shape against log-transformed centroid size. For a quantitative comparison of intra-and interspecific integration patterns, a correlation analysis was conducted on CR and EMMLi values between N. helvetica and snake-wide data, and between A. carolinensis and lizard-wide data.

Results
Morphospace PC analysis reveals that the first two PC axes together account for 54.2 and 38.2% of total skull shape variation in N. helvetica and A. carolinensis , respectively ( Fig. 2 ). In N. helvetica , PC1 is associated with the anteroposterior extent of the neurocranium, with high variation in the relative position of the jaw joint. PC2 corresponds to the mediolateral position of the maxilla, pterygoid, and the jaw joint. In A. carolinensis , shape changes along PC1 comprises relative anteroposterior length of the rostrum and occiput region, mediolateral extent of basioccipital, and the anteroposterior position of the quadrate relative to the rest of the skull. PC2 broadly correlates with the relative length and depth of the skull.

Allometry
Allometry is a significant factor in both species ( P < 0.05). In N. helvetica and A. carolinensis , allometry accounts for 8.2 and 24.9% of their respective total skull shape variation.

Intraspecific cranial integration
Based on CR and EMMLi analyses, N. helvetica and A. carolinensis exhibit shared as well as distinct patterns in cranial integration ( Table 3 ; Fig. 3 ). In N. helvetica , we observe strong integration (in upper 20% of range in CR and ρ values) between the occipital elements (supra-otoccipital, basioccipital, and occipital condyle) and jaw joint on the quadrate bone. The frontal and parietal bones exhibit elevated CR and ρ values (CR = 0.80; ρ = 0.41). The jaw joint on the quadrate has the greatest within-region correlation ( ρ = 0.98) of the cranial bones. The lowest within-region correlation value was demonstrated by the basioccipital ( ρ = 0.60). With respect to CR and ρ values, the covariation between the adjacent regions supra-otoccipital and occipital condyle is the strongest (CR = 0.96; ρ = 0.77). Both the frontal-parietal and pterygoid-palatine pairs exhibit moderate level of integration (near the middle range of CR and ρ values). The relative strength of cranial integration among partitions in the posterior skull is generally reduced in the EMMLi analysis, however the occipital condyle remains a major conduit in the posterior skull. In the viscerocranium, relatively strong integration is found between premaxillafrontal (CR = 0.66; ρ = 0.39) and premaxilla-nasal (CR = 0.72; ρ = 0.37). Similar to N. helvetica , A. carolinensis also shows strong between-region integration among occipital and palate elements ( Table 3 ; Fig. 3 ). For example, the occipital condyle and the basioccipital show the highest Table 3. Between-partition covariance ratios (upper triangle), between-partition correlations (bottom triangle), and within-partition correlations (diagonal) across cranial regions. Values calculated using modularity.test in the geomorph R package and EMMLi, respecti vel y. Greater ρ and CR values denote greater level of integration between pairs of elements. Darker shades of grey on off-diagonal elements denote relati vel y stronger integration, where five degrees of shades were applied based on range of ρ and CR values divided into five equal bins.

Natrix helvetica
Cranial integration in Natrix and Anolis 9 CR value for between-region integration, with a relatively high ρ value (CR = 0.90; ρ = 0.40). As in N. helvetica , the EMMLi analysis shows the greatest withinregion integration for the jaw joint on the quadrate ( ρ = 0.96). Likewise, A. carolinensis shares moderate to relatively strong integration between the frontal and parietal bones (CR = 0.784 in ≥60 percentile of the range in value; ρ = 0.31 in > 40 percentile). In contrast to N. helvetica , A. carolinensis shows moderate to relatively strong integration between the rostral bones, including the premaxilla and maxilla (CR = 0.80 in ≥80 percentile; ρ = 0.36 in > 60 percentile), and between the pterygoid and palatine bones based on the value ( ρ = 0.47; > 80 percentile). Taken together, A. carolinensis displays elevated levels of integration between the premaxilla and maxilla; frontal and parietal; pterygoid and palatine; and occipital regions (supra-otoccipital, basioccipital, and occipital condyle). In addition, jaw joint and occipital condyle as well as pterygoid and jaw joint exhibit heightened degree of integration as in N. helvetica .

Comparison to interspecific integration pattern
The comparison of the intraspecific pattern of integration in N. helvetica to the results of a phylogenetically corrected snake-wide analysis  shows key similarities. Notably, the presence of moderate to strongly integrated "occiput" and "palate" regions, with modular viscerocranium is com-mon between the N. helvetica and snake-wide datasets.
In addition, the frontal and parietal bones show relatively strong between-region integration in N. helvetica (CR = 0.80; ρ = 0.41) and the interspecific data (CR = 0.60; ρ = 0.39). Furthermore, the interspecific (snake-wide) data showed the jaw joint on the quadrate to have strong within-region integration ( ρ = 0.99) and the basioccipital to have the lowest within-region integration ( ρ = 0.61). This outcome is replicated in N. helvetica, where the jaw joint on the quadrate has the highest within-region integration ( ρ = 0.98) and the basioccipital has the lowest within-region integration ( ρ = 0.60). When compared to the interspecific lizard data, the pattern of cranial integration within A. carolinensis is largely congruent, with some exceptions. Notably, the intraspecific pattern shares moderately to strongly integrated "rostrum," "skull roof," "palate," and "occiput" regions that were reported for lizardwide analysis of cranial integration . Between-region integration from the lizard interspecific dataset and A. carolinensis data indicate weak integration between the premaxilla-squamosal, nasaloccipital condyle, frontal-squamosal, frontal-palatine, and squamosal-palatine bones. The maxilla ( ρ = 0.54) and supra-otoccipital ( ρ = 0.54) present the lowest within-region integration for the lizard-wide dataset. Anolis carolinensis exhibits similar within-region integration for the maxilla ( ρ = 0.45) and supra-otoccipital In the top row, the thickness of lines and diameter of circles indicate the strength of correlation ( ρ) between and within partitions, respecti vel y, from EMMLi analysis. In the bottom row, the thickness of lines is proportionate to CR values. Abbreviations: pmx, premaxilla; ns, nasal; jg, jugal; mx, maxilla; fr, frontal; par, parietal; sq, squamosal; st, supratemporal, qd, quadrate (jaw joint); so, supra-otoccipital; bo, basioccipital; pt, pterygoid; pl, palatine; oc, occipital condyle.
( ρ = 0.52), with the nasal bone also showing weak within-region integration ( ρ = 0.50). The jaw joint on the quadrate shows the strongest within-region integration interspecifically ( ρ = 0.98) and intraspecifically ( ρ = 0.96). Meanwhile, we observe relatively weak integration between the occipital condyle and the jaw joint (quadrate) across "lizards" that is not reflected in the intraspecific A. carolinensis analysis. The A. carolinensis specimens exhibit a strong integrating signal between the occipital condyle and the quadrate (CR = 0.83; ρ = 0.50; > 80 percentile), whereas the interspecific lizard data have a relatively low level of integration (CR = 0.689; ρ = 0.36).
Quantitative comparison of intra-and interspecific patterns in cranial integration supports the results described above and provides new insights. First, correlation analysis on CR and ρ values show that there is a broad correspondence in between-region integration when comparing intra-and interspecific datasets. While deviations from the best fit regression line are apparent ( Fig. 4 ), the relative degree of integration is strongly correlated between N. helvetica and snakewide datasets ( R 2 ρ = 0.723, P < 0.0001; R 2 CR = 0.512, P < 0.0001). Likewise, the correlation between A. carolinensis and lizard-wide dataset with respect to CR and ρ values is statistically significant ( R 2 ρ = 0.243, P < 0.0001; R 2 CR = 0.122, P = 0.0017). As reflected in the results of the correlation analysis, the bivariate plots of intra-and interspecific integration ( Fig. 4 ) depict tighter correspondence in CR and ρ values for N. helvetica and snake-wide datasets, compared to those of A. carolinensis and lizard-wide data. Moreover, the plots reveal pairs of skull regions that deviate from the association between intra-and interspecific patterns. For example, N. helvetica exhibits greater ρ values for premaxilla-frontal integration than the snake-wide pattern, whereas snakes show relatively stronger nasalmaxilla, frontal-quadrate, and basioccipital-pterygoid integration than N. helvetica ( Fig. 4 a, b). Among lizards, A. carolinensis presents stronger degree of integration between the occipital condyle with the premaxilla, jaw joint, and pterygoid than in lizard-wide data, whereas frontal-occipital, pterygoid-palatine, and supra-otoccipital-occipital condyle integration is relatively stronger across lizards than within A. carolinensis ( Fig. 4 c, d).

Intraspecific cranial integration
Compared to a previous study on cranial integration in N. natrix and N. tessellata ( Andjelkovi ć et al. 2017 ), our results for N. helvetica show both similarities and differences. For the regions characterized here, the previously published dataset also showed statistically significant covariations between the braincase (partitioned as single anatomical unit) and quadrate, as well as a strong association between the maxilla and palatine bones. The maxilla, pterygoid, and the jaw joint are highly mobile elements in snake skulls. As such, their functional Bivariate plots of between-region integration of interspecific against intraspecific data. (a) ρ values for snake-wide and Natrix data; (b) CR values for snake-wide and Natrix data; (c) ρ values for lizard-wide and Anolis data; and (d) CR values for lizard-wide and Anolis data. Black lines and grey bands indicate best fit regression lines and 95% confidence intervals, respecti vel y. Text labels denote the pair of regions associated with the integration value (some labels in dense clusters of points were removed to improve legibility). Abbreviations: pmx, premaxilla; ns, nasal; jg, jugal; mx, maxilla; fr, frontal; par, parietal; sq, squamosal; st, supratemporal, qd, quadrate (jaw joint); so, supra-otoccipital; bo, basioccipital; pt, pterygoid; pl, palatine; oc, occipital condyle.
properties may underlie the strong integration patterns across these regions. For instance, PC2 axis is associated with the relative position of these mobile elements, suggesting that these regions undergo coordinated changes in shape and relative position. The more moderate integration between the palatine and pterygoid in Natrix and snake-wide data, as replicated in our dataset, could be due to the low contact point between these two bones in snakes. In contrast to Andjelkovi ć et al.'s (2017) , results from other Natrix species, our N. helvetica data lack evidence of strong integration between the maxilla, quadrate, and supratemporal elements. This result could represent genuine differences in cranial integration between Natrix species. However, these discrepancies could also be due to the differing landmarking ap-proach and metrics for evaluating phenotypic integration. Andjelkovi ć and colleagues employed only fixed landmarks when characterizing skull shape and identified integrated pairs of regions based on statistical significance using the RV coefficient ( Escoufier 1973 ;Klingenberg 2013 ). Results of integration analysis are known to be influenced by the type of landmarks (e.g., inclusion of semi-landmarks), number of landmarks, and specimen sampling (e.g., Adams 2016 ;Cardini, 2019 ). However, the general principle that separated elements of the rostrum and palate in snakes are more weakly integrated than more connected elements, such as the neurocranium, remains true across both studies.
The integration pattern found here for A. carolinensis also shows a mixture of congruent and contrasting results with a previous comprehensive study on cranial integration in Anolis species ( Sanger et al. 2012 ). Using two-dimensional geometric morphometric data, Sanger and colleagues (2012) found that the "tripartite" model of integration is supported for A. carolinensis , where the rostrum, orbit, and braincase form integrated regions. Although the difference in landmarking scheme and sampling preclude direct comparison of results, the integration patterns observed are generally consistent. Although our results lack evidence for an "orbit" module, the braincase, including the occiput, and the rostrum are strongly integrated. It is worth noting that our A. carolinensis sampling consists predominantly of individuals from Alachua County in Florida, with the exception of two specimens (FLMNH 2990-8,101848). While there is a possibility that these two individuals may be members of populations with different cranial integration patterns, the morphospace ( Fig. 2 b) shows that their skull shapes are well within the variation exhibited by other specimens. Therefore, we consider the cranial integration pattern presented here for A. carolinensis to be a reliable, at least for the pattern exhibited by populations in Alachua County.

Comparison with interspecific integration
The intraspecific patterns of cranial integration within A. carolinensis and N. helvetica resemble the interspecific patterns observed across lizards and snakes, respectively. These two distinct levels of analyses present a shared set of highly integrated regions, including the "skull roof" (frontal, nasal), "palate" (pterygoid, palatine), and "occiput" (supra-otoccipital, basioccipital, and occipital condyle). Moreover, A. carolinensis and lizard-wide analyses show a more integrated rostrum than N. helvetica and other snakes. In contrast, N. helvetica specimens have a more modular rostral region, in agreement with the interspecific snake-wide pattern . This difference in the rostral region between snakes and lizards could be related to the level of articulation among rostral elements. In snakes, rostral elements, including the premaxilla, nasal, and maxilla bones, do not closely articulate with one another as in lizards.
The overall congruence in integration patterns of N. helvetica and A. carolinensis to the snake-and lizardwide analyses is notable. This outcome suggests that the processes that shape trait covariation at the population level extends to large-scale evolutionary processes. Previous studies that have examined intraspecific and interspecific patterns of integration and modularity generally corroborate this finding. For instance, the skull of ring-necked parakeet supports the same nine-module hypothesis observed across crown birds ( Mitchell et al. 2021 ), and interspecific and static patterns of cranial in-tegration are similar in caecilians ( Marshall et al. 2019 ). In lacertids, datasets capturing static and interspecific variation supported different integration hypotheses ("anteroposterior" and "neurodermatocranial", respectively; Uroševi ć et al. 2019 ). However, these modularity hypotheses are quite similar in the partitioning of landmarks into units, implying that intra-and interspecific patterns of cranial integration match to a large extent.
Consequently, phenotypic covariation within species could be used to predict covariation across clades, and vice versa, thereby unifying micro-and macroevolutionary dynamics. This conclusion that integration pattern and modules are conserved across multiple hierarchical levels seemingly contradicts the observation that pattern of modularity changes with cranial shape across Anolis species, potentially due to differences in diet ( Butler & Losos 2002 ;Sanger et al. 2012 ). However, these contrasting implications are not necessarily incompatible. We acknowledge and fully expect that the pattern and strength of phenotypic covariation evolves along and differ across lineages. In fact, there are differences in intraspecific and interspecific patterns of cranial integration, as noted above, that could signify how the sampled species have deviated from the overall clade-wide patterns. Importantly, this study demonstrates that the strongly integrated elements are broadly maintained across micro-and macroevolutionary timelines, with limited number of exceptions (e.g., highly integrated quadrate-occipital in A. carolinensis that is absent in lizard-wide pattern). Individual lineages are still able to modify the degree of integration from this conserved pattern.

Allometry and sexual dimorphism
Allometry is often a strong integrating force in biological structures ( Zelditch & Fink 1995 ;Klingenberg 2013 ;Bright et al. 2016 ). Both N. helvetica and A. carolinensis show statistically significant allometric signal in their skull shapes, with greater effect of size in the latter. Although some studies on modularity and integration have performed analyses on allometry-corrected data ( Sanger et al. 2012 ; Uroše vi ć et al. 2019 ), we found that analysis of allometry-corrected data strengthens the relative covariation and correlation of other pairs of elements (Fig. S1). While these results may include biologically informative information, the greater number of similarly integrated pairs of elements suggests that removing the overall allometric signal, known as a strong integrating factor, may dilute the original integration signal ( Zeldith & Fink 1995 ;Klingenberg 2009 ; but see Uroševi ć et al. 2020 ). In addition, keeping the allometric signal maintains the equivalency of our intraspecific analysis with the squamate-wide analysis. In the top row, the thickness of lines and diameter of circles indicate the strength of correlation ( ρ) between-regions and within-regions, respecti vel y, from EMMLi analysis. In the bottom row, the thickness of lines is proportionate to CR values. Abbreviations: pmx, premaxilla; ns, nasal; jg, jugal; mx, maxilla; fr, frontal; par, parietal; sq, squamosal; st, supratemporal, qd, quadrate (jaw joint); so, supra-otoccipital; bo, basioccipital; pt, pterygoid; pl, palatine; oc, occipital condyle. Therefore, we did not focus on reporting the results from the allometry-corrected analysis.
Sexual dimorphism is known to occur in species of Natrix ( Andjelkovi ć et al. 2016 ) and Anolis ( Schoener 1967 ;Losos 2011 ). Ecological diversity is produced by sexual differences in Anolis species as males are typically seen on larger and higher perches compared to females and due to differences in the ecological niche of males and females ( Schoener 1967 ;Butler et al. 2007 ). This difference in size across species and between sexes may contribute to changes in modularity and integration regarding functional aspects such as locomotion and feeding. Unfortunately, we lacked sex information on sampled A. carolinensis specimens, preventing a proper assessment of sex-related differences in shape and integration. While we do not see a visual indication of cranial shape dimorphism in the morphospace for A. carolinensis ( Fig. 2 ), allometry accounts for nearly a quarter of the skull shape variation. As such, known sexual size dimorphism in Anolis is expected to be associated with skull shape differences, which has been demonstrated previously ( Herrel et al. 2007 ). Whether this shape difference would lead to divergence in integration patterns is unclear. For N. helvetica , which sex information was known ( Table 1 ), we observed weak support for sexual shape dimorphism in the skull ( R 2 = 0.082; P = 0.041). This result counters a previous study showing lack of sexual shape dimorphism in the dorsal aspect of the head in N. helvetica ( Tamagnini et al. 2018 ), but is concordant with known sexual shape dimorphism reported for N. natrix and N. tessellata , where different sets of bones were sexually dimorphic between these two species ( Andjelkovi ć et al. 2016 ). Despite the presence of sexual shape dimorphism, separate CR and EMMLi analyses on female and male datasets demonstrate equivalent patterns, with the exception of stronger relative integration between the frontal and parietal bones in males ( Fig. 5 ; Table S1). While there may be subtle differences in cranial integration between sexes (e.g., relatively more integrated palatine-pterygoid and frontal-parietal in males), the results and conclusions of this study remain broadly consistent.

Conclusions
In this study, we use a high-density geometric morphometric approach to document the pattern of cranial integration in N. helvetica and A. carolinensis . We then compare these intraspecific patterns to snake-and lizard-wide cranial integration, respectively, from a previous study using a comparable landmarking scheme. We find that intraspecific patterns of cranial integration in N. helvetica and A. carolinensis species share the following integrated regions: skull roof (frontal, parietal), occiput (supra-otoccipitalbasioccipital-occipital condyle-quadrate), and palate (palatine-pterygoid). Anolis carolinensis shows a more integrated rostrum (premaxilla-nasal-maxilla) than in N. helvetica specimens. Overall, the intraspecific patterns of covariation in N. helvetica and A. carolinensis generally match the snake-and lizard-wide interspecific patterns, respectively, with few exceptions that likely represent lineage-specific deviations from clade-wide trends. Therefore, cranial integration within species may direct clade-wide morphological evolution to a large and detectable degree; thus, unifying micro-and macroevolutionary trends, and potentially process, in phenotypic evolution.