Evolution of ossification sequences in salamanders and urodele origins assessed through event-pairing and new methods

Ossification sequences of the skull in extant Urodela and in Permo-Carboniferous Branchiosauridae have already been used to study the origin of lissamphibians. But most of these studies did not consider some recent methods developed to analyze the developmental sequences within a phylogenetic framework. Here, we analyze the ossification sequences of 24 cranial bones of 23 extant species of salamanders using the event-pairing method. This reveals new developmental synapomorphies for several extant salamander taxa and ancestral sequences for Urodela under four alternative reference phylogenies. An analysis with the 12 bones for which ossification sequence data are available in urodeles and in the branchiosaurid Apateon is also performed in order to compare the ancestral condition of the crown-group of Urodela to the sequence of Apateon. This reveals far more incompatibilities than previously suggested. The similarities observed between some extant salamanders and branchiosaurids may result from extensive homoplasy, as the extreme variation observed in extant Urodela suggests, or be plesiomorphic, as the conservation of some ossification patterns observed in other remotely related vertebrates like actinopterygians suggests. We propose a new, simpler method based on squared-change optimization to estimate the relative timing of ossification of various bones of hypothetical ancestors, and use independent-contrasts analysis to estimate the confidence intervals around these times. Our results show that the uncertainty of the ancestral ossification sequence of Urodela is much greater than event-pairing suggests. The developmental data do not allow to conclude that branchiosaurids are closely related to salamanders and their limited taxonomic distribution in Paleozoic taxa precludes testing hypotheses about lissamphibian origins. This is true regardless of the analytical method used (event-pairing or our new method based on squared-change parsimony). Simulations show that the new analytical method is generally more powerful to detect evolutionary shifts in developmental timing, and has lower Type I error rate than event-pairing. It also makes fewer errors in ancestral character value or state assignment than event-pairing.