Evolutionary convergence of muscle architecture in relation to locomotor ecology in snakes

Abstract The epaxial muscles in snakes are responsible for locomotion and as such can be expected to show adaptations in species living in different environments. Here, we tested whether the structural units that comprise the superficial epaxial muscles (semispinalis‐spinalis, SSP; longissimus dorsi, LD; iliocostalis, IC) were different in animals occupying similar habitats. To do so, we analyzed and compared the muscle architecture (mass, fiber length, and physiological cross‐sectional area) of the superficial epaxial muscle segments in snakes that differ in their habitat use (e.g., arboreal, terrestrial, and aquatic). Our results showed that arboreal species have on average longer muscles and tendons spanning more segments likely important during gap bridging. Moreover, aquatic snakes show relatively heavier semispinalis‐spinalis muscles with a greater cross‐sectional area. The longissimus dorsi muscles also showed a greater cross‐sectional area compared with terrestrial and especially arboreal snakes. Whereas the more strongly developed muscles in aquatic snakes are likely associated with the dense and viscous environment through which they move, the lighter muscles in arboreal snakes may provide an advantage when climbing. Future studies comparing other ecologies (e.g., burrowing snakes) and additional muscle units (e.g., multifidus; hypaxial muscles) are needed to better understand the structural features driving variation in locomotor performance and efficiency in snakes.

undulation is used in most environments, concertina, rectilinear locomotion, and sidewinding are more specialized locomotor modes used in specific contexts (Jayne, 2020). Previous studies have suggested that the most superficial epaxial muscles, namely semispinalis-spinalis (SSP), longissimus dorsi (LD), and iliocostalis (IC) muscle groups are the main drivers of locomotion (Gans, 1986;Gasc, 1981;Gasc et al., 1989;Jayne, 1988aJayne, , 1988bMoon & Gans, 1998;Mosauer, 1932). These muscles represent roughly 65% of the axial muscle mass (Ruben, 1977) and their cross-section scales isometrically with overall body mass (Moon & Candy, 1997). Yet, variation in the mass and length of these muscles has been documented and has been suggested to be related to the locomotor environment and predation mode (Herrel et al., 2011;Jayne, 1982;Lourdais et al., 2005;Moon, 2000;Penning, 2018;Penning & Moon, 2017;Ruben, 1977;Young, 2010). However, only a few studies have quantitatively explored differences in the muscle architecture of the epaxial musculature (Jayne, 1982;Penning, 2018) and even fewer have done so in a comparative phylogenetic framework. Here, we test whether the structural units that make up the superficial epaxial muscles (SSP, LD, IC) have evolved convergently in species subjected to similar locomotor constraints. We analyzed and compared the muscle architecture (mass, fiber length, and physiological cross-sectional area) of the muscle segments of the superficial epaxial muscles (SSP, LD, IC) in snakes that differ in their habitat use (i.e., arboreal, terrestrial, and aquatic). We specifically focus on aquatic species given that (1) a significant number of species from different families have convergently evolved an aquatic or semi-aquatic lifestyle (Murphy, 2007;Pauwels et al., 2008) and (2) the strong physical constraints (drag) associated with the aquatic medium (i.e., viscosity and density) that will impact the kinematics and energetics of locomotion (Vogel, 1994). We predict that aquatic species should have stronger muscles to move the body against a dense fluid.
Conversely, we predict arboreal species to have longer and lighter muscles with longer tendons allowing them to bridge gaps in the complex three-dimensional environment they live in (Jayne, 1982;Jayne & Riley, 2007;Lillywhite et al., 2000).
Species were selected to encompass a diversity of life-styles across the different families. We were specifically interested in adaptations to an aquatic life-style and focused our sampling on this ecology. We classified semi-aquatic species as aquatic in our analyses as we had too few semi-aquatic species to treat them as a separate group.

| Dissection
All individuals were weighed (± 1 g; Kern PNJ). We measured snoutvent length and tail length with a string and a ruler to the nearest mm. Body diameter at mid body was measured using manual Vernier calipers (Hilka Tools 150 mm, ± 0.02 mm). Three sections of equal size (from 1 to 10 cm depending on the size of the snake) were analyzed given that snakes are known to exhibit longitudinal variation in the anatomy of the epaxial muscles (Pregill & Pregill, 1977). We dissected the epaxial muscles at (1) the anterior body at about one fourth of the total length of the snake, (2) midbody, and (3) the posterior body just anterior to the cloaca.
We measured the longitudinal extent (number of vertebrae, including the vertebrae of origin and insertion for muscles or tendons inserting directly on the vertebrae) and the size (in cm, with a ruler) of the upper (i.e., anterior tendon of the semispinalis-spinalis; posterior tendon of the iliocostalis) tendons, the muscle bellies, and the lower (posterior tendon of the semispinalis-spinalis; anterior tendon of the iliocostalis) tendons of the semispinalis-spinalis, the longissimus dorsi, and the iliocostalis muscles ( Figure 1). Note that the posterior tendons of the semispinalis-spinalis and the iliocostalis are attached to the longissimus dorsi. By consequence, the posterior (lower) tendon of the semispinalis-spinalis is also the upper tendon of the longissimus and the posterior (upper) tendon of the iliocostalis is also the lower tendon of the longissimus dorsi.
Next, muscles were removed and weighed on a precision balance (Mettler AE100, precision ± 0.0001 g). We then measured the insertion angle of the fibers on the tendons by using a binocular scope (Leica Wild M3Z) with camera lucida (Motic MLC-150C) and a protractor to obtain the pennation angle. Next, we submerged the muscles in a 30% aqueous nitric acid solution for 24-48 h. After digestion of the connective tissue, fibers were teased apart and the nitric was replaced by a 50% aqueous glycerol solution to stop further digestion. The fibers were then viewed, redrawn under a binocular scope (Leica Wild M3Z) with camera lucida (Motic MLC-150C), photographed and measured in ImageJ (Abramoff et al., 2004).
Muscle cross-sectional area was measured by first calculating muscle volume from muscle mass using a muscle density of 1.06 gcm −3 (Mendez & Keys, 1960) and then dividing this by fiber length while correcting for pennation angle by multiplying this by the cosine of the pennation angle.

| Statistical analyses
All data were Log 10 -transformed before analysis. First, we performed a descriptive analysis on the segment lengths to explore variation in the extent of the muscle-tendon units by exploring variation in the number of vertebrae spanned by each muscle and tendon. Next, we ran four separate principal component analysis (PCA) in R (packages FactomineR, FactoExtra) on the muscle architecture data. One for the overall data including snout-vent length, body mass, mid-body diameter, and the overall mass and length of the three muscles. Next, we ran separate PCAs in each of the three body regions (anterior, mid-body, and posterior) including absolute values of the mass, fiber length, and physiological cross-sectional area of each muscle segment. We cropped the phylogeny provided by Pyron et al. (2013) to include only the species included in our data set ( Figure 2) and used it to create a phylomorphospace (Revell & Collar, 2009) in Phytools (Revell, 2012). We next calculated the phylogenetic signal in our data. To do so we calculated a multivariate Blomberg's K (Kmult) in the package Geomorph (Adams & Otárola-Castillo, 2013).
We subsequently performed a MANCOVA (multivariate analysis of covariance) on the muscle architecture values (fiber length, muscle length, muscle mass, PCSA) for the three body parts (anterior, mid-body, posterior) separately to test whether species with different locomotor ecologies differed in their muscle architecture. Given that significant differences were observed, we then ran univariate analyses of variance to test which variables differed between groups. Next, we calculated unstandardized residuals of a regression of each muscle variable on snout-vent length and used them as input for ANOVAs coupled to Tukey posy-hoc tests. For variables showing significant differences, we TA B L E 1 Species used, including the number of individuals dissected, the ecological classification used, the family, snout-vent length (SVL), body mass, and mid-body diameter.

F I G U R E 1
Photograph of the axial muscles in Pantherophis guttatus. 1, upper tendon of the semispinalis-spinalis; 2, lower tendon of the semispinalis-spinalis; and 3, semispinalis-spinalis muscle-tendon unit. The semispinalis-spinalis is visible at the bottom, the longissimus in the middle and the iliocostalis at the top of the picture. Anterior is to the right.
inspected the mean residuals by ecological group. Finally, we performed a phylogenetic MANCOVAs in mvMORPH (Clavel et al., 2015) and phylogenetic ANCOVAs using phytools (Revell, 2012) to test whether our results held when taking into account the phylogenetic relationships. To explore which groups differed from one another we ran PGLS regressions using the package Caper in R (Orme et al., 2013) and used the residuals (independent of variation in overall size) as input for ANOVAs. When significant differences between groups were detected, we inspected the means of the residuals per ecological group to assess the directionality of these differences.
All analyses were performed in R (R Core Team, 2020) and α was set at 0.05.

| RESULTS
Quantitative data on muscle segment lengths and muscle architecture for each individual are presented in Tables S1 and S2.

| Variation in muscle origin and insertion
In arboreal species, the lower and upper tendons as well as the contractile part of the semispinalis-spinalis and iliocostalis muscles generally span more vertebrae (Table S1,

| Variation in muscle architecture
The principal component analysis performed on the overall data set including all body parts (Table S2) extracted two axes jointly explaining over 95% of the variation in the data set. The first F I G U R E 2 Phylogenetic relationships between the species included in this study (based on Pyron et al., 2013). Ecological groups are indicated by colored circles next to the species' names.
component is a descriptor of overall size with larger species with heavy muscles having positive loadings (Table S3). The second axis differentiates between species that are overall longer versus more robust (Table S3).
The PCAs performed on the muscle architecture data for the three body parts separately show similar results (Figure 4).
Overall, the first two axes explain between 84.4% and 87.9% of the variation in the data set ( Table 2). The first axis differentiates species with heavy epaxial muscles with a greater cross-sectional area such as Agkistrodon piscivorus or P. regius that cluster to the right of the plot from more slender species that have a less developed (i.e., lower mass and cross-sectional area) musculature like P. guttatus or Leptophis ahaetulla to the left. The second axis is mostly determined by variation in muscle fiber length with species with longer fibers such as Natrix natrix clustering towards the positive side of the axis.

| Ecological differences
Traditional MANCOVAs showed a significant effect of snout-vent length on muscle architecture for all three body regions (Table S4).  (Tables S6 and S7).
The phylogenetic MANCOVAs showed a significant impact of body size on muscle architecture (Table 3)  dorsi were the only variables that differed between ecological groups (Table 4). Post-hoc analysis on the residual data derived from the PGLS regressions showed differences between all ecological groups in the mass of the semi-spinalis-spinalis and longissimus dorsi in the anterior part of the body. Moreover, the PCSA of the longissimus dorsi was different between aquatic species on the one hand and arboreal and terrestrial species on the other hand (Table S8). At mid body, aquatic snakes were again different from other ecological groups in the mass of the semispinalis-spinalis.
However, for the longissimus at mid-body aquatic snakes differed only from that in arboreal species (Table S8). Finally, at the posterior body, aquatic snakes differed again from the two other groups in both the mass and the PCSA of the longissimus dorsi. Overall, aquatic species had heavier muscles with a greater force generating capacity.

| Phylogenetic signal
Blomberg's multivariate K showed moderate but non-significant phylogenetic signal in the dataset for each of the body regions analyzed

| DISCUSSION
In our data set, the phylogenetic signal in the muscle architecture data was not significant. This is quite unexpected as axial muscle anatomy has been previously used as a systematic character in snakes (Mosauer, 1935). Mosauer (1935) recognized three distinct types characterizing boids, colubrids, and vipers with little variation within each group despite the fact that he had examined species with very different ecologies. However, later authors (Auffenberg, 1966;Gasc, 1981) have suggested this to be oversimplified and the studies of Ruben (1977) and Jayne (1982) clearly demonstrated variation in the semi-spinalis-spinalis system in relation to ecology.
Our results similarly suggest that habitat use, and more specifically locomotion in different habitats, is an important driver of muscle anatomy in snakes. These results echo findings for other taxa. For example, Omura (2019) demonstrated that the abdominal muscles were better developed in terrestrial salamanders as they function to retain the visceral mass against gravity. Conversely, in aquatic salamanders, the lateral hypaxial muscles that contribute to body flexion are more strongly developed (Omura, 2019;Omura et al., 2014Omura et al., , 2015. In our data set, aquatic snakes stand out from terrestrial and especially arboreal snakes in having generally shorter, heavier, and more forceful epaxial muscles, particularly the longissimus dorsi. These results are in line with previous results for snakes showing that the musculoskeletal system of aquatic snakes is quite different from terrestrial and arboreal ones. For example, the costocutaneus muscles that help flatten the body during swimming have been shown to have a more dorsal insertion in fully aquatic snakes (Voris & Jayne, 1976).
Moreover, seminal work by Jayne (1982) showed that arboreal snakes had longer tendons in the semispinalis-spinalis compared with aquatic snakes, similar to our observations. The elongation of the tendons of the epaxial muscles likely brings a mechanical advantage by increasing the muscle lever arms and by consequence the cantilever ability of the body (Jayne, 1986), which may provide an advantage during gap bridging (Jayne & Riley, 2007). Our data also show that the muscle belly of the iliocostalis muscle is also longer in arboreal snakes suggesting that arboreal species may be characterized by an overall elongation and reduction in weight of the axial muscle system. In contrast, aquatic species appear to have shorter tendons and muscles and more robust and heavier muscles.  (Auffenberg, 1966;Gasc et al., 1989;Jayne, 1982;Moon & Gans, 1998;Pregill & Pregill, 1977;Ruben, 1977). For example, the semispinalis and longissimus muscles are also important for striking in snakes as they are the two main extensors of the spine (Young, 2010). These muscles may therefore be involved in different functions and as such may be subject to trade-offs.
Moreover, it has been suggested that the semispinalis-spinalis and iliocostalis muscles may function to modulate the activity of the more posterior longissimus and as such may show different constraints (Moon & Gans, 1998). Previous studies have shown that muscle activity propagates from the front to the back of the body and that the longissimus dorsi is activated later than other muscles (Gasc et al., 1989;Jayne, 1988aJayne, , 1988bMoon & Gans, 1998). This suggests functional differences between the different axial muscle groups.
Interestingly, not all body parts showed the same signal with the posterior part showing differences only in the longissimus dorsi. Auffenberg (1958) discussed columnar variation and showed that the anterior muscle segments are generally shorter but did not discuss differences between species with different locomotor ecologies. Unfortunately, few authors have discussed the variation in muscle architecture along the vertebral column in detail (but see Nicodemo, 2012;Pregill & Pregill, 1977). The different signal observed in our data suggests different functional roles for different body parts and suggest that mainly the anterior and mid-body axial muscles may be involved in generating the power for locomotion.
The posterior body may be involved in grasping (arboreal species), controlling the amplitude of the rear of the body during undulations, or in other functions beyond locomotion (copulation, housing of the reproductive organs) that may constrain the development of the axial muscles.
Overall, our results show that aquatic species have shorter and heavier muscle segments than arboreal or terrestrial species. This makes intuitive sense as they have to move through a denser and more viscous environment (Vogel, 1994). Moreover, there is a fundamental difference between muscle activity during aquatic and terrestrial lateral undulation (Jayne, 1985(Jayne, , 1988a(Jayne, , 1988b. Whereas during both swimming and terrestrial locomotion all muscles (SSP, IC, and LD) are activated simultaneously, the wave of activation propagates faster down the body than the kinematic wave during swimming (Jayne, 1985). This suggests that the functional demands on the axial muscles may differ for species that spend a greater proportion of their time swimming.
Although our results suggest a strong relationship between axial muscle anatomy and the locomotor environment in snakes, it would be important to enlarge the data set. Both adding additional species as well as species with different ecologies including fossorial species would be important. Similarly, exploring whether a similar signal is present in other muscles or not would be important. Muscles like the multifidus, hypaxial muscles, or levatores costae are known to play an important role during locomotion and need to be explored.
Finally, our study considers only a few individuals per species and thus neglects potential intraspecific variability and sexual dimorphism (Bonnet et al., 1998;Penning, 2018

ACK N OWLED G M ENTS
We would like to thank E. Porquet for help with the dissections and two anonymous reviewers for helpful suggestions on an earlier version of the manuscript. This study was supported through funding provided by the Agence Nationale de la Recherche to the project 'Dragon2' (ANR-20-CE02-0010).

DATA AVA I L A B I L I T Y S TAT E M E N T
All data are available in the paper and in the supplementary information.