Intra-specific variation and allometry of the skull of Late Cretaceous side-necked turtle Bauruemys elegans (Pleurodira, Podocnemididae) and how to deal with morphometric data in fossil vertebrates

Background Previous quantitative studies on Bauruemys elegans (Suárez, 1969) shell variation, as well as the taphonomic interpretation of its type locality, have suggested that all specimens collected in this locality may have belonged to the same population. We rely on this hypothesis in a morphometric study of the skull. Also, we tentatively assessed the eating preference habits differentiation that might be explained as due to ontogenetic changes. Methods We carried out an ANOVA testing 29 linear measurements from 21 skulls of B. elegans taken by using a caliper and through images, using the ImageJ software. First, a Principal Components Analysis (PCA) was performed with 27 measurements (excluding total length and width characters; =raw data) in order to visualize the scatter plots based on the form variance only. Then, a second PCA was carried out using ratios of length and width of each original measurement to assess shape variation among individuals. Finally, original measurements were log-transformed to describe allometries over ontogeny. Results No statistical differences were found between caliper and ImageJ measurements. The first three PCs of the PCA with raw data comprised 70.2% of the variance. PC1 was related to size variation and all others related to shape variation. Two specimens plotted outside the 95% ellipse in PC1∼PC2 axes. The first three PCs of the PCA with ratios comprised 64% of the variance. When considering PC1∼PC2, all specimens plotted inside the 95% ellipse. In allometric analysis, five measurements were positively allometric, 19 were negatively allometric and three represented enantiometric allometry. Many bones of the posterior and the lateral emarginations lengthen due to increasing size, while jugal and the quadratojugal decrease in width. Discussion ImageJ is useful in replacing caliper since there was no statistical differences. Yet iterative imputation is more appropriate to deal with missing data in PCA. Some specimens show small differences in form and shape. Form differences were interpreted as occuring due to ontogeny, whereas shape differences are related to feeding changes during growth. Moreover, all outlier specimens are crushed and/or distorted, thus the form/shape differences may be partially due to taphonomy. The allometric lengthening of the parietal, quadrate, squamosal, maxilla, associated with the narrowing of jugal and quadratojugal may be related to changes in feeding habit between different stages of development. This change in shape might represent a progressive skull stretching and enlargement of posterior and lateral emargination during ontogeny, and consequently, the increment of the feeding-apparatus musculature. Smaller individuals may have fed on softer diet, whereas larger ones probably have had a harder diet, as seen in some living species of Podocnemis. We conclude that the skull variation might be related to differences in feeding habits over ontogeny in B. elegans.


INTRODUCTION Principal Component Analysis and fossil sampling bias
Paleontological data are intrinsically scarce (Strauss, Atanassov & Oliveira, 2003;Hammer & Harper, 2006), leading to incomplete data sampling. This limitation impacts several approaches in paleontological studies, especially inter-specific variation analyses. Although there are some methods proposed to deal with missing entries in fossil quantitative datasets (e.g., Norell & Wheeler, 2003;Strauss, Atanassov & Oliveira, 2003), sometimes the study relies on an exploratory evaluation of general structure and Principal Component Analysis (PCA) is commonly used for this purpose.
That being said, we consider that all turtle specimens found at the Tartaruguito site might correspond to subadults to adult ages, and it is reasonable to assume that all B. elegans individuals collected in the Tartaruguito site might have belonged to a single population (agreeing with Henriques et al., 2002;Henriques et al., 2005;Henriques, 2006;Romano & Azevedo, 2007). Indeed, as suggested by Romano & Azevedo (2007), this single population would consist on different generations of turtles' corpses grouped in the same locality. One might consider that size differences might be due to sexual dimorphism (R Hirayama & S Thomson, pers. comm., 2015), in which the females would be larger and have more posteriorly extended carapaces than the males. However, sexual dimorphism on podocnemidid turtles can be assessed only on shell shape and our data is based mostly on isolated skulls (see Material and Methods). As a consequence, although it is possible that sexual dimorphism may affect measurements captured in this study, we did not consider it, given the lack of evidence to assume such outcome. Also, to our knowledge, skull shape differences related to sexual dimorphism has never been described to podocnemidid turtles yet. Moreover, Romano & Azevedo (2007) were not able to reject the single population hypothesis using shell measurements (from both plastron and carapace) in a morphometric approach neither to describe sexual dimorphism in the data, concluding that the differences were due to ontogenetic variation. Therefore, we highlight that we are assuming the population definition of Futuyma (1993), as taken on by Romano & Azevedo (2007), that a population is a conjunct of semaforonts temporally connected, i.e., a sequence of individuals from different generations, and limited in a restricted space; in this case, the Tartaruguito site. By assuming this, we explicitly follow Hennig's (1966) semaphoront concept, on which a species is modifiable (i.e., not strictly typological) and represented by a sequence of generations.

Objectives
Efforts to study fossil materials may be hampered by difficulty in accessing foreign collections. It can narrow and even preclude their studies. In addition, given the missing data problem inherent to fossil record, the way one treats the missing entries in morphometric studies can affect the results and conclusions. Regarding the use of caliper or ImageJ in taking measurements, here we tested both approaches by taking linear measurements for morphometric studies based on photographs (e.g., Bailey, 2004) and also evaluated how different approaches designed to deal with missing data can impact results of exploratory statistical procedures and data interpretation by comparing two different substitution algorithms of missing entries. These procedures are exemplified using a real paleontological data set and with paleobiological inferences. Considering the case-study, we explored the variation in skulls among individuals of Bauruemys elegans from different ages and generations. Also, we describe the differences in skull morphology along the ontogeny of the species and discuss the probable consequences of such variation to the diet preferences changes along the growth.

Sample and characters
Twenty-one skulls of Bauruemys elegans were examined in this study, including the type series plus nineteen topotypes: AMNH-7888, LPRP0200, LPRP0369, LPRP0370,  (Fig. 2) that decompose the overall shape of the skull in order to take measurements between two landmarks. Since most of the specimens have deformation and breakage, we could not perform a geometric morphometric analysis using the landmarks because the taphonomical bias would incorporate error to the analysis of form and shape. Thus, we used the landmarks to set up 29 traditional morphometric characters that correspond to a linear measurement between two landmarks (all characters are described in Table 1). Also, the use of landmarks to set up the measurements is useful to maintain the same anatomic references for all characters in each specimen, since the landmarks enable a better description of morphological variation and establishment of linear measurements, as performed by Romano & Azevedo (2007) with shell morphometric characters. All measurements were taken on the same side of the skull (right side) unless the characters could not be measured due to deformation or breakage. We are aware that deeper structures (z-axis) can influence the straight line between two landmarks in 2D images and used ImageJ version 1.47 (Rasband, 1997) to take the measurements after comparing its accuracy with the caliper (Mariani & Romano, 2014). This procedure was necessary because we obtained photos of skulls in dorsal, ventral, and lateral views housed in foreign collections and did not collect the described measurements (see Table 1) using caliper in such specimens because they were analyzed prior to this study. The error test between measurements taken using caliper and ImageJ using part of the sample are  Table 2 for vectors description). Abbreviations: bo, basioccipital;bs, basisphenoid;ex, exoccipital;fpp, foramen palatinum posterius;fr, frontal;ju, jugal;mx, maxilla;op, opisthotic;pa, parietal;pal, palatine;pf, prefrontal;pm, premaxilla;po, postorbital;pt, pterygoid;ptp, processus trochlearis pterygoidei;qj, quadratojugal;qu, quadrate;sq, squamosal;so, supraoccipital;vo, vomer. Skull lineation from Gaffney et al. (2011, p.72). described below. We followed the bone nomenclature of Parsons & Williams (1961) and extended by Gaffney (1972) and Gaffney (1979) (see all abbreviations after Conclusion topic).

Preliminary analysis: caliper vs. ImageJ
Before carrying out others statistical analyses, we compared the same characters data set (Data S1) of a sub-sample by using two different approaches (=treatments): measurements taken using caliper and measurements taken using photographs via ImageJ. This comparison was necessary in order to evaluate whether or not the two measurements methods are significantly different. Then, we performed an One-way Analysis of Variance (ANOVA) comparing the 29 measurements in 12 specimens (LPRP0200, LPRP0369, LPRP0370, MN4322-V MN4324-V, MN6750-V, MN6783-V, MN6786-V, MN6787-V, MN6808-V, MN7017-V, and MN7071-V). Two groups of variables were established: measurements taken directly from specimens using caliper (preliminary data set 1) and the same characters taken from photographs of the same specimens using ImageJ (preliminary data set 2). All characters taken using photographs/ImageJ that did not show significant differences to their correspondents taken by caliper were used on the subsequent statistical analyses of form and shape differences among the sample of Bauruemys elegans. By doing that, the sample was increased without including error and incomparable characters (i.e., by using different measurement techniques).
We found most of the measurements do not differ statistically (p > 0.05) between the two treatments (caliper and ImageJ; Table 1). However, one measurement, length of maxilla (LMX), had statistical difference (p = 0.017) between the treatments, because the maxilla is a curved structure and thus the landmarks are in different positions (LM 24 is deeper and farther from the camera in relation to LM11) in relation to the plane the picture was taken. Given that no statistical differences were found in almost all characters, ImageJ could be an economic and time-saving tool for morphometric analyses from photographs (2D), and could be applied by scientists at distant institutions.
The study in situ of the material is preferable, although pictures are an economic alternative in where cases one is not able to handle the material. We must aware that one have to choose one of the two treatments to construct a morphometric matrix, otherwise it will be composed of values obtained by two diffent methods.

Univariate, multivariate and allometric analyses
Three analyses using the complete sample were carried out: (1) a descriptive statistics (mean, standard deviation, median, variance, maximum and minimum values) of all characters (Data S2), (2) an allometric analysis of length and width characters correlating them to total length and width measurements (Data S3), and (3) a multivariate non-parametric exploratory statistics via Principal Component Analysis (PCA). The latter was divided into two different PCAs: (3.1) using 27 characters from the raw data matrix (total length and width characters were excluded in this analysis; Data S4-because PCA is sensitive to large variations in the original measurements), and (3.2) using 27 characters that correspond to the proportions of each character from the raw data (i.e., original measurements) represented by its length or width characters divided by each individual total length or width (e.g., the length of MCZ4312 postorbital divided by the total length of the skull of this specimen; see complete data in Data S5). All statistical analyses were performed using the software PAST version 3.05 (Hammer, Harper & Ryan, 2001).
In the allometric analysis (analysis 2, Data S3), all characters were previously logtransformed and a linear regression was carried out separately for length and width characters, using the least square fitting approach for residuals. We established the allometries by considering the regression's slope, i.e., the coefficient a, as following: positive allometry (a > 1), negative allometry (1 > a > 0), enantiometry (a < 0), and isometry (a = 1).
In the first PCA approach (3.1) we excluded total length and width characters because of their high influence on the PCA result, since higher values compose the majority of the summarized variance in PC's (Mingoti, 2013), and because of the redundancy between these measurements and the others. We also assessed differences by applying two different substitution algorithms for missing data in PAST, using the default ''mean value imputation'' option (i.e., missing data are replaced by the column average), and the alternative ''iterative imputation'' option, which computes a regression upon an initial PCA until it converges to missing data estimations, replacing missing data by such estimations (Ilin & Raiko, 2010). The latter is recommended and, after comparing both results, we selected it (see Data S3 to visualize PCA results computed using PAST's default option approach). The second PCA (3.2) was conducted to remove the effect of size (=growth) and perform an exploratory analysis of the shape alone. Six specimens were removed from this analysis because they were broken and the total length or width measures were not measurable.
The univariate analysis was made in order to quantify and describe the variation of the characters set in Bauruemys elegans skull, using the assumption that the sample is representative of a single population. The linear regression analyses allowed us to make inferences about osteological shape change related to size change, i.e., related to growth, by assuming that bigger specimens are older than smaller ones. This approach is, therefore, a study of allometry (sensu Huxley & Teissier, 1936;Huxley, 1950;Gould, 1966;Gould, 1979;Somers, 1989;Futuyma, 1993) and the assumption of correlation between size and aging is based on continuous growth to be common on extant turtles (Klinger & Musick, 1995;Shine & Iverson, 1995;Congdon et al., 2003). Since the use of a parametric statistic was infeasible due to the nature of the sample (i.e., a small dataset that do not show homoscedasticity and normality), the PCAs were used to search for a structure of the data that matches to the pattern found by Romano & Azevedo (2007) using shell characters (i.e., all individuals plotted inside the 95% confidence ellipse). If the pattern observed is similar to previous morphometric and taphonomic inferences, then the variation is not enough to assume that the sample represents different populations of Bauruemys elegans or a different species (see 'Taxonomic considerations on the sample'). In other words, since a parametric test is not feasible with statistical confidence, the lack of structure in the PC plots were herein interpreted as a failure to falsify the single population hypothesis. All principal components were, therefore, analyzed but we present only those with higher variance.

Descriptive analysis
The results of the descriptive statistics are summarized in Table 2. As expected, values of total length and width (TLS and WLS) were the most variable among all measurements, because the variation scale in these characters is greater than in others measurements. Characters of the bones forming the upper temporal fossa (i.e., PA, QJ, SQ, QU and OP) had great variation, with the parietal being the most variable in length (SD = 6.45) and the least variable in width (SD = 2.94), whereas quadratojugal obtained the smallest variation in length (SD = 2.38) and the greatest in width (SD = 4.03). Among the characters of the bones forming the lower temporal fossa (i.e., JU, MX, PO, PT and PAL), the variation Table 2 Descriptive statistics of all data. Descriptive statistics of the three sorts of characters analyzed (total length and width, comprised measurements, and proportions of the measurements), including mean values (Mean), median values (Median), standard deviation values (SD), number of entries (N), and maximum and minimum values (Max-Min). All measurements are expressed in millimeters, except unscaled proportions between two measurements.  in length was in general greater than in width. Postorbital and maxilla had almost the same variation in length (SD = 4.12 and SD = 4.11, respectively); WPO had the smallest variation within the group of bones forming the lower temporal fossa (SD = 1.83); and the stretch of the maxilla had the greatest variation (SD = 7.63) of all characters measured. Characters of the other bones had smaller values than the aforementioned bones, with the exception of WPO which was smaller than LFR (SD = 2.08), LVO (SD = 1.95), LBO (SD = 2.12),WFR (SD = 1.88) and WBS (SD = 2.19).

Raw data
Replacing missing data with mean values. By using the ''mean value imputation'' approach, a total of 70.32% of the variance was comprised by the first three principal components (PC1 = 42.15%; PC2 = 16.82%; PC3 = 11.35%), so that the others were less significant for the analysis by following the broken stick model, and are not presented. We interpreted that PC1 variation is due to size variation because an approach using all characters has shown a similar plot (see Fig. 6A). PC2 and PC3 seems to represent shape differences between individuals. In all PC individual projections (Figs. 6A and 6B) most of specimens were included inside the 95% confidence ellipse. Two exceptions are MCZ4123 and MN7071-V, which have not been included in the ellipse when PC1 vs. PC2 were considered (Fig. 6A); also the former was outside the ellipse in PC2 vs. PC3 scatter plot (Fig. 6B), indicating form differences of these specimens. However, both specimens have suffered different degrees of crushing due to taphonomic bias and that is likely the reason for this result.  In PC1' loadings (Table 3) Replacing missing data with regression estimation. The alternative missing data approach (i.e., ''iterative imputation''; Fig. 6C) generated two principal components which comprised 88.96% of the total variance (PC1 = 53.01%; PC2 = 35.95%). In contrast with the previous approach, PC1 was interpreted as representing shape and PC2 reflected size variations. In addition, all specimens were included inside the 95% ellipse in PC1vs.PC2 scatter plot. The specimen MN7017-V, interestingly, was excluded from the ellipse when considering PC2 vs. PC3, but the percentage of variance represented by PC3 is too low (PC3 = 3.28%) to assume any difference from the others individuals. We agree with Ilin & Raiko (2010) and prefer to choose the iterative imputation approach for dealing with missing entries (see discussion on 'The single population hypothesis'). Then, discussions concerning the form variation in our data are related to PCA analysis using iterative imputation.

Shape characters (proportions)
Replacing missing data with mean values. When applying ''mean value imputation'', 53.99% of the variance were comprised by the first two principal components (PC1 = 35.29%; PC2 = 18.70%), both corresponding to shape, as all units of measurements were removed through the division of characters before carrying out the analysis. All specimens were comprised into the 95% confidence ellipse (Fig. 7A).
The  Table 4 for all loading values). It is interesting to note that most of highly-related proportions were in reference to bones associated either with feeding apparatus (squamosal, parietal, quadratojugal and jugal) or catching food and trituration surface (maxilla).
Replacing missing data with regression estimation. The ''iterative imputation'' substitution model of missing data explained 77.35% of the variance comprised by two principal components (PC1 = 45.49%; PC2 = 31.86), both representing shape. All specimens were included in the confidence ellipse (Fig. 7B), thus shape differences do not indicate possible different populations or species.  Table 4). These loadings represent shape changes in regions of the skull that are associated with muscles' attachments as well as trituration surfaces (see below).

The single population hypothesis
In this section, we discuss the single population hypothesis considering two fronts, one underlied on the taphonomy of the Tartaruguito locality, and another on the possibility of  the skull variation represent one or more specimens of species Roxochelys wanderleyi in the sample, a shell-only species also found at the site.

The depositional context at the "Tartaruguito" site
The depositional environment at the Pirapozinho site is well-known from previous studies, which point out to seasonal floods in which turtles might have gathered in water bodies for foraging, followed by droughts that caused their death (Soares et al., 1980;Fulfaro & Perinotto, 1996;Fernandes & Coimbra, 2000;Henriques et al., 2002;Henriques et al., 2005;Bertini et al., 2006;Henriques, 2006). This is, consequently, a case of several seasonal non-selective death events, with individuals representing semaphoronts connected temporally (between generations), thus comprising a single population (agreeing with Futuyma's, 1993 population definition and used by Romano & Azevedo, 2007). We failed to disprove the null hypothesis that all individuals belong to a same population of Bauruemys elegans, agreeing with Romano & Azevedo (2007) conclusion using post-cranium data.

Taxonomic considerations on the sample
Many skulls sampled show taphonomic effects, such as cracks and crushing (Fig. 8). For instance, MN7071-V is notably the largest specimen of the sample and is represented in the uppermost positive side of the size-related PC2 axis (Fig. 6C). Although it is indeed a big specimen, it was clearly a taphonomic effect (crushing) that caused it to be larger than it really was (Fig. 8A). On the other hand, Bertini et al. (2006) indicated that turtle bodies have suffered little transportation or crushing in Tartaruguito site. We agree with this taphonomical interpretation of the site, as most specimens do not show huge breaks (Figs. 8B and 8D) that could cause misinterpretation of the morphometric results (the case of MN7071-V is an exception in our sample in this respect). Another aspect is the presence of polymorphism in B. elegans. Romano (2008) presented an unusual carapace for the specimen MN7017-V, as having a seventh neural bone, differing from the diagnostic number of six neurals for this species, and with the diagnostic four-squared second neural bone not contacting first costals (Suárez, 1969b;Suárez, 1969c;Kischlat, 1994;Gaffney et al., 2011). The morphometric analysis performed by Romano (2008), using only shells, did not reveal significant statistical differences between MN 7017-V and other B. elegans specimens. We have included the MN7017-V skull in our analysis, and there was no variation to state anything apart from Romano's (2008) conclusion that it is probably a polymorphic B. elegans specimen (Fig. 6C). Still, we reevaluated this skull and found the diagnostic characters for B. elegans. Therefore, all skulls included in our study belong to the same species (i.e., B. elegans).
Among the five valid fossil turtle species found throughout the Bauru Basin, only two have been collected at the Pirapozinho site so far (Romano et al., 2013). The first is B. elegans, which is recognized by both skull and shell materials; the second is Roxochelys wanderleyi Price, 1953 based only on shell material (Price, 1953;De Broin, 1991;Oliveira & Romano, 2007;Romano & Azevedo, 2007;Gaffney et al., 2011;Romano et al., 2013). So far, none R. wanderleyi with skull-shell associated body parts were collected, and thus we cannot claim that the skulls found at Tartaruguito site belong to this species until a skull-shell R. wanderleyi specimen be found, since all skulls analyzed here can be safely identified as belonging to B. elegans.

Ontogenetic changes in B. elegans skull
Once we have assessed that all specimens belong to the same species and are likely from the same population, we are able to discuss the skull variation in the sample assuming as due to inter-populational variety. For the sake of organization, we divided the discussion into two parts, based on the anatomical regions of the turtle skull: upper temporal fossa and lower temporal fossa, following Schumacher (1973), Gaffney (1979) and Gaffney, Tong & Meylan (2006). We have chosen this organization because the bones we found most associated with the principal components in the two PCA analyses constitute these two regions and are generally involved in aspects of the feeding mechanisms of turtles, either as muscles attachments or forming triturating surfaces.

Bones of the upper temporal fossa and skull roofing
The temporal emargination of podocnemidid turtles is formed by the dorsal, horizontal plate of the parietal, the quadratojugal and the squamosal, with no contribution of the postorbital (Gaffney, 1979;Gaffney et al., 2011). This region (and bones) is associated with the origin of the adductor muscle fibers (m. adductor complex; Figs. 9A and 9B) (Schumacher, 1973;Werneburg, 2011;Werneburg, 2012;Jones et al., 2012;Werneburg, 2013), which run through cartilago transiliens of the processus trochlearis pterygoidei of the pterygoid and insert at the coronoid process of the lower jaw (Schumacher, 1973;Gaffney, 1975;Gaffney, 1979;Lemell, Beisser & Weisgram, 2000;Werneburg, 2011). These muscles promote the closure of the mouth, thus it is reasonable to associate the attachment surface to bite force and the latter to the prey hardness. Yet on the ventral flange of the squamosal originates the m. depressor mandibulae (Schumacher, 1973;Gaffney, Tong & Meylan, 2006;Werneburg, 2011;Fig. 9B), which causes the abduction (=opening) of the mandible. The variation in this area of the skull in turtles was a matter of some studies (e.g., Dalrymple, 1977;Claude et al., 2004;Pfaller, Gignac & Erickson, 2011), which indicated allometric ontogenetic growth patterns of the bones in these regions. These authors were able to identify a high correlation with the increasing of muscle mass and shift in feeding features (Dalrymple, 1977;Pfaller et al., 2010;Pfaller, Gignac & Erickson, 2011). Moreover, there are changes in skull shape associated to the aquatic environment and foraging strategies, as suggested for emydid and testudinoid turtles by Claude et al. (2004). Although these studies focused on hide-necked turtles, the same morphoecological patterns can be applied to side-necked turtles, since there are habitat occupation similarities between side-necked and hide-necked turtles with implications to the skull morphology due to morphofunctional constraints (Schumacher, 1973;Lemell, Beisser & Weisgram, 2000), besides the adaptive selection regarding fresh water feeding strategies (see Lauder & Prendergast, 1992, Aerts, Van Damme & Herrel, 2001and Van Damme & Aerts, 1997 for feeding strategies in freshwater turtles).
The high variance and positive allometric growth of the parietal (LPA: a = 0.38; WPA: a = 0.32), quadratojugal (LQJ: a = 0.16; WQJ: a = −0.06) and squamosal (LSQ: a = 0.30) lead to an increase in temporal emargination and, consequently, a greater area for attachment of the m. adductor mandibulae externus. The consequence of this would be the generation of large forces and high velocities during the fast closing phase of an aquatic feeder, as seen in Pelusios castaneus (Lemell, Beisser & Weisgram, 2000), and even a more powerful bite for crushing harder prey, as seen in Sternotherus minor (Pfaller, Gignac & Erickson, 2011). In addition, the lenghten of the squamosal would allow a greater insertion area of the m. depressor mandibulae and muscles of the hyobranchial apparatus (e.g., m. constrictor colli) (Schumacher, 1973;Gaffney, 1979;Claude et al., 2004;Gaffney et al., 2011;Werneburg, 2011). The m. depressor mandibulae is useful for an increased gape opening speed and the hyobranchial apparatus musculature is involved in backwards water flow generation by the lowering of the hyoid apparatus, two characteristics well reported for other pleurodire turtles as Chelodina longicollis, Chelus fimbriatus and Pelusios castaneus (e.g., Van Damme & Aerts, 1997;Aerts, Van Damme & Herrel, 2001;Lemell, Beisser & Weisgram, 2000;Lemell et al., 2002). Moreover, Claude et al. (2004) demonstrated that aquatic turtles with suction feeding mode possess longer skulls than terrestrial turtles, the squamosal being most prominent bone involved in this elongation and functionally related to the style of prey capture (=suction) as a support for mandible and hyoid muscles.
Also, Gaffney et al. (2011), in a comparison with other podocnemidid turtles, indicated B. elegans as having a ''skull relatively wide and flat'' (p. 12), which could be observed by the increasing of some bones, specially the postorbital (Figs. 3G and 4H), parietal (Figs. 3A and 3J), quadratojugal (Figs. 3I and 4F) and jugal (Figs. 3C and 5B). Comparing the postorbital allometry (better discussed below) with those of the bones in contact with it in the skull roof (frontal, parietal, jugal and quadratojugal; Gaffney et al., 2011), we observe an influence of the positive growth of the former into the others, leading to flattening and widening of the skull.
In a study assessing the bite performance in turtles, Herrel, O'Reilly & Richmond (2002) suggested that a higher skull is efficient in promoting stronger bite forces, specially in species which feed on hard prey, but they also pointed out that additions in bite forces may be achieved by ''getting longer and larger'' skull with no increasing in skull height. Thus, in addition to provide gains in muscle attachment area, by the growing of parietal, quadratojugal and squamosal, leading to a longer skull, a stronger bite and possibly a change in diet along the ontogeny. Also, the allometric growths of most of skull bones, particularly the positive allometry of the postorbital, suggests a more roofed skull in larger adults of B. elegans bigger adults. Given the allometric patterns aforementioned, B. elegans may have had a wide and flat but long skull, which would have compensated the loss of muscle volume and attachment area caused by widening and flattening the skull (Herrel, O'Reilly & Richmond, 2002). Correlations between a more emarginated skull and increases in the volume of the adductor muscles were also explored in a cranial evolutionary framework of stem-turtles by Sterli & De la Fuente (2010).
At last, Gaffney, Tong & Meylan (2006) and Gaffney et al. (2011) scored a character based upon the contact between quadratojugal and parietal bones (char. 13 of Gaffney, Tong & Meylan, 2006;char. 5 of Gaffney et al., 2011). They also state that this contact is present in Hamadachelys + Podocnemididae clade, with a large quadratojugal (state 1), in contrast to most of other Pelomedusoides (state 0: contact absent, as seen in Pelomedusidae, Araripemydidae and many bothremydids (e.g., Kurmademydini, Cearachelyini and Bothremydini); state 2: contact present with small quadratojugal in some Taphrosphyini, Bothremydidae). Indeed B. elegans possess a large quadratojugal, which means that the reduction of the postorbital evolved after Bauruemys node of divergence, as confirmed in performed cladistic analyses. However, we found a greater increasing (positive allometry) of the two measurements of the postorbital and this might have influenced the growth of parietal and quadratojugal, as well as the jugal (see below), so that the state 1 seen in B. elegans is possibly a consequence of allometric changes. This is easily seen when comparing the enatiometry of the width of the quadratojugal (WQJ: a = −0.06) and the slight increasing in the length of this bone (LQJ: a = 0.16) with the postorbital measurements. It also could have influenced the growth of the parietal, but to a lesser extent, as seen in the allometries of this bone (LPA: a = 0.38; WPA: a = 0.32).
When comparing the stem-Podocnemidinura species (i.e., Brasilemys, Hamadachelys) and stem-Podocnemididae (e.g., Bauruemys, Peiropemys, Pricemys and Lapparentemys), with the crown-Podocnemididae (i.e., Podocnemis lineage + Erymnochelys lineage) (Gaffney et al., 2011;Fig. 10), it is clear that an increasing in the parietal-quadratojugal contact has occurred along the podocnemidid lineage, and consequently led to a more roofed and less emarginated skull. We suggest that in B. elegans the small contact is due to the positive growth of the postorbital resulting in a more emarginated skull than other podocnemidids, Within Pan-Podocnemididae, the contact PA-QJ is present and the skull roofing increased from a less roofed condition, found in Brasilemys josai and Hamadachelys, to a continuous increasingly growing well roofed condition within Podocnemididae, exemplified by Bauruemys elegans, Lapparentemys vilavillensis and Podocnemis unifilis, up to a fully roofed morphology in Peltocephalus. Cearachelys placidoi and T. congolensis modified from Gaffney, Tong & Meylan (2006); Brasilemys josai redrawn from Lapparent de Broin (2000); all others skulls modified from Gaffney et al. (2011). as described by Gaffney et al. (2011). Yet within crown-Podocnemididae this bone suffered the opposite effect (i.e., small growth), showing variations in size and even being absent in some species (e.g., Podocnemis sextuberculata, Ruckes, 1937;Gaffney, 1979;Gaffney et al., 2011), though the emargination is still great. On the other hand, in the Erymnochelys lineage the postorbitals are large but the quadratojugal and parietal are large as well, leading to a greater contact between these bones and a well-roofed but less emarginated skull, being a reversion in Bairdemys venezuelensis and B. sanchezi within the Erymnochelys lineage (Gaffney et al., 2011). Therefore, the increase or decrease in the temporal emargination within Podocnemididae could be due to variation of allometric patterns in bones that form the skull roof, particularly the postorbital, quadratojugal, and parietal, among different lineages. Given that observation, we speculate that characters related to the form of the aforementioned bones (postorbital, squamosal, and parietal) are potencially more susceptible to homoplasy.

Bones of the lower temporal fossa
The lower adductor chamber in Pelomedusoides is formed externally and laterally by the jugal and quadratojugal, with the addition of the maxilla in some cases (e.g., Podocnemis spp. and Bairdemys sanchezi). The well developed cheek emargination, found in most but not all podocnemidid turtles (the exceptions are all within the Erymnochelys lineage but Bairdemys spp., Cordichelys antiqua and Latentemys plowdeni), is also part of the adductor chamber (Gaffney, 1979;Gaffney, Tong & Meylan, 2006;Gaffney et al., 2011). Internally and medially, the postorbital, the jugal, and the pterygoid compose the septum orbitotemporale, partially separating the fossa orbitalis from the fossa temporalis; along with the palatine, they aid to support the processus trochlearis pterygoidei of the pterygoid (Gaffney, 1975;Gaffney, 1979;Gaffney, Tong & Meylan, 2006). There is a passage medially to the process of the pterygoid and the septum orbitotemporale, running from the fossa orbitalis to the fossa temporalis, the sulcus palatinopterygoideus. The palatine and pterygoid form the floor of its passage, whereas the parietal, postorbital and frontal limit its upper portion. In this region, the m. adductor mandibulae fibers run through the processus trochlaris pterygoidei, and the m. adductor mandibulae internus (i.e., m. pterygoideus and pars pseudotemporalis; Fig. 9B) mostly originates throughout the pterygoid and parietal bones (Schumacher, 1973;Lemell, Beisser & Weisgram, 2000;Lemell et al., 2002;Werneburg, 2011). The m. adductor mandibulae internus fibers are involved in the jaw-closure system by generating counter forces (protraction) to the m. adductor mandibulae externus (retraction) (Schumacher, 1973;Lemell, Beisser & Weisgram, 2000;Lemell et al., 2002;Figs. 9C and 9D).
Variation of the upper temporal fossa has been studied in different turtles, such as various trionychids (Dalrymple, 1977) and Chelydra serpentina (Herrel, O'Reilly & Richmond, 2002). However, few studies report the variation of the lower adductor chamber, although both the upper and lower temporal fossa are anatomically and functionally coupled (Schumacher, 1973). Dalrymple (1977) identified a positive allometry in the width of the ''temporal passageway'' in trionychids. This area is related to the cryptodire pulley system (i.e., a processus trochlearis formed by the quadrate and opisthotic) and is analogous to the pleurodire pterygoid process, and thus can be comparable functionally (Gaffney, 1979). Herrel, O'Reilly & Richmond (2002) concluded that the increase of the bite force in turtles is due to either the increased height of the skull, leading to a more open angle of the processus trochlearis in relation to skull longitudinal axis, or to enlargement (in width and length) of the skull, because it allows more area for muscle attachment and volume. We observed the same pattern of growth change in B. elegans, as evidenced by the positive allometry of the parietal, postorbital, palatine and pterygoid bones. Other features were observed by Dalrymple (1977) in trionychids (e.g., height and width of the supraoccipital crest, lengthening of the squamosal crest and a development of a horizontal crest in the parietal) and were correlated to changes in skull shape with a shift in feeding habits, from softer to harder preys as individuals age. Again, this seems to be the case in B. elegans, as evidenced by the positive allometry of the squamosal and parietal bones.
The bones that mainly compose the skull rostrolaterally and the lateral emargination revealed a correlated allometric shape shift. Even so, jugal and maxilla showed small allometric variation (Figs. 4B, 4C, 6A, and 6B). The reduction of the jugal (WJU: a = −0.23) and quadratojugal (WQJ: a = −0.06) along with the small growth of the maxilla (WMX: a = 0.19) demonstrate a decrease in height at the anterior portion of the skull. Because of the contact between jugal and quadratojugal with the postorbital (and its increase; see previous topic), we suggest that the latter would possibly has affected the growth of the former bones. Moreover, the strong development of the postorbital would ultimately affect the width of the maxilla, which in turn would also affect the jugal. In contrast, the lengthen of this bone would be less affected (LMX: a = 0.39). In addition, there is a considerable increment in the stretch of maxilla (SMX: a = 0.70) (Fig. 3H) leading to a broader rostrum. Yet this could allow a greater area for crushing, as observed by Kischlat (1994) for B. elegans, but also related to ontogenetic growth (T Mariani, pers. obs., 2016). All these allometric changes indicate that B. elegans owns a more flattened and wider skull (Gaffney et al., 2011), which could have allowed greater bite forces generation (Herrel, O'Reilly & Richmond, 2002).
There are other morphological implications in which the lower adductor chamber bones are involved and that are worth discussing. As previously pointed out, three bones compose the septum orbitotemporale: pterygoid, jugal, and postorbital (Gaffney, 1979;Gaffney, Tong & Meylan, 2006). Together with the palatine, these three bones provide support for the processus trochlearis pterygoidei, where upon runs the tendon that connect the m. adductor externus complex into the lower jaw (Schumacher, 1973;Gaffney, 1975;Gaffney, 1979;Lemell, Beisser & Weisgram, 2000;Gaffney, Tong & Meylan, 2006;Werneburg, 2011). Nearby the process, many muscle fibers originate or cross towards their insertions points (Schumacher, 1973;Werneburg, 2011). The temporal emargination at the upper adductor chamber becomes more emarginted during growth. As a consequence, the attachment area for m. adductor mandibulae externus increase during aging, potentially generating stronger bite forces. The consequence of this temporal emargination indentation is that the trochlear process would must become more robust to support higher forces. We interpret that the positive allometries of pterygoid (LPT a = 1.37), postorbital (LPO a = 1.25 and WPO a = 1.36), and palatine (LPAL a = 1.11) could be a response to this robustness of the trochlerar process during growth. In other words, they would act together by giving more resistance to the area in which the high forces created by the m. adductor mandibulae externus are applied. Gaffney (1979) suggested this robustness occurs because muscle volume increase and, consequently, higher bite forces, so these three bones would reinforce the septum orbitotemporale in order to support and do not break when muscles are contracted. In addition to such reinforcement, the growth of palatine could be associated with a larger area for crushing preys such as mollusks and crustaceans, as pointed out by Kischlat (1994).
The m. adductor mandibulae internus and m. adductor mandibulae posterior (Fig. 9B), which originate at the quadrate, prootic, pterygoid, palatine, postorbital and the descending process of the parietal (Schumacher, 1973;Werneburg, 2011), are important during the jawclosure phase. The importance of these muscles has been debated for early tetrapods with flat skull and aquatic lifestyle (e.g., Temnospondyli and Lepospondyli;Frazzetta, 1968), in which the internal muscle might have assumed the main function of closing the jaw (Werneburg, 2012). This also occurs in turtles with flat skulls and with poorly developed crista supraoccipitalis (e.g., Chelidae; Werneburg, 2011;Werneburg, 2012). However, B. elegans does not have a skull as flat as chelids, but has a long supraoccipital bone as well as a greater temporal emargination (Gaffney et al., 2011), indicating more area and volume available to m. adductor mandibulae externus (Dalrymple, 1977;Sterli & De la Fuente, 2010). The mechanical effects of adductor muscles upon the lower jaw during food capture has been demonstrated in some turtles (Schumacher, 1973;Lemell, Beisser & Weisgram, 2000;Lemell et al., 2002;Pfaller, Gignac & Erickson, 2011). These studies agree that besides acting to close the mouth, the m. adductor mandibulae internus executes counter protraction forces to the m. adductor mandibulae externus retraction forces, while m. adductors mandibulae posterior produce medial forces (Figs. 9C and 9D). The contraction of all these muscles together avoid displacements of the mandible and reduce stresses at the articulation (Schumacher, 1973;Lemell, Beisser & Weisgram, 2000;Lemell et al., 2002). The positive allometries of the bones of the lower adductor chamber of B. elegans, therefore, may reflect greater resistance for a more robust musculature of m. adductor mandibulae internus and m. adductor mandibulae posterior in response to higher forces created by external adductors. Besides, these muscles also play the main role in feeding, as proposed for aquatic feeders (Frazzetta, 1968;Werneburg, 2012), in addition to a larger area between the two tips of the maxilla (i.e., SMX a = 0.70) and a flattened skull.

Feeding changes over ontogeny in B. elegans
Changes in skull shape may be due to habitat differences in which terrestrial turtles (e.g., testudinids) possess higher and shorter skulls while aquatic turtles (e.g., emydids) have flatter and longer skulls (Claude et al., 2004). The changes in skull shape of turtles along ontogeny have been assessed in living species (Dalrymple, 1977;Pfaller, Gignac & Erickson, 2011). Generally, it is supported that a diet shift occurs from small soft prey to bigger harder ones, in association with higher, larger and more robust skulls. These, in turn, are more suitable for crushing clams and/or to capture fishes by having a greater gape. The overall aquatic morphology comprises adaptations to suction feeding, which was also discussed by Herrel, O'Reilly & Richmond (2002), and could be the case of B. elegans. Firstly because taphonomic studies at Pirapozinho site suggested a riverine ephemerous system (Soares et al., 1980;Fulfaro & Perinotto, 1996;Fernandes & Coimbra, 2000;Henriques et al., 2002;Henriques et al., 2005;Bertini et al., 2006;Henriques, 2006) and fossils that experienced little transportation (Bertini et al., 2006), thus it is more likely that B. elegans was a semi-aquatic turtle, similar to the extant freshwater turtles. Secondly, the general pattern observed revealed form and shape changes in both temporal and lateral emargination (upper and lower adductor chamber, respectively): as a whole, B. elegans skull seems to become more emarginated, flattened and longer as it grows, according to the skull shape for aquatic turtles found by Claude et al. (2004), and indicating greater area and volume for muscles attachment. In addition, the deeper temporal emargination of B. elegans indicates a greater increase in muscle volume (Kischlat, 1994), thus leading to a stronger bite force (Sterli & De la Fuents, 2010). This leads us to interpret such changes as related to a shift in diet as individuals grow instead of a shift in habitat. Malvasio et al. (2003) described diet changes in Podocnemis expansa, P. unifilis and P. sexturberculata due to aging, concluding that the latter is a carnivore species, whereas the two former are omnivorous. Whereas P. expansa changes its diet becoming more herbivorous, P. unifilis remains more balanced with similar ingestion of vegetables and meat (Malvasio et al., 2003). Kischlat (1994) suggested that B. elegans might have fed of hard preys and, given the several mollusk and crustacean species described for the Pirapozinho site (Dias-Brito et al., 2001), it might have composed the diet of B. elegans. In this context, we agree with Kischlat (1994) and suggest that smaller juveniles individuals might have fed on less hard and small food items (e.g., snails and small fishes) whereas bigger old specimens fed on harder and larger preys, such as crustaceans and bigger mollusks.
Although there is a possibility that size differences could be due to sexual dimorphism as aforementioned (see Introduction, 'Geological settings and taphonomic context of the Tartaruguito site'), we were not able to assume such assumption. Furthermore, if there is size-related dimorphism, it would imply on potential diet differentiation between adults male and female of B. elegans. Since we were not able to determine size-related sexual dimorphism, such a statement is merely speculative. Romano & Azevedo (2007) (for shell material), our data did not show enough morphometrical variation to suggest population differences among our sample. Therefore, we did not have evidence to disprove that the ''Tartaruguito'' site is composed of a single population of B. elegans. However, it is feasible to assume that different generations of individuals were crowded in this locality by the accumulation of corpses due to several drying events as previously suggested by Henriques et al. (2005) and Henriques (2006). Since none B. elegans hatchling were found in the ''Tartaruguito'' site until now, it might have been preferentially a freshwater foraging area.

As in
As regards to the morphometric data, the observed variation and allometries in the skull bones, mainly the PA, QJ, SQ, QU, PO, JU, MX, PAL and PT, as well as PCAs loadings, reflect shape differences in both upper and lower adductor chambers. We interpret this allometric variation as an indicative of more area attachment and resistance for stronger adductor muscles, which are accompanied by changes in diet during aging, from softer to harder prey, as seen in living turtles species.
In regard to the use of images for carrying out morphometrics studies, we conclude that the use of calipers can be replaced by softwares that work on images. ImageJ is a useful and time-saving tool for this matter. However, one needs to beware when measuring straight lines between landmarks that are located in different depths, which result in angled lines against the projection orthogonal plane. In attention to this detail will lead to assess lower values for a given measurement than its real size.
In regard to the approaches applied to our data to deal with missing entries in the matrix (i.e., mean value or iterative imputation), both were useful for answering the questions we raised (i.e., the single population hypothesis), though little different results were obtained (few specimens out of 95% confidence ellipse in mean value approach in contrast with none specimen out of ellipse in iterative imputation approach). However, we recommend the iterative imputation as the most appropriate approach to deal with missing data in paleontological studies on the basis of the statistical assumptions it was developed (a sample-based regression for characters estimation) and the more conservative results.

AMNH
American