Common Patterns of Skull Bone Fusion and Their Potential to Discriminate Different Ontogenetic Stages in Extant Birds

The degree of sutural closure between bones generally allows for the classification of skeleton maturity in tetrapods. In mammals, the sutural closure of skull bones was previously used as proxy to evaluate the ontogenetic stage of single individuals. However, due to temporal variation, this process can be only applied among mammalian subclades, but not for all mammals in general. In contrast, the process of sutural closures in bird skulls could be a more reliable ontogenetic proxy for this clade as adult birds commonly show a generally high degree of bone fusion. To test this, we studied the process of sutural closure in ontogenetic series of 18 extant bird species regarding the presence of an ontogenetic signal and compared the results with changes in skull size and proportions. Univariate analyses indicate that bone fusion happens faster in altricial than in precocial birds. However, the use of PCoA and multivariate regressions reveal that the skull bone fusion follows a common pattern among birds and thus can be used as proxy to identify different ontogenetic stages. In general, the process of sutural closure spreads from posterior to anterior and from ventral to dorsal. In contrast, skull measurements reflect rather interspecific allometry than ontogeny. The used of bone fusion as proxy will help to better identify and compare different stages of maturation in birds, including historical material from osteological collections.

As adult birds show in general an extreme degree of bone fusion in the skull (Jollie, 1957;Plateau and Foth, 2020), the process of sutural closure may provide another potential proxy for the identification of postnatal ontogenetic stages (Rager et al., 2013;Bailleul et al., 2016). Several case studies already described patterns of sutural closure in a handful of birds, e.g., Gallus gallus (Jollie, 1957), Macronectes giganteus (Piro and Hospitaleche, 2019), Aptenodytes forsteri (Sosa and Hospitaleche, 2018), Casuarius casuarius (Green and Gignac, 2020), Myiopsitta monachus (Carril et al., 2020), and Dromaius novaehollandiae (Bailleul et al., 2016). According to Bailleul et al. (2016), the progressive nature of the sutural closure, allows identifying different ontogenetic stages. In addition, in their macroevolutionary study, Plateau and Foth (2020) found a common pattern of skull bone fusion among extant birds, but taking only two ontogenetic stages per species into account. Nevertheless, the results of the latter two studies indicate a strong correlation between the degree of skull bone fusion and the ontogenetic stage.
To further explore this process, we compared the sutural closure of skull bones among the ontogenetic series of 18 extant bird species. Using univariate and multivariate analyses, we identified major patterns of sutural closure common in birds and compared them to morphometric measurements regarding the presence of an ontogenetic signal. Because the developmental mode (precocial vs. altricial) of birds relates to noticeable differences in growth rate, egg size and hatchling size, we tested if a similar pattern can be found in skull bone fusion.

Taxon Sampling and Data Collection
For a total number of 18 extant bird species (Supplementary File 1), we sampled multiple specimens of at least three individuals with different ontogenetic stages, which were determined based on the degree of fusion in the skull. The species were grouped according to their (a) taxonomy [in Palaeognathae (number of specimens: n = 7), Inopives (n = 79), Strisores (n = 4), Colombaves (n = 7) and Aeqorlitornithes (n = 34) following Prum et al., 2015] and (b) developmental modes [precocial (n = 29) and altricial (n = 102) summarizing the different subcategories from Botelho and Faunes, 2015].
Following Bailleul et al. (2016), a matrix documenting the contacts between skull bones was created for each bird specimen, identifying ontogenetic sutural changes based on the following character states: All neighboring bones with a minimal physical suture, but a distinct opening along the suture line (0); neighboring bones that are partially in physical contact, but with minor openings along the suture line (1); neighboring bones that are in physical contact having a continuous suture (2); neighboring bones are completely fused to each other (3) (Figures 1A,B). Due to the interspecific variation of the pterygoid-palatine complex in modern birds (Bellairs and Jenkin, 1960;Plateau and Foth, 2021), only the contacts between the palatine and premaxilla and maxilla were taken into account, while contacts that include the pterygoid and vomer were ignored. In total, we identified 35 sutural characters changing during ontogeny (Table 1 and Supplementary Table 1).
To compare the pattern of bone fusion between different species, we calculated the average of the sutural closure (ASC) from the ontogenetic matrices of each bird specimen by calculating the sum of character states divided by the total number of sutural characters (Wilson and Sánchez-Villagra, 2009;Bailleul et al., 2016). Juvenile birds will have a mean value closer to 0, while the adult mean value will be closer to 3. Based on the ASC values, which represent a continuous range from 0.84 to 2.95, we established a third grouping defining three equally ranged ontogenetic stages [Stage 1 (n = 18), Stage 2 (n = 29) and Stage 3 (n = 84)].
In addition, we created a second matrix of 14 linear skull measurement estimating morphometric variation through ontogeny (Supplementary Table 2 and Supplementary  Figure 1). Due to the large interspecific variation in the skull size, ranging from 30.39 mm (Apus apus) to 247.2 mm (Ciconia ciconia) in length, all measurements were standardized as follows, minimizing the effect of size. For each specimen, all measurements were divided through their geometric mean (Mosimann, 1970;Claude, 2008), which is the nth root of the product of n measurements and used as a general proxy for size (Claude, 2008). Afterward, the size corrected measurements were log e -transformed. To compare ontogenetic size variation between different bird species, the geometric mean of each specimen of a species was expressed as the ratio of the geometric mean from the largest specimen. This size proxy is herein called relative size percentage (RSP).

Univariate Analyses of Bone Fusion
In order to test the relationship between skeletal maturation and size for each bird species, we expressed the values of the average of the sutural closure as ratio from the most mature specimen of each species (RASC) and applied an ordinary least squares regression (OLS) between RASC and relative size percentage (RSP), using the lm function in R package stats 4.0.3 (R Development Core Team, 2011). The resulting slopes represent the rate of fusion in relation to size. In the next step, the slopes of precocial and altricial birds were compared using a non-parametric Mann-Whitney test, which estimates  Bailleul et al. (2016): Neighboring bones with a minimal physical suture, but a distinct opening along the suture line (0); neighboring bones with partial physical contact, but with minor openings along the suture line (1); neighboring bones that are in physical contact having a continuous suture (2); neighboring bones are completely fused to each other (3). All specimens not in scale.
whether two independent samples are taken from populations with equal medians (Hammer and Harper, 2006). This procedure was performed in the program PAST 4.03 (Hammer et al., 2001). In addition, we compared the relationship between RASC and RSP of all specimens with OLS, considering the precocial and altricial mode. Equality of regression slopes between both groups were tested with an ANCOVA using the lstrend function of the R package Emmeans 1.5.4 (Russell, 2021).
To test how the rates of fusion change through evolution, we mapped the slopes as continuous characters on to phylogenetic tree generated from http://birdtree.org/ (Jetz et al., 2012), using the Ericson et al. (2006) topology (Supplementary File 3). For that, we downloaded 9.999 phylogenetic trees and averaged the nodal ages of the single trees using the consensus.edges function of the R package phytools 0.7-70 (Revell, 2012). The slopes were traced onto the phylogeny, using the ace function of the R package ape 5.4-1 (Paradis, 2012), which uses a maximumlikelihood (ML) model under Brownian motion with a constant rate of diffusion. The phylogenetic strength of the fusion rate was estimated, using the phylosig function in phytools based on lambda optimization (Pagel, 1999).
To test if the fusion rate of the skull during ontogeny is dependent from the growth physiology of the birds, we obtained the growth rates of 16 bird species from Starck and Ricklefs (1998a), Goonewardene et al. (2003) and Cooper (2005). The relationship between fusion rate and growth rate was tested with help of Phylogenetic Independent Contrasts (PIC) (Felsenstein, 1985), using the pic function from the Ape package 5.4-1 (Paradis, 2012) and phylogenetic Generalized Least Squares (pGLS) (Adams, 2014), using the gls function (in combination with a phylogenetic tree and corBrownian argument) from the R package nlme 3.1-152 (Pinheiro et al., 2021).

Multivariate Analyses of Bone Fusion and Skull Measurements
In order to visualize the ontogenetic shape variability of bird skulls, we performed a principal coordinate analysis (PCoA) for two datasets: (a) all ontogenetic characters (including ASC characters) and (b) all size-corrected measurements (Dataset 2). Each dataset was transformed into a distance matrix (Wills et al., 1994), using the Gower coefficient (Gower, 1971). Following Caillez (1983), the PCoA was run with a transformation exponent c of 2, reducing the "horseshoe" effect (Podani and Miklos, 2002). Both steps were done, combining the vegdist (Vegan package 2.5-6; Oksanen et al., 2020) with the pcoa function of the Ape package in R. Finally, the number of relevant principal coordinates from both analyses was estimated using

Spenial-Prearticular 35
Names of the sutures correspond to the two bones in contact.
Numbering of the sutures correspond to the order in dataset 1.
the broken-stick-method (Frontier, 1976;De Vita, 1979;Jackson, 1993). To test the statistical degree of overlap between the different groupings, we used the most relevant PCOs of each PCoA (see above) to perform a permutational multivariate analysis of variance (perMANOVA) and a linear discriminant analysis (LDA). The perMANOVA is a non-parametric test that explores for significant differences in the distribution of groups in morphospace based on permutation (Anderson, 2001). The analysis was done in PAST, using a Euclidean distance matrix and 9,999 permutations. The p-values were Bonferroni corrected afterward, which multiplies the original values with the number of comparisons (Hammer and Harper, 2006). The LDA reduces a multivariate data set to a smaller set of dimensions by maximizing the separation between two or more groups (Fisher, 1936;MacLachlan, 2004). The method was performed using the lda function from the R package MASS (Ripley et al., 2013), with and without leave-one-out cross-validation, testing the accuracy of the predictions and prevent overfitting. Based on the PCoA results we estimated the size disparity of each ontogenetic stage in R, using the sum of variances as metric. All calculations are based on the first three PCO scores for Dataset 1 and the first six PCOs for Dataset 2, which summarize significant variation within the two datasets based on the broken stick method (Frontier, 1976;De Vita, 1979;Jackson, 1993). Differences in disparity between the ontogenetic groups in both datasets were tested for significance applying a permutation test with 10,000 replicates, using a modified R script of Brusatte et al. (2014), which takes sample size difference between the different groups into account ; see also Foth and Joyce, 2016). This script estimates whether a certain group had more or less total disparity than the other by keeping sample size constant and by shuffling taxa randomly between the groups being compared.
In addition, all relevant PCOs were transformed into regression scores using the RegScore function of the R package Morpho v2.8. (Schlager et al., 2020). The regression scores were tested against RASC and RSP with help of OLS regression. Using ANCOVA (see above), we tested if the slopes of precocial and altricial birds differ from each other or not.
As stated above, the ontogenetic stages are defined by the ASC using all bone contacts of the skull. However, to identify the ontogenetic stages of different individuals, it would be useful if one could define the stages based on a subset of contacts that are easy to evaluate. According to Green and Gignac (2020), the skeletal maturity of palaeognathous birds can be identified by the degree of ossification in the orbital region between the interorbital septum, mesethmoid, frontals, laterosphenoids, and basiparasphenoid complex, which represent only eight out of 35 contacts from Dataset 1. Thus, we repeated the regression analysis of RASC against RSP as well as the PCoA, perMANOVA and LDA with the scheme of Green and Gignac (2020) for all birds sampled and compared the outcome with the results of original analyses of Dataset 1 (Supplementary Table 4).

Analysis of Ontogenetic Fusion Patterns
In order to analyze the sequence of ontogenetic fusions in bird skulls, we ordered all specimens according to their ontogenetic groups based on the ASC and applied them to a hypothetical ontogram with three ontogenetic stages. Dataset 1 and the ontogram were loaded into the program Mesquite 3.51 (Maddison and Maddison, 2019) and the sequence of ontogenetic changes were analyzed by tracing the character history.

Relative Average of the Sutural Closure
For each species, the slopes of the OLS regressions of the RASC against RSP shows a positive correlation. Due to small sample sizes, the p-values of some species are not significant, but their R 2 scores higher than 0.57 (or 0.58 if both parameters are log e -transformed), still indicating a correlation (Supplementary Figure 2). The majority of birds show positive allometry, meaning that the fusion of bones happens faster than the growth of the skull. However, Struthio camelus has an almost isometric fusion pattern (slope: 0.933; p-value: 0.164), while it is negative allometric for Dromaius novaehollandiae (slope: 0.893; p-value: 0.408).
Based on the ANCOVA, the slopes do not indicate any taxonomic distinction between most bird clades, except for Inopinaves (slope: 2.17; p-value < 0.001), which are significantly different from Aequorlitornithes (slope: 1.41; p-value < 0.001) and Palaeognathae (slope: 0.92; p-value: 0.024) ( Supplementary  Table 9). However, the insignificance found between the other taxonomic group could be related to their small sample sizes (see above), and do not necessarily reflect a uniform development of this trait. In contrast, the slopes of precocial birds (slope: 1.18; p-value < 0.001) are significantly smaller than those of altricial ones (slope: 2.04; p-value < 0.001) (Figures 2A,B). The latter result is further supported by the Mann-Whitney test (Uvalue: 0.00; z-value: 2.92; p-value: < 0.05). The single slopes, which represent the rate of fusion, do not contain a phylogenetic signal (Pagel's λ: 0.97, p-value: 0.128) (Figures 2C,D). While Palaeognathae possess a slower rate than Neognathae, it is shared more randomly in the latter group (Supplementary Figure 3).

PCOA and Multivariate Regression of Fusion Patterns
For the dataset containing the ontogenetic sutural characters (Dataset 1), the first three relevant axes of the PCoA explain   (Figure 3A). Based on the perMANOVA all groups overlap with each other, except for Palaeognathae and Inopinaves. This result is supported by the LDA (total error of correct identification: 35.1%), which could only identify species from the Inopinaves with a high certainty (error of correct identification: 8.0%) (Figure 3A). Applying the ASC grouping on the PCoA results, we found the three ontogenetic stages (1, 2, and 3) separated from each other along PCO1, while they overlap with each other along PCO2 (Figure 3B). While stage 1 (Sum of Variances: 0.028) and stage 2 (Sum of Variances: 0.020) do not significantly differ from each other in size (p-value: 0.119), stage 3 (Sum of Variances: 0.006) is significantly smaller than stage 2 (p-value: 0.004). All three stages are found to be significantly different in the perMANOVA and the LDA (total error of correct identification: 0.8%) (Figure 3B). For the different developmental modes (using five groups), all morphospaces show a strong overlap along PCO1 and PCO2, in which their distribution along PCO2 follows the gradation of developmental spectrum from precocial to super-altricial to a certain degree ( Figure 3C). Based on the perMANOVA all developmental modes overlap with each other, except for the precocial birds. The strong overlap is also detected in the LDA (total error of correct identification: 48.9%) (Figure 3C). In contrast to the perMANOVA only the super-altricial birds can be clearly identified (error of correct identification: 2%). If the four groups are combined into two groups, precocial and altricial birds are clearly separated from each other based on the perMANOVA (F-value: 4.325; p-value: 0.025), while LDA identifies both groups with an error of 22.1% correctly (see Supplementary Tables 6, 7).

PCOA and Multivariate Regression of Linear Measurements
For the second dataset that include only measurements, the first six relevant PCO axes explain about 52.85% of the total variation. Applying taxonomic grouping, PCoA results indicate a certain resolution between the different clades. While the morphospace of Aequorlitornithes occupies most of the positive side of the range of PCO1, the majority of Inopinaves plots on the negative side, overlapping with the morphospace of Colombaves and Palaeognathae. The Strisores are isolated from the other groups plotting at the negative margin of PCO1.
In contrast, all clades overlap with each other along PCO2. While Aequorlitornithes and Inopinaves occupy the total range of PCO2, Strisores, Colombaves and Palaeognathae are located on the negative side ( Figure 5A). Based on the perMANOVA and LDA (total error of correct identification: 3.1%) results, all taxonomic groups are found significantly different ( Figure 5A). For the ASC grouping, the morphospaces the three different ontogenetic stages overlap with each other along PCO1 and PCO2, showing no significant differences in size, although a increasing trend is still recognizable (sum of variances: Stage 1: 0.023; Stage 2: 0.024; Stage 3: 0.028: all p-value > 0.05) ( Figure 5B). Based on perMANOVA, the morphospace position of stage 1 and 3 are significantly different. This is supported by the LDA (total error of correct identification: 26.0%) which correctly identifies stage 3 from other two with a small percentage of error (error of correct identification: 8.0%) ( Figure 5B). However, the other two stages cannot be distinguished from each other. Applying developmental modes (using five groups), altricial birds occupy almost the complete range of PCO1, partly overlapping with all the other groups. While altricial and super-altricial birds, which plot on the negative side of the axis, show a strong overlap, precocial and semi-precocial birds are separated from each other along PCO1. All groups show a strong overlap with each other along PCO2 ( Figure 5C). Based on the perMANOVA, all developmental modes are significantly different from each other except for precocial and semi-altricial birds. This is mainly supported by LDA (total error of correct identification: 6.1%), which found a certain error in the discrimination of altricial and semi-altricial birds. If only two developmental modes are applied, both groups are separate from each other according to perMANOVA (F-value: 22.760; p-value: < 0.001) and LDA (total error of correct identification: 7.6%) (Figure 5C and see Supplementary Tables 9, 10). The Regression Scores, summarizing the variation of the first six PCO axes, correlates with RASC (slope: 0.17; p-value < 0.001), but on a low level (R 2 : 0.31) (Figure 4B). The correlation is driven by PCO2 to PCO6. There is no difference between precocial (slope: 0.17; p-value < 0.001) and altricial (slope: 0.16; p-value < 0.001) birds. When tested against RSP, the Regression Scores shows a weak correlation (slope: 0.34; p-value < 0.001; R 2 : 0.27) (Figure 4D). This correlation is driven by PCO2 to PCO4 and PCO6. In contrast to RASC, the correlation with RSP shows a significant difference between precocial (slope: 0.18; p-value < 0.001) and altricial (slope: 0.37; p-value < 0.001) birds (see Supplementary Table 11).

Major Ontogenetic Changes of Sutural Characters
During the first ontogenetic stage, the majority of skull bones possess physical contacts with minor openings along the suture line (state 1). However, the contacts between maxilla and nasal, premaxilla and mesethmoid, mesethmoid and laterosphenoid, laterosphenoid and parasphenoid, supraoccipital and parietal, and squamosal and exoccipital are still represented by minimal physical sutures with distinct openings (state 0). In contrast, the contacts between maxilla and palatine, jugal and quadratojugal and angular and splenial have a continuous suture (state 2), while the contacts between basioccipital and basisphenoid, basioccipital and exoccipital, basisphenoid and exoccipital, basisphenoid and parasphenoid, and articular and angular are already fused in the first stage (state 3) (Figure 6).
For the second ontogenetic stage, the majority of bone sutures do not change their morphology. Only in 10 out of 35 contacts, distinct changes are traceable in comparison to stage 1. The contact between maxilla and jugal develops into continuous sutures (state 2), while many bones from the postrostral region (i.e., squamosal and exoccipital, squamosal and laterosphenoid, parietal and exoccipital, exoccipital and supraoccipital, laterosphenoid and parasphenoid, and mesethmoid and parasphenoid) and lower jaw (i.e., dentary and splenial, articular and surangular, and articular and prearticular) become fully fused (state 3). Furthermore, we detected seven more contacts, where changes of the sutural morphology are not resolved (Figure 6), indicating an interspecific variation in the fusion process.
From ontogenetic stage 2 to 3 the remaining bone contacts become fully fused (state 3), with a handful of exceptions. The contacts between premaxilla and nasal, and the nasal and mesethmoid develop into continuous sutures (state 2) while the contact between maxilla and jugal (see above) does not change from stage 2 to 3. The sutural morphology between dentary and angular remains unresolved (Figure 6).

DISCUSSION
Studying the ontogenetic process of bone fusion in the skull of 18 extant bird species, we show that the degree of sutural closure correlates with size and allows for the determination of different ontogenetic stages. While the fusion of bones is not synchronous within different skull regions, the process of fusion itself happens usually in a fast manner. In this context, altricial and precocial birds show a significant difference in the rate of fusion. By contrast, on an interspecific level, linear skull measurements are not suitable to determine different ontogenetic stages, but instead contain a strong taxonomic signal for the identification of major bird clades, when analyzed with a multivariate approach.

Bone Fusion
Based on our results, the degree of fusion between different skull bones (expressed as ASC) can be used as reliable proxy for identifying different ontogenetic stages in birds. Thus, our results support the case study of Bailleul et al. (2016) for Dromaius novaehollandiae. Likewise, bone contacts in crocodylian skulls follow an ontogenetic trend, but in an opposite direction, meaning that most bone sutures become wider during growth, while bone fusion is only observed between two frontals and parietals (Bailleul et al., 2016). This is different from the situation in mammals, where ontogenetic fusion of bones is common but varies both on the interspecific and intraspecific level and cannot be used as general proxy for ontogeny-only as specific proxy for particular subclades (Rager et al., 2013). Only the bones from the occipital region seem to share a common pattern among mammals (Rager et al., 2013) even if for instance in Cingulata some adults can still have visible sutures (Le Verger et al., 2020). Thus, although all three amniote groups show bone fusion in the skull, the strong ontogenetic signal detected in this study are most likely linked to high degree of fusion in adult crown birds (Plateau and Foth, 2020), which correlates with a general decreasing trend of sutural disparity. The process of bone fusion seems to be highly constrained in crown birds and part of their ornithuromorph stem affecting also the postcranium (Wang et al., 2017), while the process is more variable in non-ornithuromorph archosaurs, making ontogenetic discrimination more difficult (Irmis, 2007;Hone et al., 2016).
Our results also found some taxonomic traits that allow distinguishing between different bird clades. However, based on the PCoA results, the ontogenetic signal (PCO1) is more than six times stronger than the taxonomic signal (PCO2) and the fusion of bones reflects much better the ontogenetic allometric variation (against RSP) than the interspecific allometric variation (against the geometric mean). Like morphometric measurements (see Köppl et al., 2005), bone fusion can only identify young juveniles prior to skull obliteration, which occurs fairly early in bird life. However, progressive closure of the skull sutures allows a more accurate identification of different ontogenetic stages than skull measurements (see below), highlighting the important nature of connectivity between bones as a morphological property to identify ontogenetic stages in birds (Plateau and Foth, 2020).
As several studies have previously pointed out that the incomplete ossification of the orbital region of birds allowed for the identification non-mature specimens (Pycraft, 1900;Kesteven, 1942;Maxwell, 2008Maxwell, , 2009Green and Gignac, 2020), we re-analyzed the data using a subsample of eight bones from this particular region and compared the results with the complete set of characters. Showing similar results for PCoA, perMANOVA, and LDA than the original dataset (Supplementary Tables 6, 7), the orbital ossification is indeed a reliable proxy for bird ontogeny. As many ornithological collections cannot provide exact information on the age of their skeletal specimens (see above), estimating the ASC of the orbital region allows for a quick determination of different developmental stages in specimens that did not reach skeletal maturity.

Precocial Versus Altricial Ontogeny
The ontogenetic fusion pattern also contains a developmental signal that allows distinguishing between precocial and altricial modes and their subcategories (see e.g., Portmann, 1935;Starck and Ricklefs, 1998b;Botelho and Faunes, 2015), which variation is also capture by PCO2. As the precocial-altricial spectrum follows the phylogeny of birds (Starck and Ricklefs, 1998b;Dyke and Kaiser, 2010;Botelho and Faunes, 2015), this relationship can be expected. Having a greater slope between RASC and RSP, skull bone fusion happens faster in altricial birds than in precocial ones. This faster development could be linked to a higher growth rates occurring in altricial birds (Starck, 1993;Starck and Ricklefs, 1998c). However, a comparison between growth rates and rates of fusion was only able to detect a significant correlation between both parameters when not corrected for phylogeny. Probably related to the small sample size, the correlation disappeared when the phylogenetic relationship was taken into account. In addition to the present putative correlation, altricial birds produce also larger eggs than birds with precocial developmental mode when compared to the female body mass (Dyke and Kaiser, 2010). Due to an isometric relationship between hatchling and egg size, the relative hatchling size of altricial birds (in comparison to adult size) is also bigger than in precocial birds (Deeming and Birchard, 2006). The latter relationship was also detected in our data, as the RSP of the juvenile altricial birds is higher than those of juvenile precocial birds. Thus, these comparisons indicate a potential causal linkage between relative egg size, relative hatchling size, growth rate and rate of fusion with respect to the developmental mode. Although this relationship might be true for most birds, at least the precocial kiwi (Apteryx spp.) clearly represent an exception, as this bird produces the largest egg in comparison to its body mass (Abourachid et al., 2019) but needs multiple years to become fully grown (Bourdon et al., 2009;Heck and Woodward, 2021). However, when comparing only those precocial and altricial birds that fall in the same range of RSP, bone fusion is still slightly faster in altricial birds, but the rate of bone fusion is not significantly different anymore. This is consistent with our finding that most bone fusion happens late in stage 2 and 3 (see below) and the differences in the original slopes relate to the relative hatchling sizes of precocial and altricial birds.

Linear Measurements
In contrast to the fusion pattern, the set of linear measurements do not contain a clear ontogenetic signal when studied in an interspecific context with a multivariate approach. Only adult birds (Stage 3) could be correctly identified with certainty, while juvenile (Stage 1) and subadult (Stage 2) birds cannot be distinguished from each other. In contrast, the linear measurements possess a strong signal for taxonomy and different developmental modes of the precocial-altricial spectrum, which are both captured by PCO1. Thus, similar to the bone fusion dataset, taxonomy and developmental mode are captured by the same axis, which might indicate a strong phylogenetic relation (Starck and Ricklefs, 1998b;Dyke and Kaiser, 2010;Botelho and Faunes, 2015;Prum et al., 2015). Reducing the second data set to only postrostral measurements, the taxonomic and developmental signals are retained but become much weaker (Supplementary Tables 9, 10). Thus, the classification of both groupings seems to rely on the relative length of the beak. This agrees with the results of Marugán-Lobón and Buscalioni (2003), who found that precocial birds possess meso-to longirostral skulls, while the skulls of altricial birds tend to be meso-to brevirostral. When tested against the geometric mean, the linear measurements contain a strong allometric signal that is primarily expressed in PCO1 (Supplementary Table 11). As this coordinate captures variation in the length of the cranium, mandible and beak, the increase in size could be correlated with a phylogenetically driven elongation of the skull, particularly the beak. This pattern of relative elongation of the rostrum is known as cranial evolutionary allometry (CREA) and is present in both mammals (Radinsky, 1985;Cardini, 2019) and birds (Bright et al., 2016;Linde-Medina, 2016;Tokita et al., 2017). In comparison, the linear measurements reflect interspecific allometric variation (against the geometric mean) much better than ontogenetic allometric variation (against RSP). However, on an intraspecific level, the analysis of morphometric measurements is still an appropriate approach to document changes of size and proportions during ontogeny (see e.g., Carrier and Leon, 1990;Köppl et al., 2005;Bennett, 2008;Hanai et al., 2021).

Ontogenetic Changes of Sutural Characters
Our results indicate a generally fast process of fusion, in which bone contacts change directly from open sutures (state 0 or 1) to fully obliterated (state 3). Bone fusion appears first in the posterior part of the lower jaw and the basicranium around the foramen magnum. During the next stage, fusion appears in the anterior part of the mandible and the temporal region, integrating ethmoidal, sphenoidal and occipital bones. At last, fusion appears in the facial region, zygomatic arch, skull roof and around the splenial. Similar patterns have been previously documented for Dromaius novaehollandiae (Palaeognathae), Gallus gallus (Galloanserae), Aptenodytes forsteri, Macronectes giganteus (both Aequorlitornithes) and Myiopsitta monachus (Inopinaves) (Jollie, 1957;Bailleul et al., 2016;Sosa and Hospitaleche, 2018;Piro and Hospitaleche, 2019;Carril et al., 2020), indicating a highly constrained development among birds. Based on this comparison, we see a common fusion trend from posterior to anterior and from ventral to dorsal. Likewise, the skull ossification, taking place during embryogenesis and early postnatal ontogeny, follows a certain sequential pattern. Based on embryological studies (Maxwell, 2008(Maxwell, , 2009Carril and Tambussi, 2017;Arnaout et al., 2021) most bones of the upper and lower jaw, zygomatic arch and the pterygoidpalatine complex ossify prior to bones of the skull roof, occipital and sphenoid region, while the mesethmoid is the last element that becomes bone. Thus, the sequences of prenatal ossification and postnatal bone fusion are in part reversed, when compared to each other. Nevertheless, the majority of fusion events happen ontogenetically late, in stage 2 and 3. Indeed, many altricial birds in our sample still showed a relevant number of sutures while they have already reached 90% of their adult RSP.
Beside this common developmental pattern, we still found some temporal variation in the fusion process between birds during the intermediate second stage. This affected a small number of contacts in the cranium (i.e., premaxillary and mesethmoid, nasal and frontal, parietal, and squamosal) and lower jaw (i.e., dentary, and angular, dentary and surangular, and angular and prearticular), and highlight a certain degree of interspecific variation, which is also reflected in the taxonomic signal of PCO2 in the Dataset 1. On the one hand, this variation is partly caused by the arbitrary division into three ontogenetic stages, which was limited by the sampling of early and intermediate ontogenetic stages for some taxa. Applying a scheme with more stages could potentially lead to a better temporal resolution for the bone contacts in question. On the other hand, this pattern could be also an indication for heterochronic shifts in the fusion pattern within crown birds (see Lecointre et al., 2020). This has to be studied in more detail in the future.

CONCLUSION
In contrast to mammals, the degree of sutural closure between different skull bones allows to identify common ontogenetic patterns across all birds and not only particular groups. The degree of fusion correlates with relative skull size, supporting the case study by Bailleul et al. (2016) on Dromaius novaehollandia. While this study focused primarily on the complete skull, subsampling only bones from the orbital region as suggested by Green and Gignac (2020), lead to similar results. The sutural closure happens relatively fast during ontogeny and shows a common directional trend from posterior to anterior and ventral to dorsal. Despite these patterns, the developmental mode has a strong influence on the speed of bone fusion, which is faster in altricial birds, mirroring patterns of relative egg sizes, relative hatchling sizes and the growth rates within the precocial-altricial spectrum. In contrast, linear skull measurements do not contain a clear ontogenetic signal but allow for a better identification of taxonomical groups, when analyzed on an interspecific level with a multivariate approach. This comparison indicates that at least in birds the connectivity of bones is an important morphological property for the identification of ontogenetic stages beyond the intraspecific level.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because the specimens investigated for this study were all housed in osteological collections of various natural history museums in Europe. Consequently, the specimens were already dead when we studied them.