Signatures of echolocation and dietary ecology in the adaptive evolution of skull shape in bats

Morphological diversity may arise rapidly as a result of adaptation to novel ecological opportunities, but early bursts of trait evolution are rarely observed. Rather, models of discrete shifts between adaptive zones may better explain macroevolutionary dynamics across radiations. To investigate which of these processes underlie exceptional levels of morphological diversity during ecological diversification, we use modern phylogenetic tools and 3D geometric morphometric datasets to examine adaptive zone shifts in bat skull shape. Here we report that, while disparity was established early, bat skull evolution is best described by multiple adaptive zone shifts. Shifts are partially decoupled between the cranium and mandible, with cranial evolution more strongly driven by echolocation than diet. Phyllostomidae, a trophic adaptive radiation, exhibits more adaptive zone shifts than all other families combined. This pattern was potentially driven by ecological opportunity and facilitated by a shift to intermediate cranial shapes compared to oral-emitters and other nasal emitters.

Curves with equidistance sliding semi-landmarks: 1. From L27 to L28, along the dorsal midline of the cranium (blue) 2-3. From L29/30 to L31/32, along the dorsal profile of the zygomatic arch (red) [4][5]along Table 1: Number of critical principal components per anatomical structure (cranium or mandible) and per clade (all bats or phyllostomids), with (pPCA) or without (PCA) phylogenetic correction. PC scores from these axes were used in all subsequent macroevolutionary model fitting analyses.

Methods:
Geometric morphometric methods are currently intolerant of missing data. We follow the approach of Arbour and Brown 1 to evaluate the effect of missing landmark coordinate data and specimen incompleteness on the resulting morphospaces of the cranium and mandible. Namely, we evaluated whether it was better to estimate data missing from incomplete specimens, or to exclude those specimens entirely. While both datasets have a low overall proportion of missing data, this is spread across a large number of specimens and landmarks (e.g., ~40% of crania and 78% of cranial landmarks have at least 1 missing value). Therefore, the exclusion of missing data would substantially limit our dataset, especially for rare taxa with lower sample sizes. To test for the impact of missing data estimation, we used the following steps (and illustrated in workflow below). 1) We selected from our data a subset representing all specimens with complete landmark sets as our training dataset (dataset COM). We assume that this dataset better reflects the variance associated with landmark positions across bats than either datasets with missing data estimated, or incomplete specimens excluded. 2) From a chosen proportion of specimens, we randomly sampled landmarks based on the distribution of missing landmarks from our real incomplete specimens, and removed the corresponding coordinate data. Missing values were subsequently estimated using either "reflected relabeling" or BPCA (see below), and this formed the dataset EST. 3) We then excluded all simulated incomplete specimens (dataset EXC). 4) We evaluated the change in the variance structure of the dataset imparted by estimation or exclusion. We used the R function "procrustes" (from the package "vegan") to quantify the fit between the eigenvectors of COM vs. EST and COM vs. EXC as calculated by the "Procrustes sum of squares" (PSS). Lower values of PSS indicate a stronger fit between datasets.
For step 2, we first exploited bilateral symmetry across the cranium and mandible. Missing bilaterally symmetrical landmarks can be estimated using "reflected relabeling" 2 , wherein the specimen is mirrored and aligned with the (overlapping) original landmarks using Procrustes superimposition; missing landmarks on either side are then imputed from the mirrored specimen. We applied reflected relabeling to the landmark data using the R function "flipped" from the package LOST 3 . For those missing landmarks that could not be estimated using reflected relabeling (3.2% and 0.4% of all specimen cranial and mandibular landmarks and semi-landmarks, respectively), we employed a Bayesian Principal Component Analysis approach 4 , using the R function "MissingGeoMorph" from the package "LOST" 1,3,5 .
The preceding four steps were carried out on each of the cranium and mandible datasets across a range of values for the proportion of incomplete specimens, based on that observed in each of our datasets (2% to 20% of mandibles and 30 to 50% of all crania). We simulated 500 datasets with missing data per value of the proportion of incomplete specimens. This helps to prevent the evaluation of axes with eigenvalues lower than from random datasets, and less likely to represent useful shape variation. We also carried out PCA both with (pPCA) and without (PCA) phylogenetic correction. To improve comparability between PSS values across PCA and estimation methods (and because the EXC datasets contains fewer eigenvectors than the COM or EST datasets), we used the eigenvectors for the number of critical axes observed in the main results (3 PC axes from the cranium, 4 PC axes from the mandible).
Additionally, we tested the relative impact of estimating missing data by exploiting bilateral symmetry from each individual specimen, or by comparisons with the full dataset. For each of the cranium and mandible we created two datasets with missing landmarks, a) one with landmarks missing from only one side (to evaluate reflected relabeling, see methods) and b) one with landmarks missing from either side (to evaluate BPCA, see methods). We selected BPCA for non-paired missing landmarks as it has shown to reliably estimate missing data from a variety of multidimensional morphometric datasets, across a range of sample sizes and taxonomic scopes 1,5 .

Results:
Across both the cranium and mandible, the impact of missing data estimation was less than the impact of excluding incomplete specimens at approximating the shape variance associated with all fully complete specimens. At low proportions of missing specimens, the PSS generated by the EST and EXC mandible datasets overlapped (Supplementary Figure 3), however at 2% incomplete mandibles (which showed the greatest overlap) the number of individual datasets for which the PSS of EXC was lower than the PSS of EST was only 0.2%. Additionally, we would expect the results of EXC and EST to converge as the proportion of incomplete specimens decreases, as both datasets contain almost all specimens present in COM. The median value of PSS using BPCA was consistently lower than the EXC dataset for the cranial dataset. However, the absolute range of PSS also overlapped when missing data in the cranial dataset was estimated using BPCA. Similar to the mandible dataset, across most individual datasets estimation was preferred over exclusion (e.g., at 40% incomplete specimens, only 2.8% of simulated datasets showed lower PSS values for EXC than EST, and at 30% incomplete specimens, only 3.8% of simulated datasets showed lower PSS values for EXC than EST). "Reflected relabeling" performed better than BPCA across both datasets, likely as a result of better incorporating intra-individual and intra-specific shape variation. A large proportion of missing data in our original datasets represented landmarks that could be estimated by exploiting bilateral symmetry (see main text methods), and our initial selection of damaged skulls was generally biased towards such specimens. Thus, our analysis of the impact of BPCA is conservative, as markedly fewer landmarks needed to be estimated this way than represented in our simulated datasets.
These results are consistent with previous analyses of the impact of incomplete specimens in geometric morphometric analyses, which showed that for most estimation approaches it was generally better to estimate missing data than exclude incomplete specimens. Supplementary Figure 2: Results of missing data analyses of cranium morphology. The fit between the eigenvectors, of either PCA (a and c) or pPCA (b and d), of the COM dataset (complete specimens) versus either the EST dataset (missing data estimated) or the EXC dataset (incomplete specimens excluded) were determined using Procrustes-sum-of-squares (PSS) from the R function "procrustes". Missing data was estimated using reflected relabeling (a-b) or BPCA (c-d). Lines indicate the median and shaded area indicate the 95% range of 500 simulated incomplete datasets. See source data file.
Supplementary Figure 3: Results of missing data analyses of mandible morphology. The fit between the eigenvectors, of either PCA (a and c) or pPCA (b and d), of the COM dataset (complete specimens) versus either the EST dataset (missing data estimated) or the EXC dataset (incomplete specimens excluded) were determined using Procrustes-sum-of-squares (PSS) from the R function "procrustes". Missing data was estimated using reflected relabeling (a-b) or BPCA (c-d). Lines indicate the median and shaded area indicate the 95% range of 500 simulated incomplete datasets. See source data file.