Material & Methods

For resolving SSD and SShD, we investigated 227 crocodile newt specimens of the genera Echinotriton and Tylototriton housed in natural history collections (Supplementary Table S1) including the following: E. andersoni from Okinawa Island, Tylototriton asperrimus, T. himalayanus, T. kweichowensis, T. shanjing, T. shanorum, T. taliangensis, T.uyenoi and T. verrucosus . The selected species represent all major clades of crocodile newts comprising the different mating modes i.e., showing a circle dance or applying an amplexus (Pogoda et al., 2020). Following species are generally thought to be circle dancers: E. andersoni, T. asperrimus, T. kweichowensis, T. shanjing , while the others are regularly observed to apply an amplexus. Tylototriton asperrimus represent the only member of the subgenus Yaotriton with a sufficient sample size of both sexes whereas unfortunately not enough female specimens of other species were available in natural history collections due to a heavy male biased field sampling during the breeding season. To access osteology for SD analyses, specimens were µCT-scanned. CT-scans were carried out either with a Bruker SkyScan1272 with the software NRecon (Bruker CT) for reconstructions or within the X-ray imaging laboratory at the Institute for Photon Science and Synchrotron Radiation, Karlsruhe Institute of Technology (KIT) employing a microfocus x-ray tube (XWT-225, X-RAY WorX, Garbsen, Germany) and a flat panel detector (XRD 1621 CN14 ES, PerkinElmer, Waltham, USA) in combination with a custom designed mechanical sample manipulator. For the CT scans made at KIT, Octopus 8.6 (Inside Matters, Gent, Belgium) was used to perform the tomographic reconstruction. The scan resolution was either 20.1 (SkyScan) or 21.3 µm (KIT-custom build scanner).
To catch the entire shape variation of the cranium 45 three-dimensional (3D) landmarks were digitized and for the analysis of the humerus shape six fixed landmarks and 50 semi-landmarks in three curves were digitized (Fig.1). Prior to landmark digitization, potential error in setting landmarks was validated by digitizing one specimen five times and five additional specimens of the same species to compare consistent placement by Procrustes distance of the respective mean shapes. Landmark digitization was carried out by one author with the software Checkpoint v.2019.03.04.1102 (Stratovan Ltd.). Geometric morphometrics was performed in R version 3.6.3 (R Development Core Team, 2019) using the packages geomorph v.3.2.1, RRPP v. 0.5.2 and Morpho 2.8 (Schlager, 2017, Collyer & Adams, 2018, Adams et al., 2019). Complete landmark configurations are a prerequisite for GM analyses. Hence, missing landmarks (due to anomalies or injuries) were first estimated by thin plate spline approach implemented in the function ‘estimate.missing’ (Gunz et al., 2009). Semi-landmarks in the humerus dataset were equally spaced along the digitized curve. Variation due to location, rotation and scale was removed by a generalized Procrustes alignment (GPA) using the function ‘gpagen’ (Rohlf & Slice, 1990). In the humeri dataset, semi-landmarks were simultaneously slided using minimized bending energy (Bookstein, 1997a, Perez et al., 2006). As asymmetry was not in the scope of this study, bilateral landmarks in cranial landmark configuration were symmetrized by averaging left and right landmark pairs. Skulls and humeri were analysed further in the same approach. A principal component analysis (PCA) on Procrustes coordinates was performed and plotted to investigate general shape variation. To account for size, we used logarithm of centroid size (CS), which represents a measure of size in GM (Bookstein, 1997b, Zelditch et al., 2012).
A full factorial model design including species and mating mode was precluded by the model system as each species comprises only a single mating mode. Thus, several Procrustes ANOVAs had to be performed to investigate all potential sources of morphological variation. First, a Procrustes ANOVA as implemented in the function ‘procD.lm’ with size, species and sex including all interactions was performed. Allometry between sexes was not different, indicated by non-significant interaction between sex and CS. Thus, we explored allometric shape change by another Procrustes ANOVA including CS and sex only. To test whether mating mode affects the pattern of SD, we first ran a model design including only sex, mating mode and its interaction, and second a model including these factors plus CS as covariate. To explore different patterns of SShD, we performed a trajectory analysis to visualize shape change directions between species and performed a group mean prediction with 95%-confidence intervals for males and females in each species, implemented in the function ‘predict.lm.rrpp’ in the RRPP package. Sexual size dimorphism patterns between species and mating modes were estimated by a Procrustes ANOVA of species and sex on CS and, in a second one, sex and mating mode as variables. The function ‘pairwise’ was used to reveal which groups were different. According to the model, a grouping variable of sex with species or mating mode was used. Significance testing was performed using Residual Randomization by 10.000 random permutations (Collyer et al., 2015, Collyer & Adams, 2018). Shape changes were visualized as TPS-grids by warping the mean shape by thin-plate spline approach with the function ‘plotRefToTarget’.