High-density three-dimensional morphometric analyses support conserved static (intraspecific) modularity in caecilian (Amphibia: Gymnophiona) crania

Analyses of phenotypic integration and modularity seek to quantify levels of covariation among traits to identify their shared functional, developmental and genetic underpinnings (‘integration’), which may delineate semi-independent subsets of highly integrated traits (‘modules’). Existing studies have focused mainly on mammals or model organisms, limiting our understanding of factors that shape patterns of integration and modularity and their importance for morphological evolution. We present the first study of static (intraspecific) integration and modularity in caecilian crania, using dense surface sliding semi-landmarks to quantify cranial morphology in Boulengerula boulengeri and Idiocranium russeli. Eleven hypotheses of modular organization were compared with a likelihood approach and best-fitting models confirmed with covariance ratio analysis. Allometric corrections and subsampling analyses demonstrated the robustness of results. Allometry had a substantially larger influence on the results of landmarkonly analyses relative to analyses incorporating semi-landmarks. Idiocranium russeli displayed significantly higher variation than B. boulengeri, but they had similar 12and 13-module patterns, respectively, suggesting that caecilian crania are highly modular and that the modularity pattern is largely conserved across these species, despite their divergent morphologies and > 175 Myr of independent evolution. As in previous mammalian studies, our caecilian study suggests that cranial modularity patterns are conserved within major vertebrate clades.


INTRODUCTION
The concepts of integration and modularity have been explored within a diverse range of disciplines for many years and have become a major focus of study in evolutionary biology in recent decades (Goswami, 2006a;Klingenberg, 2008Klingenberg, , 2014. The two concepts are inherently linked, with 'morphological integration' referring to the cohesion (covariation) among traits attributable to shared ontogeny or function, whereas 'modularity' refers to the idea that these traits can be grouped into semi-autonomous sets of highly correlated characters ('modules'), which are strongly related internally while being relatively independent from each other (Olson & Miller, 1958;Wagner, 1996;Winther, 2001;Hallgrímsson et al., 2002;Klingenberg, 2004Klingenberg, , 2008Klingenberg, , 2014Goswami, 2006a;Parr et al., 2016). Morphological integration is often categorized according to the factors expected to be influencing trait correlations, for example: 'genetic integration', the coinheritance of traits; 'developmental integration', the shared developmental influences on traits; and 'functional integration', the physical interactions between traits required by a shared functional task (Goswami et al., 2014). However, both integration and modularity within a structure are likely to result from a combination of genetic, developmental and functional interactions (Cheverud, 1982;Wagner, 1996;Wagner & Altenberg, 1996;Zelditch et al., 2004;Goswami, 2006b;Porto et al., 2009;. Modular organization has been identified at multiple levels and in diverse systems, from gene regulatory networks to colonies (Klingenberg, 2004;Clune et al., 2013;Laurin, 2014), with several studies focusing on identification and quantification of the relationships between genetic, developmental and functional modules (Goswami et al., 2009). The identification and definition of 'modules' in statistical terms varies, and here we use this term to indicate highly integrated regions in which the correlations among traits within regions are greater than those between regions. As detailed further below, the hierarchy of correlations among regions also provides rich information on cranial integration beyond any single definition or statistic describing those relationships.
Caecilians (order Gymnophiona) are a monophyletic group of limbless, elongate, superficially snake-like vertebrates constituting one of the three orders of modern amphibians. The 212 currently recognized species of extant caecilians are classified into ten families and found in mainly tropical regions of Central and South America, Sub-Saharan Africa, the Seychelles and South and Southeast Asia Kamei et al., 2012). With the exception of the few aquatic and semi-aquatic species of the Typhlonectidae (Taylor, 1968;Nussbaum & Wilkinson, 1995), caecilians are terrestrial as adults and predominantly fossorial (Taylor, 1968;Himstedt, 1996;Gower & Wilkinson, 2005;Wilkinson, 2012). Caecilians have a generally cryptic nature, which, along with their (at least historically) relative inaccessibility in the tropics and subtropics, has contributed to their diversity and ecology being notoriously understudied and often poorly understood (Gower & Wilkinson, 2005;Wilkinson, 2012;Sherratt et al., 2014). Although it is widely accepted that all adult caecilians practise some head-first burrowing, morphology and observations in both the field and the laboratory suggest that species vary from being dedicated burrowers to having relatively more surface-active lifestyles (Nussbaum, 1977;Ducey et al., 1993;Gower et al., 2004;Wilkinson, 2012). Despite the constraints on morphological diversification of the skull that might be expected to be imposed by the functional demands of headfirst burrowing, caecilian crania have been shown to be diverse and variable in shape and composition of elements (Taylor, 1969;Nussbaum, 1977;Wake, 2003;Sherratt et al., 2014). In addition to burrowing, caecilian crania are important in feeding, gas exchange, locomotion, sensory perception and, potentially, social interactions; their structure is thus expected to be shaped by many competing functional demands.
Previous work has suggested that patterns of cranial modularity are highly conserved in other vertebrate clades, such as therian mammals (Goswami, 2006a). A six-module pattern has been reported for placental mammals (Cheverud, 1982;Goswami, 2006a;Porto et al., 2009;Goswami & Finarelli, 2016), and Felice & Goswami (2018) found a seven-module structure in the avian skull that showed some concordance with mammalian structures. Studies of modularity in amphibians have generally found support for more highly integrated crania, for example, two or four developmental modules in salamanders (Laurin, 2014). A study on frogs (Vidal-García & Keogh, 2017) found no strong support for a modular structure, and Simon & Marroig (2017) found a range of best-supported modularity patterns across the toad species of the Rhinella granulosa complex. Previous interspecific studies on modularity in caecilians found support for independence of the snout from the rest of the cranium (Sherratt, 2011) and for a highly partitioned, ten-module structure consisting of frontal, parietal, dorsal nasopremaxilla, palatal nasopremaxilla, maxillopalatine (combining lateral, interdental and palatine shelf of the maxillopalatine), occipital (combining occipital and occipital condyle regions), quadrate-squamosal (combining jaw joint, squamosal and quadrate regions), ventral os basale-vomer, pterygoid and stapes modules (Bardua et al., 2019). There are no existing studies examining intraspecific (static) modularity in caecilians. Hence, we aim to fill this gap, and specifically, we aim to investigate whether or not the same patterns are found intraspecifically as have been identified by Bardua et al.'s (2019) cladewide work in caecilians, and whether the same pattern of modularity is conserved across two different species, Boulengerula boulengeri and Idiocranium russeli. We also aim to investigate whether the degree and pattern of cranial integration and modularity restricts or facilitates variation in traits. Under the constraint hypothesis, wherein high integration is hypothesized to restrict shape change in specific directions, modules displaying relatively high integration (high withinmodule trait correlations) would exhibit relatively lower levels of variation than observed in more weakly integrated regions (Goswami & Polly, 2010;Goswami et al., 2014;Felice & Goswami, 2018). In contrast, under the facilitation hypothesis, high integration of traits facilitates shape change, such that modules with high within-module correlations would exhibit high variation (Parr et al., 2016;Randau & Goswami, 2017;Larouche et al., 2018).
In this study, we quantify the magnitude of integration in two African caecilian species from the families Herpelidae and Indotyphlidae . The herpelid B. boulengeri Tornier (1896), is endemic to the Usambara Mountains, Tanzania (Gower et al., 2004;Loader et al., 2011), whereas the indotyphlid I. russeli Parker is endemic to Cameroon (Gower et al., 2015). These species are believed to have shared a common ancestor > 175 Mya (San Mauro et al., 2014: supplementary data). Idiocranium russeli is one of the smallest caecilians, with a maximal total length of ~170 mm (Gower et al., 2015), compared with ~280 mm in B. boulengeri (Gower et al., 2004). Idiocranium russeli is considered to be miniaturized (Wake, 1986;Theska et al., 2018), which may be associated with substantial reductions in the skull, and possibly increased variation, especially in anterior elements, such as the frontal and nasopremaxillary bones (Hanken, 1984). Although both species in this study are oviparous direct developers (San Mauro et al., 2014), B. boulengeri hatchlings have vernal teeth (sensu San Mauro et al., 2014), indicating the presence of maternal dermatophagy (skin-feeding behaviour) (Kupfer et al., 2006;San Mauro et al., 2014). The incorporation of a dermatophagic stage might impact upon the ontogenetic trajectory of the skull, with potential consequences for adult cranial integration and variation of dentigerous regions involved in dermatophagy. For example, Müller (2007b: chapter 4) found that the dermatophagic Boulengerula taitanus appears to exhibit accelerated development of the premaxilla relative to the non-dermatophagic Gegeneophis ramaswamii and Hypogeophis rostratus, and that this tooth-bearing bone remains relatively larger in B. taitanus compared with G. ramaswamii throughout ontogeny. Boulengerula boulengeri can be found by digging in soil but are rarely found in pitfall traps or observed on the surface (Gower et al., 2004), which, together with their gross morphological features (including 'bullet-shaped' heads, reduced eyes covered with bone and strongly stegokrotaphic skulls) and studies of diet, has led to this species being interpreted as leading a predominantly dedicated burrowing, endogeic lifestyle (Jones et al., 2006). In contrast, in I. russeli the eye is small and covered with skin but not bone (Parker, 1936), and the skull can be interpreted as being relatively loosely constructed (Wake, 1986) and less 'bullet-shaped' (Parker, 1936) compared with B. boulengeri (Fig. 1). These features, together with some evidence that I. russeli is most readily found in very wet soils (Gower et al., 2015), suggest that I. russeli is a less dedicated burrower than B. boulengeri and is restricted to softer soils. Here, we investigate whether cranial integration and modularity are conserved across B. boulengeri and I. russeli or whether differences in lifestyle and morphology have resulted in divergences in patterns of cranial integration and variation between these species.

SpecimenS
This study is based on 41 specimens (22 males, 19 females) of B. boulengeri and 36 specimens (13 males, 18 females, five unknown) of I. russeli in the collections of the Natural History Museum, London (NHM) and University of Cambridge Museum of Zoology (full specimen details are presented in Supporting Information, Table S1). The heads of whole, alcoholpreserved specimens were imaged using a Nikon Metrology X-Tek HMX ST 225 micro-CT scanner at the NHM. Details of scan parameters are presented in the Supporting Information (Table S2). The three-dimensional (3D) volumes (tomographs) of 22 B. boulengeri specimens and four I. russeli specimens had previously been digitally dissected to create 3D isosurface models of cranial bone using VG Studio MAX v.2.0 (Volume Graphics, 2001), as described by Sherratt (2011); the remainder of the tomographs were post-processed using Avizo v.9.3 (FEI, Hillsboro, OR, USA). These latter 51 surface reconstructions were simplified using 'Quadric Edge Collapse Decimation' in MeshLab (Cignoni et al., 2008) to reduce data volume without compromising quality. All reconstructed surface meshes were digitally dissected in Geomagic Wrap 2014/2015 (http://www.geomagic.com; Accessed December 2018) to remove any remaining vertebrae and mandible elements, leaving only the bony elements of the upper skull. Given that this study was focused only on the cranial elements, removal of these Downloaded from https://academic.oup.com/biolinnean/article-abstract/126/4/721/5320147 by university of adelaide user on 13 February 2020 other parts to isolate the area of interest helps with accurate placement of landmarks and semi-landmarks on regions in articulation or close contact with noncranial elements. Three B. boulengeri and 12 I. russeli scans were mirrored using the 'mirror' function in Geomagic owing to damage on the right-hand side of the cranium, with all subsequent processing and analyses performed on the right side of the cranium only. All scans were then treated in Geomagic to fill surface holes and remove highly creased edges and spikes, which may interfere with the semi-automatic morphometric data capture approach used in this study.

Shape analySeS
The crania of both I. russeli and B. boulengeri are each composed of ten different types of cranial elements. We identified 17 individual cranial regions (Table 1), either single cranial elements or divisions of larger bones (e.g. lateral and ventral surfaces of the maxillopalatine), which could be identified and segregated consistently in every specimen. Bones were divided so as to represent potentially divergent functional regions of the individual elements where these could be differentiated using anatomical structures, such as tooth rows and muscle attachment ridges. These 17 regions were almost identical to those defined by Bardua et al. (2019), which will ultimately allow comparison of patterns across caecilians at both the static and evolutionary scales. The only difference between the regions in this study and those of Bardua et al., (2019) is the addition of a separate mesethmoid in our study, because this region was variably present across specimens in the clade-wide study, necessitating it to be grouped with the frontal region a priori. A 3D surface-based landmarking procedure similar to that of Dumont et al. (2016), Felice & Goswami (2018) and Bardua et al., (2019) was used to quantify the external surface morphology of the cranium. The surface meshes were first imported into IDAV Landmark Editor v3.6 (Wiley et al., 2005). Sixty-eight Type I landmarks (located at the discrete juxtaposition of two or three bones) and Type II landmarks (located at the maxima of curvatures) (Bookstein, 1991) (Fig. 1A) and 340 sliding semi-landmarks over 72 3D curves (Gunz et al., 2005) were digitized on the right side of all B. boulengeri specimens; 66 landmarks ( Fig. 1B) and 336 sliding semi-landmarks over 70 3D curves were digitized on the right side of all I. russeli specimens (Supporting Information, Tables S3 and S4). Landmarks and semilandmarks were the same across species except where different approaches were necessary to represent their morphologies adequately. Accordingly, two additional landmarks and two additional curves were digitized on B. boulengeri specimens to represent the shape of the maxillopalatine. Non-corresponding points were removed from the specimens before interspecific comparisons, resulting in the removal of two landmarks and 50 semi-landmarks on curves from B. boulengeri specimens, and 35 semi-landmarks on curves from I. russeli. The R package 'Morpho' (Schlager, 2017) in R v.3.3.2 (R Core Team, 2017) was used to re-sample semilandmarks along curves, resulting in 822 and 807 semi-landmarks equidistantly placed on curves for B. boulengeri and I. russeli, respectively (Supporting Information, Table S4). Identical landmarks and curves, along with 736 surface semi-landmarks (Supporting Information, Table S5), were placed onto one mesh for each species (Fig. 2), serving as 'templates', from which these surface semi-landmarks were then projected onto each specimen. To achieve this, atlases were generated with the 'createAtlas' function in 'Morpho' using each template, and these atlases were then used in the function 'placePatch' in 'Morpho' to project the surface semi-landmarks onto each specimen in a semiautomatic procedure (Fig. 2). Inflate values (inflate the semi-landmarks along the normals of the deformed atlas to ensure they stay on the surface of the target mesh) used in the 'placePatch' function were 0.3 and 0.05 for the dorsal and ventral surfaces, respectively, with the argument 'relax.patch' set to false. The semilandmarks were then slid using the function 'slider3D' from 'Morpho' (with three iterations and step size of 0.1) to minimize the bending energy of a thin-plate spline (Gunz et al., 2005;Dumont et al., 2016;Schlager, 2017), in order to generate dense and comparable 3D shape coordinates for large areas that lack clear suture lines or other clearly defined, homologous landmarks.
The final template for B. boulengeri contained 1626 points: 68 landmarks, 822 sliding semi-landmarks on 72 curves, and 736 sliding semi-landmarks for surfaces representing 17 regions (Fig. 3A). The final template for I. russeli contained 1609 points: 66 landmarks, 807 sliding semi-landmarks on 70 curves, and 736 sliding surface semi-landmarks (Fig.  3B). We were unable to find a set of parameters that correctly projected all surface semi-landmarks on every specimen, necessitating the removal of 53 and 140 surface semi-landmarks from the B. boulengeri and I. russeli datasets, respectively; 18 of these surface semi-landmarks occurred in both species, such that 175 surface semi-landmarks were removed from all specimens for interspecific comparisons (Supporting Information, Table S6), along with the aforementioned non-corresponding landmarks and semi-landmarks. The resulting datasets were used for all subsequent analyses unless otherwise stated.
The function 'mirrorfill' in the R package 'paleomorph' v.0.1.4 (Lucas & Goswami, 2017) was used to mirror the landmarks and semi-landmarks onto the left side of the cranium to improve the Procrustes alignment (Cardini, 2017). Generalized Procrustes analysis (Gower, 1975) was performed using the 'gpagen' function in the R package 'geomorph' v.3.0.5 (Adams & Otárola-Castillo, 2013;Adams et al., 2017), which removed all non-shape aspects of the data, specifically rotation, translation and isometric size (Rohlf & Slice, 1990;Dryden & Mardia, 1998). The mirrored points were removed before analyses in order to reduce dimensionality of the dataset. The progression of the three-dimensional (3D) sliding semi-landmark procedure for Boulengerula boulengeri specimens. A, 68 landmarks and 340 sliding semi-landmarks on 72 3D curves were manually digitized on one side of every specimen. B, identical landmarks and curve sliding semi-landmarks, along with 736 surface semi-landmarks were manually placed onto a template mesh for each species. C, this template was used in a semi-automatic procedure to project the surface semi-landmarks onto the cranium of every specimen. An identical procedure was followed for Idiocranium russeli specimens. Images for steps A and B were captured in IDAV Landmark Editor. Images for step C were captured in R. Some slight inconsistencies in the appearance of the pictured skull are attributable to differences in the field-of-view perspective angle within the different software packages. Specimen shown: BMNH 2002.809.

Modularity analyses
Two methods were used to identify patterns of modularity within the cranium. First, a maximum likelihood approach with landmark correlation matrices (Goswami & Finarelli, 2016) implemented in the 'subSampleEMMLi' function in the R package 'EMMLiv2' v.0.0.3 (https://github.com/hferg/EMMLiv2; Accessed December 2018) was used to compare support for alternative hypotheses of modular organization. Eleven models of modularity were defined based on proximity and current understanding of functional and developmental connectivity in the skull of caecilians and other vertebrates (Bardua et al., 2019) (Table 2). These included a model analogous to the six-module placental mammal modular pattern (Goswami, 2006a) adapted to caecilian cranial architecture [consisting of a facial/snout module, a palate module, a jaw/ cheek module, an occipital module, a vault module and a ventral os basale (to some extent analogous to the basisphenoid region in mammals and birds) module], models adapted from (Sherratt, 2011) and variations similar to that found by Bardua et al., (2019). The final set of models ranged from a model of maximal modularity (each of the 17 cranial regions as a separate module) to one of complete cranial integration. Analyses were performed on the full dataset of landmarks and semi-landmarks, and on only the Type I and Type II landmarks, as described above. Given that the mesethmoid was delineated with only two landmarks, this region was excluded from landmark-only analyses, but was included in all other analyses. Owing to the high ratio of data to specimen number in the full dataset and the potential  Table 1. consequences of this high dimensionality (Goswami & Polly, 2010), robustness of the results was evaluated by conducting analyses of modularity with EMMLi after performing a jack-knife subsampling procedure wherein the full dataset of landmarks and semilandmarks was randomly subsampled to 10% over 100 iterations. Robustness was evaluated further by repeating analyses after correcting the data for allometric effects. Allometry was removed from all datasets by performing a Procrustes ANOVA using the 'procD.lm' function in the R package 'geomorph' v.3.0.3 (Adams & Otárola-Castillo, 2013;Adams et al., 2016), with log centroid size ('the square root of the sum of squared distances between the centroid of an object and its landmarks'; Souter et al., 2010) as a factor, and extracting the residuals. Following Felice & Goswami (2018) and Bardua et al., (2019), we compared withinand between-region trait correlations in the bestsupported model found by EMMLi to determine whether any regions could reasonably be combined into multi-region modules. Pairs of regions in which the between-region correlation (ρ value) was stronger than, or within 0.1 of, either within-region correlation were merged if they were proximal or potentially had functional or developmental associations (based on comparisons with studies of other taxa; e.g. Goswami, 2006a;Maddin et al., 2016) (Table 3). The threshold of 0.1 was set based on observations of the pattern and distribution of within-and between-region correlations. This approach, examining the output of EMMLi analysis beyond simple identification of the best-supported model, allows for a posteriori consideration of models of modular organization that were not originally tested or hypothesized, which is crucial as analyses move away from well-studied systems.
Second, covariance ratio (CR) analysis (Adams, 2016) was applied to the complete, landmark-only and allometry-corrected datasets using the 'modularity. test' function in the R package 'geomorph' v.3.0.5. Covariance ratio analysis was also undertaken for each pair of regions that were merged into modules from EMMLi analysis, to determine whether CR analysis also supported the pooling of these regions into modules. Covariance ratio analysis quantifies the degree of modularity in a hypothesized model with the CR coefficient, a ratio of the overall between-module covariation relative to the overall within-module covariation (Adams, 2016), where a coefficient of one indicates the absence of modularity (i.e. the levels of between-module covariation equal those of withinmodule covariation). Covariance ratio coefficient values between zero and one indicate the presence of a more modular structure (i.e. between-module covariation is less than within-module covariation), whereas values greater than one describe structures with greater between-module covariation compared with within-module covariation. The significance of this coefficient is assessed by a permutation procedure (1000 iterations), randomly assigning landmarks into subsets of the same cardinality as the model being assessed to generate a null distribution for comparison.

Procrustes variance analyses
The 'morphol.disparity' function in 'geomorph' v.3.0.5 was used to quantify and compare the degree of cranial shape variation within and between each species. This function estimates the shape variance as Procrustes variance (Pv), using residuals of a linear model fit. Significance was assessed through a randomized residual permutation procedure using 1000 permutations. The variance of each individual cranial region and module was quantified and compared after correcting for differences in region and module size by dividing by the number of landmarks and semi-landmarks on each. Adjusted variances were compared with within-module correlation values to ascertain whether a relationship exists between the magnitude of within-module integration and module variance. This comparison was also performed on the allometry-corrected datasets.

EMMLi ('evaluating modularity with maximum likelihood') modularity analyses
For both species, the 17-module model, with each region acting as a separate module, was best supported by the complete datasets and when subsampling 10% of the data, both with and without correction for allometry. Likewise, with the landmark-only datasets the model with the greatest division of traits, in this case a 16-module model (owing to the exclusion of the mesethmoid), was the best supported for B. boulengeri both with and without allometric correction. For I. russeli, a 16-module and a 15-module model emerged as the best-supported models for the original landmarkonly dataset, with the 15-module model being the best supported after allometric correction. This 15-module model combines the jaw joint and quadrate regions into one module, with all other regions acting as separate modules.
Comparison of the between-and within-region estimated correlations that are output from the EMMLi analysis indicated several pairings with relatively strong between-region trait correlations relative to (≤ 0.1 different from) either of the withinregion correlations. After comparing between-and within-region correlations (Fig. 4Ai, Bi; Supporting Information, Tables S7 and S8) in both species, the occipital and the occipital condyle were merged into an 'occipital module', the dorsal and palatal nasopremaxilla merged to form a 'facial/snout module', and the jaw joint, squamosal and quadrate merged into a 'jaw/cheek module' (Fig. 5). This last grouping mirrors the results of the EMMLi analysis on the landmarkonly datasets for I. russeli, wherein the jaw joint and quadrate regions were combined. Although the network graphs for both species (Fig. 4Ai, Bi) show a strong correlation between the pterygoid and this 'jaw/cheek module', the between-region correlations far exceeded the 0.1 difference threshold; therefore, the pterygoid was not included in this module. Following these steps, only one module differed between the two species, with the lateral and interdental maxillopalatine regions covarying strongly and being merged into a 'maxillary module' in I. russeli but not in B. boulengeri. The final best-fit modularity patterns for B. boulengeri and I. russeli were considered to be a 13-module (frontal, parietal, facial/snout, lateral maxillopalatine, occipital, jaw/cheek, ventral os basale, vomer, interdental maxillopalatine, palatine shelf of maxillopalatine, pterygoid, stapes and mesethmoid) and 12-module (frontal, parietal, facial/snout, maxillary, occipital, jaw/cheek, ventral os basale, vomer, palatine shelf of maxillopalatine, pterygoid, stapes and mesethmoid) model, respectively. Subsampling the complete datasets to 10% of landmarks/semi-landmarks over 100 iterations and applying the 0.1 difference threshold returned identical patterns of modularity to those observed for the full dataset ( Fig. 4A ii, Bii; Supporting Information, Tables S9 and S10). Correcting for allometry for both the complete and 10% subsample datasets also revealed near identical results (Fig. 4A iii, Biii; Supporting Information, Fig. S1Ai, Bi; Tables S11-14). Results of EMMLi analysis of the landmarkonly datasets for each species yielded patterns that are generally similar to the complete datasets in terms of between-region magnitudes of integration ( Fig. 4A iv, Biv); however, within-region correlations were generally lower and overlapped more with between-regions correlations, suggesting a much more integrated structure (Supporting Information, Tables S15 and S16). Although correcting the complete and 10% subsampled datasets for allometry had almost no effect, allometric correction had a large effect on the landmark-only dataset, supporting a more modular hypothesis similar to that suggested by the other datasets (Supporting Information, Fig. S1Aii, Bii; Tables S17 and S18). Allometry had a significant (P = 0.001) and similar influence on cranial shape variation in the complete and landmark-only datasets for B. boulengeri (R 2 = 0.285 and 0.271, respectively) and for I. russeli (R 2 = 0.250 and 0.235, respectively).
In both species, the stapes, mesethmoid and pterygoid modules had the highest (≥ 0.75) withinregion correlation, and the ventral os basale module had the lowest (≤ 0.46) within-region correlation. There appeared to be a general tendency for smaller-sized regions to exhibit higher within-region correlations relative to larger modules, although surface semilandmark density is similar across regions.

Covariance ratio analyses
The CR coefficients of the 17-module model of B. boulengeri and of I. russeli were both significant (P = 0.001) at 0.725 and 0.621, respectively, supporting a modular structure in the crania of both species. After correcting for allometry, the CR coefficients were still significant (P = 0.001), with values of 0.563 and 0.519 for B. boulengeri and I. russeli, respectively. The CR coefficients of the landmark-only dataset, excluding the mesethmoid as in the EMMLi landmark-only analyses, were high but significant at 0.988 (P = 0.002) and 0.851 (P = 0.001) for B. boulengeri and I. russeli, respectively. As in the EMMLi analyses, CR coefficients of the landmark-only analyses supported greater modularity after allometric correction, with values of 0.824 and 0.738 (P = 0.001 for both) for B. boulengeri and I. russeli, respectively.
For the complete datasets, regions that were pooled into modules after EMMLi analysis were then investigated using pairwise CR analyses, which partly supported these proposed modules. The CR values for most pairs of regions were very close to one (i.e. suggesting pooling of the regions), although all values were still found to be significantly modular (P = 0.001), i.e. suggesting that the regions should remain separate. We considered these high CR values to indicate that these pairings are in fact highly integrated (albeit not as much as the individual regions themselves), and hence support the pooling of these regions into modules. The only pairing with a CR value considerably lower than one was the jaw joint and squamosal pairing in I. russeli (0.632), which does not support the pooling of these two regions. However, the integration among the other region pairs in that module supported the combining of the jaw joint and squamosal with the quadrate. Therefore, these results generally supported the combining of these regions into the 13-module and 12-module models recovered from the EMMLi analysis for B. boulengeri and I. russeli, respectively. Covariance ratio analysis also suggested some additional groupings that were not identified in the EMMLi analyses. Although these were all supported as significantly modular (P = 0.001), CR coefficient values approached one between the parietal and occipital regions (0.888; P = 0.001) and the jaw joint and pterygoid regions (0.875; P = 0.001) in B. boulengeri, and between the frontal and parietal regions (0.873; P = 0.001) and the interdental maxillopalatine and palatine shelf regions (0.873; P = 0.001) in I. russeli. In B. boulengeri, high covariation was also found between the lateral and interdental maxillopalatine regions (0.881; P = 0.001), a relationship that was identified from the , showing the results of EMMLi analysis: (i) sampled at 100%; (ii) mean from 100 iterations subsampling the full dataset to 10%; (iii) sampled at 100% with allometric correction; and (iv) using only landmarks. The size of the circle for each unit represents the magnitude of withinmodule correlations, and the line thickness is proportional to the magnitude of the between-module correlation. The layout of each graph corresponds approximately to a cranium in right lateral view, with the anterior extreme to the right of the image. Abbreviations are as in Table 1. EMMLi analysis for I. russeli but not for B. boulengeri. Pairwise CR coefficient results are presented in Table  4. Correction of the complete dataset for allometry resulted in slightly lower pairwise coefficient values throughout the dataset; however, the regions combined into modules from the EMMLi analyses still showed the highest values, and the same overall pattern of modularity was still supported (Supporting Information, Table S19).
As in the EMMLi analyses, the results of the CR analyses of the landmark-only datasets supported a much more integrated structure than the results of the complete datasets, with most pairwise coefficient values ≥ 0.85 (Supporting Information, Table S20). However, allometric correction of the landmarkonly dataset returned a more modular structure, more similar to that of the complete datasets, with the regions combined into modules after EMMLi analysis generally showing the highest CR coefficient values (> 0.9), consistent with the suggestion that these pairings are highly integrated (Supporting Information, Table S21). Other pairings that had especially high CR coefficient values (> 1.0) in both species were the frontal and parietal (B. boulengeri, 1.192; I. russeli, 1.296), the vomer and interdental maxillopalatine (B. boulengeri, 1.056; I. russeli, 1.069), and the vomer and palatine shelf of the maxillopalatine (B. boulengeri, 1.276; I. russeli, 1.087). Other very high CR coefficient values were found in one species and not the other; for example, the occipital condyle and squamosal (1.106) and the occipital condyle and quadrate (1.086) in I. russeli, and in B. boulengeri the frontal and vomer (1.108), the parietal and occipital (1.198), the parietal and vomer (1.215), the occipital  Table 1. and vomer (1.095), the ventral os basale and vomer (1.149) and the palatal nasopremaxilla and vomer (1.219).
For each species, the Procrustes variance of each module was adjusted for the number of landmarks representing that module (Table 5), after which the variances of all multi-region modules were greater than those of any single-region modules in B. boulengeri, with the vomer having the lowest variance overall (Fig. 6A). The vomer also had the lowest variance in I. russeli, but the multi-region modules did not exhibit the highest variances; instead, the ventral os basale and parietal modules had the highest variances (Fig.  6B). Comparison of the variances of each module with their respective within-module correlations showed no obvious pattern, with the results of a linear regression showing no significant relationship in either B. boulengeri (multiple R 2 = 0.0329, adjusted R 2 = −0.0550; P = 0.553) or I. russeli (multiple R 2 = 0.0586, adjusted R 2 = −0.0356; P = 0.449) (Fig. 7). This was also the case when using the allometrycorrected datasets for each species; B. boulengeri (multiple R 2 = 0.0223, adjusted R 2 = −0.0666; P = 0.627) and I. russeli (multiple R 2 = 0.00175, adjusted R 2 = −0.0981; P = 0.897) (Supporting Information, Fig. S2).

DISCUSSION cranial modularity in caecilianS
As revealed by our analyses, cranial modularity is very similar in B. boulengeri and I. russeli, supporting patterns of 13 modules and 12 modules, respectively. This high degree of modularity is therefore plausibly conserved (at least among non-scolecomorphid teresomatan caecilians). These modularity patterns are strikingly more complex than the six-module structure identified in therian mammal crania (Goswami, 2006a) and the seven-module structure of avian crania (Felice & Goswami, 2018), discussed further below. The patterns are very similar to the highly partitioned ten-module structure recovered by Bardua et al., (2019) in their clade-wide caecilian study, with analyses of both species here and the interspecific analysis (Bardua et al., 2019) all supporting separated pterygoid, stapes, frontal and parietal modules. All three analyses also support an 'occipital module' combining the occipital and occipital condyle regions, and a 'jaw/cheek module' merging the quadrate, jaw joint and squamosal regions. There is some similarity in the maxillopalatine region, with Bardua et al.,'s (2019) interspecific analysis supporting merger of the lateral, palatine shelf and interdental regions of the maxillopalatine into one module, and our analysis of I. russeli supporting merger of the lateral and interdental regions into one module. However, there is no support for inclusion of the palatine shelf in this module in I. russeli, and no support for any merging of these regions in the analysis of B. boulengeri. The main differences between the pattern recovered by Bardua et al., (2019) and those for B. boulengeri and I. russeli reported here is that the intraspecific analyses supported merger of the dorsal and palatal nasopremaxilla regions, whereas the interspecific study supported these as separate modules; conversely, the interspecific study supported merger of the ventral os basale and vomer regions into one module, with the analyses of the individual species supporting these as separate modules. The only other difference between the pattern found by the cladewide study and those we recovered here is the addition of a separate mesethmoid module, because Bardua et al. (2019) merged this region with the frontal region a priori owing to its variable presence across the clade. Thus, the patterns of trait correlations between cranial regions were highly congruent across both intra-and interspecific studies, suggesting the conservation of modularity in caecilians at both static and evolutionary levels.
Neither EMMLi analysis of the complete and 10% subsampled datasets nor CR analysis of the complete datasets supported the combination of elements into an integrated palate or vault module in either B. boulengeri or I. russeli, although some of the relevant elements had CR values that were significant, but relatively close to one (0.87-0.88), as described above. However, the landmark-only analyses did support some combining of elements within palate and vault regions, even after allometric correction. The lack of consistent support for a highly integrated vault module is perhaps surprising given its role in support and protection of the brain (Goswami, 2006a), which might be expected to be important in species that practise head-first burrowing, particularly in the putatively more dedicated burrower, B. boulengeri. Comparison of the relationships between the potential vault module components (the frontal, parietal and mesethmoid regions) in the two species gave very similar results, with no consistently stronger support for combination of these regions into a single module for either species. In both species, there was support for a ventral os basale ('basisphenoid') module, with a role in braincase support (Goswami, 2006a), and an occipital module constituting the occipital and occipital condyle regions, which typically consist of separate ossifications that fuse early in development in most tetrapods and have a shared role in skullaxial skeleton attachment. Both species showed support for a jaw/cheek module constituting the jaw joint, squamosal and quadrate, which is likely to be critical in feeding and have a pivotal role in the dual jaw-closing mechanism unique to caecilians (Nussbaum, 1983). There was also some support for a facial/snout module through grouping of the dorsal and palatal nasopremaxilla regions. This last grouping was perhaps expected because these regions are part of the same bone; however, strong covariation between regions on the same bone was not consistently found, with the three maxillopalatine regions remaining separate in B. boulengeri and only the interdental and lateral surfaces grouping in I. russeli. This finding that single cranial bones can exhibit modular structures highlights the value of splitting individual bones into multiple regions, particularly where those bones may contribute to multiple functional structures.
The crania of these caecilian species were more modular than previous analyses that were based solely on landmarks have suggested (Sherratt, 2011). This difference is almost certainly a result of the differences in data collection approach and in the amount of morphology captured and analysed by landmark-only studies vs. the dense surface sliding semi-landmark analysis applied here. Surfacebased methods allow the complex morphology of the cranium to be more fully represented than Figure 6. Network graphs from the EMMLi analysis of the 13-module and 12-module models for Boulengerula boulengeri (A) and Idiocranium russeli (B). The size of the circle for each unit represents the magnitude of withinmodule correlations, and the line thickness is proportional to the magnitude of the between-module correlation. The layout approximately corresponds to a cranium in right lateral view. Modules are graded high (red) to low (yellow) according to their Procrustes variance values. Abbreviations are as in Table 1. traditional landmark-only approaches. The placement of landmarks is limited largely to sutures, which by definition lie at boundaries of potential modules rather than within modules. Therefore, sampling more of the morphology, in particular the surface structure and within-element variation, with semilandmarks might be expected to better reflect the biological covariance distribution of the cranium. A similar finding was reported by Parr et al. (2016), who recovered ten modules in dingo crania, a higher level of modularity than had previously been found in carnivorans (Goswami, 2006b). They attributed the additional divisions to their finer sampling of cranial morphology using a surface-based method, compared with previous carnivoran studies using only Type I and II landmarks (Parr et al., 2016). The identification of a stronger pattern of modularity in caecilians than in birds might also reflect (at least to some extent) methodological differences, despite the use of dense surface semi-landmarks in both studies, but also biological differences. Bird skulls are highly fused, and thus equally modular patterns to those analysed here cannot be tested, because fewer cranial regions can be demarcated readily in birds than in caecilians and other vertebrates. Sliding semi-landmarks are inherently dependent on each other within the curve or surface patch, and therefore consideration should be given to how the recovered patterns of integration and modularity compare to those supported by landmarkonly approaches. However, the high degree of fusion in bird skulls, which makes demarcation of more regions impossible, might also indicate that bird skulls are more integrated and less modular than those of caecilians.
Both EMMLi and CR modularity analyses of the landmark-only datasets indicated a more integrated structure in the caecilian crania than was revealed by the complete datasets. However, allometric correction had a substantial effect on the landmark-only analyses, after which results were more congruent with analyses of the full landmark and semi-landmark dataset. Allometry is a highly integrating factor (Klingenberg, 2013), which might explain why its removal results in a more modular structure. Removal of allometry from the complete and 10% subsampled datasets also decreased the magnitude of integration slightly when analysed using EMMLi, and when the complete dataset was analysed using CR; however, this effect was negligible compared with the effect on the landmark-only datasets. This is interesting given that the degree of influence of allometry on the complete and landmark-only datasets was found to be comparable, being 28.5 and 27.1%, respectively, in B. boulengeri and 25.0 and 23.5%, respectively, in I. russeli. As noted above, landmark-only analyses are likely to support higher integration across regions than analyses involving surface semi-landmarks owing to the weighting of landmarks on sutures between regions, and thus instances where the two approaches produce discordant results might accurately reflect the nuanced and complex hierarchy of integration and modularity in the vertebrate skull. Given the finding that removal of allometry has a relatively larger effect on landmark-only analyses, it might be prudent to include this correction when drawing comparisons between analyses using semi-landmarks and those using only landmarks.

procruSteS variance and integration
The finding that I. russeli displays slightly but significantly greater variation than B. boulengeri might suggest that these species have been subject to different selection pressures. Boulengerula boulengeri is likely to be a more dedicated burrower than I. russeli, with greater functional constraint on cranial architecture owing to the demands of head-first  Table 1. burrowing. The putative miniaturization of I. russeli (Wake, 1986;Theska et al., 2018) and its overall smaller size relative to B. boulengeri might also be connected causally to its greater variance. Specifically, reduction of particular elements (e.g. frontals and nasopremaxillae) is indicative of a selective regime in which greater variation in these elements might be anticipated. However, examination of the variance of individual modules showed the frontal modules to have a similar variance in both B. boulengeri and I. russeli, and the facial/snout module (the dorsal and palatal nasopremaxilla) to have a relatively greater variance in B. boulengeri. This result was particularly unexpected given the suggested differences in burrowing behaviour between the species. The more dedicated B. boulengeri might have been expected to display less variation in the modules most relevant to head-first burrowing, i.e. the facial/snout module, particularly given previous findings of extreme conservatism in snout shape in other head-first burrowers (e.g. Hipsley et al., 2016).
Interpretation of differences in the variance of the facial/snout module is further complicated by its involvement in the maternal dermatophagy of, and the associated presence of a specialized and ephemeral vernal dentition in, hatchlings of B. boulengeri (San Mauro et al., 2014). We might have expected the dentigerous modules involved in dermatophagy to display higher variance owing to the influence of this additional, specialized life-history stage on the phenotypic variation of these modules. Comparison of the dentigerous modules between the species shows only the facial/snout module (the dorsal and palatal nasopremaxilla regions) to be in keeping with our prediction, displaying relatively high variance in B. boulengeri and relatively low variance in I. russeli, whereas the maxillopalatine modules have relatively low variance values in both species.
We find a weak inverse relationship between module size and within-module correlation (= integration) in both species, with smaller modules, such as the mesethmoid and pterygoid, showing the highest integration values, and larger modules, such as the ventral os basale, showing the lowest. Although this pattern could suggest a statistical bias, for example if smaller modules were more densely sampled in the surface-based method, we tried to mitigate this effect by sampling the cranial regions proportionately to their size, i.e. with an even density across the cranium, although this can be difficult to achieve in practice. A similar pattern was also observed when analyses were limited to landmarks, supporting a biological explanation rather than a statistical one. Smaller modules perhaps show greater integration because they may be more likely to originate from a single developmental source or serve fewer functions in the cranium, or both, which may promote higher integration in these regions than across larger elements or structures, which may be derived from multiple tissue types or be involved in multiple competing functions. Multi-region modules generally had higher variance values than single-region modules, probably owing to their larger size and greater complexity in development and function than smaller modules.
The vomer module had the lowest variance in both B. boulengeri and I. russeli, indicating that this module might be under greater constraint in terms of available phenotypic variation or owing to greater stabilizing selection. There is growing evidence, including findings from recent work on I. russeli (Theska et al., 2018), that the vomer is among the first cranial bones to ossify in oviparous indotyphlids (Müller et al., 2005;Müller, 2006), whereas it has been found to ossify later in the viviparous Dermophis mexicanus (Wake & Hanken, 1982). This lends some support to the idea that differences in ossification sequences can be, at least in part, attributed to differences in life histories. Given that I. russeli and B. boulengeri are both oviparous direct developers, it may be the case that the vomer is also among the earliest bones to ossify in B. boulengeri. Thus, its low variance in both species might be explained by a reduced opportunity for variation as a result of early ossification.
Comparison of module variances and magnitude of integration showed no overall significant relationship, with some highly integrated modules exhibiting high variances, whereas others exhibited lower variances. Goswami & Polly (2010) predicted that under the 'facilitation hypothesis' (shape change is facilitated by the presence of modules) highly integrated modules would exhibit high variances and weakly integrated modules would exhibit low variances, whereas under the 'constraint hypothesis' (modular organization restricts shape change) highly integrated modules would exhibit low variances and vice versa. The lack of a clear pattern in our results might suggest that strong integration facilitates evolution in some modules and constrains it in others, or that these effects are minimal, at least at this scale of analysis and in these taxa. Given that most current understanding of phenotypic modularity comes from intensive study of a few clades, especially mammals, the identification of a conserved but highly modular pattern of cranial organization in these two species of caecilians provides important new data for further analysis of the evolution and macroevolutionary significance of phenotypic modularity in the vertebrate skull.

CONFLICTS OF INTEREST
The authors declare no conflicts of interest.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Table S1. List of Boulengerula boulengeri (top) and Idiocranium russeli (bottom) specimens included in study. Specimens with a BMNH and an MW prefix are accessioned and yet-to-be accessioned specimens in the Natural History Museum, London, respectively. Specimens with a UCMZ prefix are in the University of Cambridge Museum of Zoology. Abbreviation: FR, forest reserve. Table S2. Scan settings for each specimen as taken from text files within scanning folders. All scans were performed with 360 o rotation. Table S3. Definition of the Type I and II landmarks assigned to each specimen. Abbreviation: LM, landmark index. Table S4. Description of curves applied to each specimen. Abbreviation: LM, landmark index. Table S5. Description of surface semi-landmarks applied to each specimen. Table S6. Original surface semi-landmarks assigned to each region and the surface semi-landmarks removed before analyses. Landmark indices have been adjusted to be consistent in both species. The presence of '-' indicates that no surface semi-landmarks were removed; 'NA' indicates that this region is not applicable to this species. Table S7. Within-and between-region ρ values from EMMLi analysis of the complete dataset of Boulengerula boulengeri, without allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S8. Within-and between-region ρ values from EMMLi analysis of the complete dataset of Idiocranium russeli, without allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S9. Average within-and between-region ρ values from EMMLi analysis of a 10% subsample of the complete dataset of Boulengerula boulengeri over 100 iterations, without allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S10. Average within-and between-region ρ values from EMMLi analysis of a 10% subsample of the complete dataset of Idiocranium russeli over 100 iterations, without allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S11. Within-and between-region ρ values from EMMLi analysis of the complete dataset of Boulengerula boulengeri, with allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S12. Within-and between-region ρ values from EMMLi analysis of the complete dataset of Idiocranium russeli, with allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S13. Average within-and between-region ρ values from EMMLi analysis of a 10% subsample of the complete dataset of Boulengerula boulengeri over 100 iterations, with allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S14. Average within-and between-region ρ values from EMMLi analysis of a 10% subsample of the complete dataset of Idiocranium russeli over 100 iterations, with allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S15. Within-and between-region ρ values from EMMLi analysis of the landmarks-only dataset of Boulengerula boulengeri, without allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S16. Within-and between-region ρ values from EMMLi analysis of the landmarks-only dataset of Idiocranium russeli, without allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S17. Within-and between-region ρ values from EMMLi analysis of the landmarks-only dataset of Boulengerula boulengeri, with allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S18. Within-and between-region ρ values from EMMLi analysis of the landmarks-only dataset of Idiocranium russeli, with allometric correction; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses. Table S19. Results of applying the covariance ratio modularity analysis to the complete datasets of Boulengerula boulengeri and Idiocranium russeli with allometric correction. Bold text highlights the values ≥ 0.85; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses of the complete datasets. The upper triangle represents the results from B. boulengeri, and the lower represents results from I. russeli .  Table S20. Results of applying the covariance ratio modularity analysis to the landmark-only datasets of Boulengerula boulengeri and Idiocranium russeli. Bold text highlights the values ≥ 0.85; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses of the complete datasets. The upper triangle represents the results from B. boulengeri, and the lower represents results from I. russeli .  Table S21. Results of applying the covariance ratio modularity analysis to the landmark-only datasets of Boulengerula boulengeri and Idiocranium russeli with allometric correction. Bold text highlights the values ≥ 0.85; values marked with an asterisk indicate that these region pairings constituted the modules recovered from the EMMLi analyses of the complete datasets. The upper triangle represents the results from B. boulengeri, and the lower represents results from I. russeli. Figure S1. Network graphs of Boulengerula boulengeri (A) and Idiocranium russeli (B), showing the results of EMMLi analysis: (i) mean from 100 iterations subsampling the full dataset to 10% with allometric correction; and (ii) sampled using only landmarks with allometric correction. The size of the circle for each unit represents the magnitude of withinmodule correlations, and the line thickness is proportional to the magnitude of the between-module correlation. The layout corresponds approximately to a cranium in right lateral view. Abbreviations are as in Table 1. Figure S2. Linear regression models of the relationship between individual module Procrustes variance and within-module correlation after correction for landmark number in Boulengerula boulengeri (A) and Idiocranium russeli (B) for the allometry-corrected datasets. In both species, there was no significant relationship between the individual module variance and internal correlation values. Abbreviations are as in Table 1. File S1. B.boulengeri landmark data.csv. File S2. I.russeli landmark data.csv. SHARED DATA 3D models are available for download on www.phenome10k.org