Differences in jaw function experienced through ontogeny can have striking consequences for evolutionary outcomes, as has been suggested for the major clades of mammals. In contrast to placentals, marsupial newborns have an accelerated development of the head and forelimbs, allowing them to crawl to the mother’s teats to suckle within just a few weeks of conception. The different functional requirements that marsupial newborns experience in early postnatal development have been hypothesized to have constrained their morphological diversification relative to placentals. Here, we test whether marsupials have a lower ecomorphological diversity and rate of evolution in comparison to placentals, focusing specifically on their jaws. To do so, a geometric morphometric approach was used to characterize jaw shape for 151 living and extinct species of mammals spanning a wide phylogenetic, developmental, and functional diversity. Our results demonstrate that jaw shape is significantly influenced by both infraclass and diet, with substantial ecomorphological convergence between metatherians and eutherians. However, metatherians have markedly lower disparity and rate of mandible shape evolution than observed for eutherians. Thus, despite their ecomorphological diversity and numerous convergences with eutherians, the evolution of the jaw in metatherians appears to be strongly constrained by their specialized reproductive biology.

Script_jaw_mammals_paper.r : it is the R script containing all the analyses (can be open in R or RStudio or any text reader)

Script_bayestrait_BM.txt: it is the bayestrait script used in this study (can be open using any text reader)

Ident_tout_mean.csv: identification for the species names, infraclass, diet used in the analyses (can be open in excel, R, text reader)

means.Rdata: procrustes coordinates coordinates for the entire dataset (can be open in R and Rstudio)

means_extant.Rdata:procrustes coordinates coordinates for extant dataset (can be open in R and Rstudio)

phy.Rdata: phylogeny for the entire dataset (can be open in R and Rstudio)

phy_extant.Rdata: phylogeny for the extant dataset (can be open in R and Rstudio)

(a) Material

Our data set is composed of the mandible from 151 individuals representing 52 modern and 3 fossil metatherian species and 75 modern and 20 fossil eutherians. This sample was chosen in order to represent a wide phylogenetic breadth, as well as diversity in ecology, morphology, and function among terrestrial eutherians and metatherians. One hundred and thirty meshes were generated for this study and 21 were collected from online repositories (electronic supplementary material, table S1). All information about the specimens and the definition of their dietary and reproductive mode came from the literature (see electronic supplementary material, table S1), especially for extinct taxa. Furthermore, we consider that all metatherians share the reproductive mode observed in extant marsupials, although this is ambiguous for stem taxa (e.g. borhyaenids). Dietary categories were defined as carnivorous, lingual feeder, browser, insectivorous, grazer, tuberivores, frugivorous an omnivorous following the literature.

(b) Quantification of mandibular shape using three dimensional geometric morphometrics

A total of 16 landmarks and 98 curve sliding semi-landmarks were identified to comprehensively capture the shape of the mandible (electronic supplementary material, figure S1 and tables S2-3). All landmarks and curve semilandmarks were taken manually by the same person (A.-C. F.) using Checkpoint (Stratovan, Davis, CA, USA). In order to transform the curve semilandmarks into geometrically homologous points, a three-dimensional sliding semilandmark procedure was performed following the protocol described in previous studies. In order to slide all of the curve semilandmarks while minimizing bending energy, we used the ‘slider3d’ function from the Morpho R package. Finally, a Procrustes superimposition was performed using the ‘gpagen’ function from the geomorph R package to remove the effects of non-biological variation (rotation, translation, and isometric size).

(c) Phylogenetic tree

A phylogenetic tree was constructed to incorporate all of the species used in this study, based on a recently dated molecular tree for extant mammals. Fossils were grafted onto this tree based on recent morphological studies (see electronic supplementary material), with tip positions informed by their occurrences.

(d) Shape variation depending on infraclass and diet

Both non-phylogenetic and phylogenetic principal components analysis (PCA) were used to visualize jaw shape variation across eutherians and metatherians. To do so, we used respectively the functions ‘gm.prcomp’ from the R geomorph package, and ‘phyl.pca’ and ‘phylomorphospace’ from the R phytools package. Next, we tested for shape differences between infraclass and diet using a type II phylogenetic multivariate analysis of variance (MANOVA) in R mvMORPH package. The multivariate phylogenetic linear models were fitted with a Pagel’s lambda by penalized likelihood using the ‘mvgls’ function. Pagel’s lambda has the advantage to provide increased flexibility in estimating the error structure and it is equivalent to fitting a phylogenetic mixed model while accounting for departure from Brownian Motion. Subsequently, this model was used as input in the ‘manova.gls’ function (electronic supplementary material, table S4). The significance of each type II phylogenetic MANOVA was assessed using a Pillai statistic and 1000 permutations. Finally, to identify morphological convergence depending on diet and infraclass, we used the phylogenetic ridge regression method of the RRphylo package in R. We first performed a phylogenetic ridge regression on the tree and with the first 24 PC scores (accounting for 95% of the overall variance) using the ‘RRphylo’ function in order to obtain the branch-wise evolutionary rates and the ancestral character estimates at nodes. Next, we establish morphological convergence for each diet category using the ‘search.conv’ function under “state” cases on PC scores. From these analyses, we retrieve the mean angle (mean angle between species within the dietary category), and the p of the mean angle (significance level for mean angle).

(e) Disparity based on infraclass

To assess and compare disparity between eutherian and metatherian reproductive modes, we calculated Procrustes disparity using the ‘morphol.disparity’ and the ‘dispRity.per.group’ functions respectively in the Geomorph and the DispRity packages in R. Procrustes disparity allows to estimate the Procrustes variance per group by using the residuals of a linear model fit. Here, we used as input the Procrustes coordinates in order to calculate the sum of the diagonal elements of the infraclass group covariance matrix divided by the number of observations in the infraclass group. To test for Procrustes disparity difference between infraclass, a pairwise comparison was performed using the ‘morphol.disparity’ function in the Geomorph package in R. We further used a non-parametric Wilcoxon test to assess the significance of differences in disparity between eutherians and metatherians using the ‘wilcox.test’ function in the DispRity package in R. Disparity analyses were performed on the entire dataset (extinct and extant species) and on a dataset comprising only modern mammals.

(f) Rates of morphological evolution depending on infraclass

In order to assess and compare morphological rates of evolution between infraclass, we calculated evolutionary rates for each infraclass category for both entire and extant species data sets. To do so, the reconstructed history between eutherian and metatherian reproductive modes categories on which a state-specific Brownian motion (BMM) model were fitted was obtained through stochastic character mapping across a sample of 100 trees using the ‘make.simmap’ function and an ‘ARD’ model in the R package phytools. Model fit was performed using a state-specific Brownian motion model in the ‘mvgls’ function in mvMORPH R package 1.1.4 both using the entire data set and on the modern data set only. Finally, we estimated a branch-specific evolutionary rates and rate shifts using the variable rates model implemented in BayesTraitsV3 (http://www.evolution.rdg.ac.uk/). Shifts in the rate of continuous trait evolution (modelled by a Brownian motion process) were detected using a reversible-jump Markov chain Monte Carlo algorithm. The phylogenetic principal components (PCs) accounting for 95% of the overall variation in jaw shape were used as input (i.e., the first 33 phylogenetic PCs for the entire data set and the first 32 for the extant data set). Ten independent chains were run for 200,000,000 iterations, sampling every 10,000 iterations, and the first 25,000,000 iterations were discarded as burn-in. Trace plots were examined and two independent chains that were stationary after burn-in were kept. We assessed the effective sample size (ESS) of the posterior samples (ESS > 100) as well as the convergence of the chains, with a Gelman and Rubin’s convergence diagnostic, using respectively the functions ‘effectiveSize’ and ‘Gelman.diag’ implemented in the R coda package (see electronic supplementary material, figure S2-3 and tables S4-7). Finally, we plotted the results of the analyses on the tree using the function ‘mytreebybranch’ (https://github.com/anjgoswami/salamanders/blob/master/mytreerateplotter.R). The branch-specific average rate and the posterior probability of rate shifts were summarized from the posterior samples using the ‘rjpp’ and ‘plotShift’ functions in the btrtools R package v.0.0.0.9000 (https://github.com/hferg/btrtools/tree/master/R).