Evolution of Phenotype and Mitochondrial Genome Reveals Limbless and Body-elongated Squamates may Change Their Energy Basis for Locomotion

Zeng WANG,Wei WU,Jinlong REN,Changjun PENG,Dechun JIANG and Jiatang LI,3,4＊

1 CAS Key Laboratory of Mountain Ecological Restoration and Bioresource Utilization &Ecological Restoration and Biodiversity Conservation Key Laboratory of Sichuan Province,Chengdu Institute of Biology,Chinese Academy of Sciences,Chengdu 610041,Sichuan,China

2 University of Chinese Academy of Sciences,Beijing 100049,China

3 Southeast Asia Biodiversity Research Institute,Chinese Academy of Sciences,Yezin Nay Pyi Taw 05282,Myanmar

4 CAS Center for Excellence in Animal Evolution and Genetics,Chinese Academy of Sciences,Kunming 650223,Yunan,China

Abstract Limb reduction in Squamata present the dramatic characteristic to focus and usually accompanied with particularly morphological modifications,impacting tremendous locomotion changing and might generate different energy requirement.Herein,we combined both morphological and mitochondrial genomic data to explore the evolution of phenotypic transformation and mitochondrial genome of limbless and body-elongated squamates.We collected phenotypic measurements of 503 individuals,representing limbed or limbless taxa across all major lineages in Squamata to investigate the morphological correlations with limb-reduction.Furthermore,we provided the mitochondrial genome of the representative limbless and elongated species Dibamus bourreti (Angel,1935) to detect selective constraints on limbless clades with published mitogenomes of other squamate reptiles.Our results evidenced that body elongation had certain negative relationship with limbreduction in Squamata lineage and Lacertilia lineage (R=-0.495,P ＜ 2.2e-16;R= -0.332,P=1.1e-13,respectively),while tail length showed slight correlation in both clades (R=0.156,P=4.3e-04;R= 0.192,P=2.1e-05,respectively).Besides,detection demonstrated that ATP6 has experienced accelerated evolution among limbless lineages,suggesting selective pressure on mitogenomes may play an essential role in energy disparity for locomotion of limbed and limbless squamates.

Keywords limb-reduction,mitogenome,morphology,selective pressure,Squamata

1.Introduction

Remarkable morphological variations (i.e.the relative dimensions of the body,tail,and limbs),characteristic evolution and biotic interchange form the pattern of biodiversity and become one of the major goals in evolutionary biology (Carroll,1997;Futuyma,2005;Young

et al.

,2007;Ren

et al.

,2017;Jiang

et al.

,2019;Yu

et al.

,2021).Limb reduction in vertebrate exhibits the most intriguing trait,which repeatedly evolved in amphibians (Parra-Olea and Wake,2001;Urosevic

et al.

,2016),reptiles (Wiens

et al.

,2006;Miralles

et al.

,2015;Bergmann

et al.

,2020),and mammals (Bejder and Hall,2002;Law

et al.

,2019).Squamates,in particular,present an excellent model to study this dramatic body form transformation.Transitions from a fully pentadactyl form to an almost or completely limbless body form have independently occurred in squamate reptiles (Wiens

et al.

,2006;Brandley

et al.

,2008),and frequently accompanies with several morphological alterations like body elongation,digit loss or tail lengthening on particular taxonomic groups such as

Bachia

,

Chalcides

and

Lerista

(Gans,1975;Lande,1978;Caputo

et al.

,2000;Wiens and Slingluff,2001;Galis

et al.

,2010;Lee

et al.

,2013).However,the morphological correlations with limb-reduction remains unexplored in the whole squamates lineage.Reduction of limbs in squamates appears to induce locomotory mode changing (Gans,1962;Morinaga and Bergmann,2020).Locomotion of quadrupedal form is mainly supported by the rhythmic and coordinated movements of four limbs,whereas the propelling wave with different undulation pattern propagated along the trunk becomes the predominant mode in snake-like body forms (Gans,1986;Renous

et al.

,1998;Bergmann and Morinaga,2019;Bergmann

et al.

,2020),suggesting enormous locomotory difference might exist between limbed and limbless squamates.Furthermore,locomotion is energy-consuming,with most of the required energy supplied by mitochondria through oxidative phosphorylation (OXPHOS) (Boore,1999;Das,2006;Sun

et al.

,2011).Studies in birds and fishes indicated that locomotive patterns were closely related to the evolution of mitogenome (Shen

et al.

,2009;Sun

et al.

,2011).In this regard,we conjecture that phenotypic changes accompanied with locomotion might influence the evolution of mitochondrial genome in squamates.Dibamids are a group of elongated,fossorial,and limbless lizards (only males have flap-like traces near cloaca) living beneath leaf-litter,rocks,and rotting logs (Townsend

et al.

,2011;Quah

et al.

,2017),therefore,they arouse much interest for their typical phenotype and locomotory style.In the present study,we provided the complete mitochondrial genome and morphologic characteristics of

Dibamus bourreti

(Angle,1935) combining with other squamates to investigate the evolution of phenotypic and mitogenome in limbless and body-elongated squamate reptiles.

2.Material and Methods

2.1.Morphological data collection and analysis

An adult male of

D.bourreti

(voucher No.CIB118200) was collected from Mangshan,Hunan,China,muscle tissues were taken and preserved in 95% ethanol and then stored at -20°C in the Herpetological Museum,Chengdu Institute of Biology,Chinese Academy of Sciences,Chengdu,China.We measured snoutvent length (SVL) and tail length (TL) using a measuring tape to the nearest 1 mm.Hindlimb length (HL) (the flap structure) was obtained by a digital slide caliper to the nearest 0.01 mm.All measurements were following the description of Bergmann and Irschick (2010).Morphometric data of other squamates was collected from previous studies (Wiens

et al.

,2006;Siler

et al.

,2011;Morinaga and Bergmann,2017;Bergmann and Morinaga,2019).Totally,we gathered 503 species including limbed and limbless species (Table S1),covering all major lineages of squamates.To summarize the pattern of morphological variation,we divided our samples into four groups based on the following： limbed lizards,limb-reduced lizards,limbless lizards,and snakes.Limb-reduced lizards referred to those lacking all fingers or toes but remaining forelimb or hindlimb,while limbless lizards lacked all limbs without any vestige.A principal component analysis (PCA) was conducted using the “prcomp” function and

ggbiplot

package in R (R core team,2020 )based on the dataset.Then we tested the characteristics correlation in Squamata clade using

Hmisc

package in R (R core team,2020).The traits were defined as 0 (limbless),1 (limb-reduced),or 2 (limbed).Snakes were included in limbless group therein.Considering the morphological peculiarity of snakes,we did the correlation test in Lacertilia clade independently.

2.2.Sequencing,mitogenome extraction and annotation

Total genomic DNA was extract from muscle tissue and a pairend Illumina library was constructed with insert size ranging from 300 bp to 500 bp.The library was sequenced on Illumina Hiseq 2500 to generate raw data.Clean data was obtained after removing contaminated reads,low quality reads and reads with more than 5% ‘N’ bases of raw data and the complete mitochondrial genome of

D.bourreti

was extracted using Novoplasty (Dierckxsens

et al.

,2017) from it.The mitochondrial genome annotation web server (MITOS) was used for annotation (Bernt

et al.

,2013) and the tRNA genes were scanned by tRNAScan-SE2.0 online website (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Chan,2016).

2.3.Detecting selective pressure

Complete mitochondrial genomes were collected form GenBank for 27 species,representing disparate clades of Squamata (Table S2).The species tree was downloaded from Time Tree Database (http://www.timetree.org/).Based on the topology,the branch model of CODEML in PAML4.7 (Yang,1997) was used to evaluate the

ω

parameter,which means the ratio of nonsynonymous substitution (dN) and synonymous substitution (dS) in protein-coding genes (PCGs).We tested the null hypothesis using the null model (model=0) firstly,meaning the

ω

ratio in all branches are equal.An alternative hypothesis was subsequently detected (model=2),allowing different

ω

ratios on branches of interest.All limbless clades were regarded as the same foreground and the others were background.Finally,we calculated the significance of the null hypothesis and the alternative hypothesis using Likelihood Rate Test (LRTs).Genes with a

P

value less than 0.05 were considered as evolving with significantly faster rate in foreground branches.

3.Results

3.1.Analyses of character correlation

All species were clustered into four distinct groups in PCA analysis:limbed lizards,limb-reduced lizards,limbless lizards,and snakes (Figure 1) based on the four indices (SVL,TL,FL,and HL).The first component from the PCA result showed that PC1 explained 56.6% of variation (Table 1).Individuals loaded highly positively on PC1 represented the snake-like species,characterized by reduced limb and elongated body.By contrast,negatively loaded species were lizard-like with robust limbs.Therefore,we interpreted PC1 as the indicator for lizard-like or snake-like groups since it contrasted variables associated with body elongation (SVL with positive loadings).PC2 explained 32.6% of the variance and characterized by strongly loadings of SVL and TL.PC3 and PC4 explained low variation (Table 1) and therefore were not considered further.The correlation test between morphological traits indicated certain negative correlation between SVL and limb-reduction in both Squamata lineage and Lacertilia lineage (

R

= -0.495,

P

＜ 2.2e-16;

R

= -0.332,

P

=1.1e-13,respectively),while slight relationship existed between TL and limb-reduction in either Squamata or Lacertilia clade (

R

=0.156,

P

=4.3e-04;

R

=0.192,

P

=2.1e-05,respectively) (Table 2).

Figure 1 PCA plot of 503 squamate species based on the first two principal components (PCs).Ellipses indicate different groups with 68% confidence intervals.

Table 1 Morphometric PCA result for all variables.PCA loadings and percent explained variance for each PC are presented. n = 503.

Table 2 Correlation test of traits.

3.2.Mitogenome characterization of

We totally obtained 241.5 million reads (～36 Gbp) clean data of

D.bourreti

from Illumina Hiseq 2500 platform.The circular mitogenome of was 16 706 bp in size (Table 3,Figure 2),with base composition of 31.6% A,23.9% T,30.0% C and 14.5% G.The mitochondrial genome contained 13 protein-coding genes (

ND1-6

,

ND4L

,

COI-III

,Cyt

b

,

ATP6

,and

ATP8

),two rRNA fragments (12s and 16s rRNA),22 tRNA fragments,and a control region (D-loop) (Table 1).Most genes were coded on the H-strand except

ND6

and 8 tRNAs (tRNA-Gln,Ala,Asn,Cys,Tyr,Ser,Glu,and Pro).All the PCGs were initiated with typical start codon ATG (Chen

et al.

,2019),except for

COI

with GTG.For the PSGs,the longest gene was

ND5

(1824 bp) and the shortest was

ATP8

(168 bp).The 22 tRNAs ranged from 64 bp of tRNA-Cys to 75 bp of tRNA-Leu in length.The 12S rRNA and 16S rRNA were 948 bp and 1526 bp respectively.The mitogenome of

D.bourreti

had been deposited in GenBank under accession number MW368917.

Figure 2 Mitochondrial gene organization of Dibamus bourreti.Genes encoded by the H-strand are outside the circle,while genes inside the circle are coded by the L-strand.All tRNA genes were denoted by IUPACIUB one-letter amino acid abbreviation.

3.3.Selection on mitogenomes of limbless squamates

The topology of the phylogenetic tree reveals that limbless species were not recovered as a monophyletic group,suggesting limbreduction evolved independently several times (Figure 3).In our result,the ratio of dN/dS substitutions of all genes were lower than 1,demonstrating more nonsynonymous substitutions than synonymous substitutions (Table 4).According to LRT,the selective pressure of

ATP6

therein has significantly difference between foreground and background branches (

ω

=0.0376,0.0191,

P

＜ 0.05,respectively),suggesting that

ATP6

of limbless species (foreground branch) experienced rapid evolution compared to the background branch.

4.Discussion

Here we integrate the morphological data to elucidate the pattern of morphological variations and test the morphological correlations with limb-reduction in squamates.All species were clustered into four groups based on the first principal component in our PCA plot (Figure 1),and correlation analysis showed definite relationship between SVL and limb-reduction in both Squamata and Lacertilia lineages (Table 2).It was obvious that quadrupedal squamates varied from snake-like squamates,but the discrepancy among limbless lizards andsnakes should also be noted.Distinctive modes for elongation might explain the different clusters of limbed,limb-reduced,limbless lizards and snakes in squamates.In squamates,limbless species elongated body forms in three ways:lengthening the trunk,the tail or the both (Wiens

et al.

,2006;Brandley

et al.

,2008;Morinaga and Bergmann,2017).Snakes tended to display elongation in both trunk and tail region,whereas lizards appeared to elongate tails multiple times than trunks (Woltering,2012).However,correlation between TL and limbreduction showed slight support in either Squamata or Lacertilia in our study (Table 2),suggesting the lengthening of tails associating with limb reduction might occur in certain clades such as Anguidae (Wiens and Slingluff,2001).Furthermore,our results demonstrated the association between body elongation and limb reduction in Squamata,consistent with previous studies in amphibians (Parra-Olea and Wake,2001;Bonett and Blair,2017),reptiles (Wiens

et al.

,2006;Brandley

et al.

,2008;Grizante

et al.

,2012) and mammals (Buchholtz and Schur,2004;Buchholtz

et al.

,2007;Law

et al.

,2019),prompting that limbreduction accompanied with elongated body forms tend to be a pervasive rule in vertebrate.

Figure 3 Phylogeny of 28 squamate species included in this study.Lines marked with red represent the elongated,limbless species and black lines represent the limbed,lizard-like species， respectively.

Table 3 Annotation of whole mitogenome of Dibamus bourreti.

Table 4 The selective pressure of 13 protein-coding genes in mitogenomes.

Previous studies had put forward three locomotive patterns for squamates with different degrees of limb loss:limbed squamates used their limbs for movement,limb-reduced squamates moved with their remaining limbs alternately with body swing and limbless squamates (including limbless lizards and snakes) relied on completely body swing (Gans,1986;Gans and Gasc,1990;Gans and Fusari,1994).Distinct locomotive modes were considered to associate with energy metabolic as well as the evolutionary pattern of mitochondrial proteins (Shen

et al.

,2009;Sun

et al.

,2011;Jacobsen

et al.

,2015).Mitochondria usually provide up to 95% energy in metabolism,which was directly produced by ATP synthase of mitochondrial respiratory chain,and support for locomotion (Sperl

et al.

,2006;Shen

et al.

,2010;Mitterboeck

et al.

,2017).Researches in birds and insects prompted that selection pressure on mitochondrial genomes was closely related to locomotive abilities (Shen

et al.

,2009;Mitterboeck and Adamowicz,2013).Sun

et al.

(2011) evidenced that fishes with different locomotive patterns faced distinctive energy requirement and selective pressure on mitochondrial proteins vary among these groups,suggesting selective constrains exhibited by mitochondrial genome was due to the restriction of metabolism.By evaluating selective constrains on mitochondrial protein-coding sequences,our work showed that

ATP6

had undergone significantly accelerated evolution among limbless lineages comparing to the limbed lineages (Table 4).

ATP6

is a subunit of the F-proton channel in ATP synthase,catalyzing synthesis of ATP from ADP (Saraste,1999).Jacobsen

et al.

(2015) investigated mutations in

ATP6

of two migratory groups in eels and demonstrated that positive selection in

ATP6

involved with different locomotive patterns in two lifehistory.Likewise,studies in mammals showed that mutations within

ATP6

indicated the association with differences in metabolism and selection linked to energetics (Fontanillas

et al.

,2005;Kucharczyk

et al.

,2010).Limbless and elongated forms facilitate undulate locomotion (Webb,1982;Gans and Fusari,1994;Morinaga and Bergmann,2019) and moving through narrow tubes or cluttered habitats,experiencing less drag than their limbed counterpart (Gans,1962;Gans and Gasc,1990;Mehta

et al.

,2010;Morinaga and Bergmann,2020).Therefore,we speculate that limbless squamates perhaps receive different energetic requirement for their undulate locomotory pattern,and selective constraints on mitochondrial genes might play an important role in the evolution of locomotory patterns in squamate reptiles.But the role of

ATP6

involved in energy metabolism in different groups and whether other genes in the whole genome change to meet energetic requirement for locomotion need more evidence.

Acknowledgement

We thank Mr.Hang YANG for the fieldwork.We also thank Mr.Zhongliang PENG and Ms.Dan WANG for their help in data processing.This work was supported by the National Natural Science Foundation of China (31772434;32070410);the Second Tibetan Plateau Scientific Expedition and Research Program(STEP) (2019QZKK0501);the International Partnership Program of Chinese Academy of Sciences (151751KYSB20190024);Biological Resources Programme,Chinese Academy of Sciences (KFJBRP-017-14);and the Distinguished young scholars in Sichuan Province (2021JDJQ0002).