Sea snakes in the Hydrophis-Microcephalophis clade (Elapidae) show exceptional body shape variation along a continuum from similar forebody and hindbody girths, to dramatically reduced girths of the forebody relative to hindbody. The latter is associated with specialisations on burrowing prey. This variation underpins high sympatric diversity and species richness and is not shared by other marine (or terrestrial) snakes. Here, we examined a hypothesis that macroevolutionary changes in axial development contribute to the propensity, at clade level, for body shape change. We quantified variation in the number and size of vertebrae in two body regions (pre- and post-apex of the heart) for ~94 terrestrial and marine elapids. We found Hydrophis-Microcephalophis exhibit increased rates of vertebral evolution in the pre- versus post-apex regions compared to all other Australasian elapids. Unlike other marine and terrestrial elapids, axial elongation in Hydrophis-Microcephalophis occurs via the preferential addition of vertebrae pre heart apex, which is the region that undergoes concomitant shifts in vertebral number and size during transitions along the relative fore- to hindbody girth axis. We suggest that this macroevolutionary developmental change has potentially acted as a key innovation in Hydrophis-Microcephalophis by facilitating novel (especially burrowing) prey specialisations that are not shared with other marine snakes.

Vertebral number: samples and measurements

We sampled 275 alcohol-preserved adult specimens representing 94 species of terrestrial (n=51), fully marine (n=40) and semi-aquatic (n=3) ecologies from the Elapidae family (Supp. Table S1 of paper). Species were chosen from the Australasian region to capture the diversity in number of pre-cloacal vertebrae. Details of specimens studied and the museums from which they were sourced are given in Table S1 of the paper. Sex was not considered in this study, but tail vertebrae were omitted because they show the most sexual dimorphism. Data for marine species was taken from [12]. Counts of vertebrae were made from ventral scales for the terrestrial species, and Aipysurus + Emydocephalus species group, since they display a 1:1 ratio for pre-cloacal vertebrae [23]. Heart position was found through small ventral incisions and noted based upon vertebral number counted from behind the head to posterior-most point of heart. Other marine species were examined using X-rays since they do not display the ventral scale correspondence (details of X-ray below). Metal pins were placed into the preserved specimens at the posterior-most point of the heart and the cloaca to be visible in X-rays (example shown in Supp. Figure S1 of paper). Counts started at the first vertebra attaching to the skull and proceeded to the last vertebra anterior to the cloaca.

Digital X-rays were made using two systems: small specimens were imaged using the Faxitron LX-60 machine at University of Adelaide Health and Medical Sciences facility; and for large specimens we used a Siemens Multix Fusion Max machine at Dr. Jones and Partners Medical Imaging, Adelaide (Supp. Figure S1 of paper). Both systems save the X-rays as a DICOM format with size embedded into the file, such that measurements are scaled automatically. 

Absolute body size and number of pre-cloacal vertebrae differ greatly among the elapid species we sampled; for example, among the smallest species is a terrestrial/fossorial species, Simoselaps bertholdi, at 23.5 cm neck to cloaca with 124 pre-cloacal vertebrae, and among the largest is a terrestrial species, Oxyuranus scutellatus at 146.5 cm and 237 vertebrae. But length and number of vertebrae is not always correlated in elapid snakes [23, 24], for example, Toxicocalamus preussi is one of the smaller species at 55 cm with the most pre-cloaca vertebrae in the dataset (325). Since snakes have indeterminate growth, and museum samples may be skewed to smaller individuals, absolute body size was not considered further; instead relative vertebra size was used (actual vertebra size divided by body length, see below), which is known to remain consistent during ontogeny in some elapids [13] (but see [25, 26]). 

Intracolumnar vertebral size: samples and measurements

To capture variation in intracolumnar vertebral size we sampled one adult of each species for a subset of 61 species of terrestrial (n=27), fully marine (n=31) and semi-aquatic (n=3) elapid snakes (Table S2 of paper). Sampling one individual per species is sufficient in this instance because there is no appreciable difference in vertebral column intracolumnar profiles among adults within species (Supp. Fig. S2) nor during postnatal growth [13]. Vertebrae size was examined by measuring the length of every vertebra from the first pre-cloacal vertebra after the atlas to the vertebra anterior to the cloaca [12]. Tail vertebrae were not examined in this study because they are known to exhibit sexual dimorphism [e.g., 27]. We measured the length of each vertebra from the X-rays using the ‘multipoint tool’ in ImageJ v.1.52i [28]: landmarks were placed medially along the vertebral column at the anterior limit of the centrum of each vertebra. Coordinates (x,y) of the landmarks in millimetres were exported into the R statistical environment v.4.0.5 [29] and inter-landmark distances were calculated by applying the Pythagorean theorem between sequential coordinate points. Vertebral width could not be measured due to the changing orientation of the vertebral column (axial torsion) resulting from specimen preservation. 

Phylogenetic hypothesis

We built a consensus phylogenetic tree of elapids to perform a comparative analysis. Molecular data were obtained for 183 species of elapid using the mitochondrial 12S, 16S, ND4 and cytochrome b genes and the nuclear C-mos, RAG-1 and RAG-2 genes. The bulk of the alignment was obtained from [11] however we removed Aipysurus pooleorum as the GenBank sequences were a composite of A. pooleorum and A. foliosquama. We added Acanthophis laevis, Acanthophis pyrrhus, Hoplocephalus bungaroides, Neelaps bimaculatus, Pseudonaja inframacula, Pseudonaja nuchalis from GenBank sequences. The sea snake taxa Aipysurus tenuis and Emydocephalus orarius were added from supplementary material in [30] while new ND4 sequence was obtained for Demansia reticulata (Genbank number to be added) using the protocols in [31]. Additional taxa were aligned and checked by eye in Geneious Prime v. 2022.0.1 (https://www.geneious.com) resulting in a final alignment of 9,078 bp with 183 terminal taxa.

We initially used the partitioning scheme of [11] however due to issues with reaching convergence we used simpler models of sequence evolution. The final partitioning scheme consisted of the following partitions: i) nuclear coding regions, codons 1+2 – HKYig; ii) nuclear coding regions, codon 3 – HKYig; iii) mitochondrial coding regions, codon 1 – HKYig; iv) mitochondrial coding regions, codon 2 – HKYig; v) mitochondrial regions, codon 3 – HKYg; vi) 12S rRNA – HKYig; and vii) 16S rRNA – HKYig. Phylogeny and dates were reconstructed using BEAST v. 2.6.6 with dates calibrated using the same nodes as [11]. Clock and tree models were linked with a strict clock and a Yule tree model was selected. The Markov Chain Monte Carlo (MCMC) was run for 10,000,000 generations with trees sampled every 1000 states. Convergence was checked using Tracer v.1.7.2 [32] and a burn in of 25% and Effective Sample Size (ESS) values of >100 were reached for most parameters before burn in. Maximum clade credibility trees were produced from the remaining 7,500 trees using TreeAnnotator v.2.6.6 [33]. The tree was pruned to 90 taxa included in this study using ‘drop.tip’ function in ape R package v.5.6-2[34], where Denisonia maculata was substituted for Denisonia devisi. Taxa without a suitable substitution were omitted from comparative analyses: Brachyurophis fasciolatus, Hoplocephalus stephensi, Salmonelaps par, and Hydrophis melanosoma.

We divided the tree into four groups based upon ecological niche and phylogenetic relatedness [add citations where this info is from]: terrestrial species; the monophyletic group of marine species (Aipysurus & Emydocephalus); the semi-aquatic species (Hydrelaps darwiniensis, Ephalophis greyae, Parahydrophis mertoni); and the other monophyletic group of marine species (Hydrophis & Microcephalophis). Note that the semi-aquatic species render the marine clade paraphyletic. 

References:

12. Sherratt E., Coutts F.J., Rasmussen A.R., Sanders K.L. 2019 Vertebral evolution and ontogenetic allometry: The developmental basis of extreme body shape divergence in microcephalic sea snakes. Evol Dev 21 135–144. (doi:10.1111/ede.12284).

13. Sherratt E., Sanders K.L. 2020 Patterns of intracolumnar size variation inform the heterochronic mechanisms underlying extreme body shape divergence in microcephalic sea snakes. Evol Dev 22(3), 283-290. (doi:10.1111/ede.12328).

23. Voris H.K. 1975 Dermal scale-vertebra relationships in sea snakes (Hydrophiidae). Copeia 4, 746–757.

24. Lindell L. 1994 The evolution of vertebral number and body size in snakes. Funct Ecol 8(6), 708-719.

27. King R.B. 2008 Sexual dimorphism in snake tail length: sexual selection, natural selection, or morphological constraint? Biol J Linn Soc 38(2), 133-154. (doi:10.1111/j.1095-8312.1989.tb01570.x).

28. Schneider C.A., Rasband W.S., Eliceiri K.W. 2012 NIH Image to ImageJ: 25 years of image analysis. Nat Meth 9(7), 671-675. (doi:10.1038/nmeth.2089).

29. R Development Core Team. 2022 R: a language and environment for statistical computing. v.4.0.5

30. Nitschke C.R., Hourston M., Udyawer V., Sanders K.L. 2018 Rates of population differentiation and speciation are decoupled in sea snakes. Biol Lett 14(10), 20180563. (doi:10.1098/rsbl.2018.0563).

31. Sanders K.L., Lee M.S., Mumpuni, Bertozzi T., Rasmussen A.R. 2013 Multilocus phylogeny and recent rapid radiation of the viviparous sea snakes (Elapidae: Hydrophiinae). Mol Phylogenet Evol 66(3), 575-591. (doi:10.1016/j.ympev.2012.09.021).

32. Rambaut A., Drummond A.J. 2007 Tracer. v.1.7.2

33. Drummond A.J., Rambaut A. 2007 BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7(1), 214.

34. Paradis E., Claude J., Strimmer K. 2004 APE: Analyses of phylogenetics and evolution in R language. Bioinformatics 20(2), 289-290.

R statistical environment v.4.0.5 Libraries

• library(ape) # v. 5.6-2
• library(phytools) # 1.0-3
• library(geomorph) # 4.0.4