In our study, 12 body size related gene families and 26 body-size-related genes (a total of 101 genes, Supplementary Table 1) were selected from body size regulation signaling pathways and pertinent literatures. The 28 representative species (Supplementary Table 2) includes 12 lizards: Anolis carolinensis, Eublepharis macularius, Gekko japonicus, Hemiphyllodactylus yunnanensis, Lacerta agilis, Paroedura picta, Podarcis muralis, Pogona vitticeps, Salvator merianae, Shinisaurus crocodilurus, Varanus komodoensis, Zootoca vivipara, and 16 snakes: Crotalus tigris, Hydrophis melanocephalus, Laticauda laticaudata, Naja naja, Notechis scutatus, Ophiophagus hannah, Pantherophis guttatus, Pantherophis obsoletus, Protobothrops flavoviridis, Protobothrops mucrosquamatus, Pseudonaja textilis, Ptyas mucosa, Python bivittatus, Thamnophis elegans, Thamnophis sirtalis, Thermophis baileyi. The third-generation genomes of G. japonicus and H. yunnanensis were obtained from our laboratory (data has not been published yet), and the remained genomes were downloaded from NCBI database¹, Ensembl database², and GigaDB database³. The 101 body-size-related amino acid sequences of Homo sapiens, Xenopus tropicalis, Gallus gallus, A. carolinensis and G. japonicus were acquired from NCBI and used as queries to identify the coding sequences of body-size-related genes. TBLASTN searches were conducted against the genomes of the studied species, and the E-value was set to 1e-5. For genomes without annotation files, we predicted the candidate genes by GeMoMa v1.7.1 (Keilwagen et al., 2019) using known reference genomes with annotation. Incomplete sequences and sequences containing premature stop codons were abandoned. All sequences were removed stop codons and were aligned using PRANK v170427 (Löytynoja, 2014) at the codon level. Poor aligned and highly variable regions were removed using Gblocks v0.91b (Castresana, 2000) with relatively strict parameters (“-t = c, -b5 = h”). Information about orthologous genes of each body-size-related gene and accession numbers of genomes were provided in Supplementary Tables 1, 2. The body size attributes (mean, mode, median, etc.) of squamate suborders (amphisbaenians, lizards and snakes) are totally different (Feldman et al., 2016), so we separately performed evolutionary analyses on lizard-only and snake-only datasets. High-quality multiple sequence alignment (MSA) files were from lizard-only dataset containing orthologous genes of 12 lizard-species, and snake-only dataset containing orthologous genes of 18 snake-species.

In evolutionary biology, comparisons of non-synonymous/synonymous substitution ratios (dN/dS) are used to quantify the influence of natural selection (Kimura, 1983). Different values represent different types of selection: ω < 1, ω = 1 and ω > 1 indicate purifying selection, neutral selection and positive selection, respectively. CODEML program implemented in PAML v4.9 (Yang, 2007), was applied to detect the selective pressure by site model, branch model, branch-site model and clade model. The sequence length was an exact multiple of three, and the tree topology of lizards and snakes used as input for the CODEML were pruned from the consensus reptilian phylogenetic tree (Tonini et al., 2016; Figure 1).

Firstly, we used three evolutionary models (M0, M8a, M8) to estimate site-specific selection. Model 0 which assumed that one ratio occurred across the phylogeny was used to calculate the global ω values representing the average evolutionary rate of genes, and two pairs of site models (M8a vs. M8) were applied to identify the sites under positive selection in each body-size-related gene. The likelihood ratio test (LRT) which approximately follows a chi-square distribution was employed to compare hierarchically nested models (Bielawski and Yang, 2003). p < 0.05 represents that the two models are statistically significant. Positively selected sites were detected using Bayes empirical Bayes (BEB) analysis with posterior probabilities (PP) of >0.8 (Yang, 2005). Fisher’s Extract Test was used to determine whether there was a significant association between the number of selected genes in lizards and that in snakes or not.

Then, in order to test branch-specific evolutionary rates, we performed branch models with a two-ratio model which allowed ω values to vary between the foreground and background branches. The two-ratio model is compared to a one-ratio model that enforces the same ω values for all branches. According to the body size of squamates from Feldman et al. (2016) (Supplementary Table 2), we applied K-means cluster analyses to separately classify lizards and snakes into two categories (large-bodied and small-bodied). S. merianae, V. komodoensis in lizards, and P. bivittatus, O. hannah in snakes were relatively large-bodied species. In the lizard-only and snake-only datasets, the terminal branches of large-bodied species were set as the unified foreground branches, and the remaining lineages were background branches. Large-bodied species were compared against small-bodied species, providing a clear understanding of body-size-related genes subject to different selective pressure in two categories. Spearman’s correlations between the estimated ω values for branches (large lizards, small lizards, large snakes, and small snakes) were implemented in IBM SPSS Statistics for Windows, Version 25.0 (IBM Corp., Armonk, NY, United States).

To investigate whether episodic positive selection was limited to a subset of sites along particular lineages, branch-site models were conducted which allowed variation in ω values across branches and sites in the phylogeny. Besides the unified branches of large-bodied species, some particular lineages of interest were also set as the foreground branches (Supplementary Table 3). BEB analyses were also used to count the posterior probability (PP > 0.8). For each body-size-related MSA, the alternative model Ma was compared with the null model Ma0 which ω value was fixed to 1.

To further test for the selective pressure acting on large-bodied species relative to other small-bodied species, Clade model C (CmC) were utilized in the lizard-only and snake-only datasets. CmC assuming that a class of sites evolve conservatively along the phylogeny, and others evolve contrastly between foreground and background branches, was compared against the null model M2a_rel (nearly neutral).

In view of the large number of candidate genes and relatively few species, all results of branch model, branch-site model and clade model were corrected for multiple testing bias using Benjamini–Hochberg method implemented in q-value package (R v4.0.3) (R Core Team, 2020). Q-value indicates the false discovery rate (FDR) using the critical value 0.05, which is calculated on the basis of p value.

To explore the potential association between the evolutionary rate of body-size-related genes and phenotypes (body length, body mass), phylogenetic generalized least squares (PGLS) regression were implemented in the R package ‘caper’ (Orme et al., 2018). The time-calibrated phylogenies were acquired from the consensus tree of reptile (Tonini et al., 2016). The phylogenetic signal (lambda, λ) was estimated by the maximum likelihood (ML) method, which can vary from 0 to 1. λ = 1 indicates that observed data obey Brownian motion model of evolution and have a strong phylogenetic signal; λ = 0 means that the phylogeny occurs independently.

As the method suggested by Montgomery et al. (2011), the evolutionary rate (root-to-tip ω) were calculated by free-ratio model implemented in the CODEML program of PAML v4.9. The root-to-tip ω of each species is the average of ω values from the last common ancestor (LAC) to each terminal branch, which removes temporal effects on ω values and reflects the average evolutionary rate of genes. If dN or dS is less than 0.0002 leading to ω values extremely big or small, it is labeled as an outlier ‘N/A.’ The phenotypic data of lizards and snakes studied were obtained from the publication (Zhang et al., 2005; Tonini et al., 2016; Xing et al., 2018). Moreover, in order to improve normality, all phenotypic data and root-to-tip ω values were log₁₀-transformed for subsequent analyses.

To address the impact of ecological variables (microhabitat, climatic niche, range size, rate of range expansion) on the evolutionary rate of body-size-related genes, we first performed univariate PGLS regression analyses separately on the lizard-only and snakes-only datasets (Pagel, 1999). Then, variables significantly associated with body-size-related genes (p < 0.05) were selected to conduct a series of phylogenetic multiple regression analyses, considering all possible combinations of these variables to build candidate models. The models with lowest Akaike Information Criterion (AIC) values were considered as the best-fitting model overall (Burnham and Anderson, 2002). Finally, for the best-fitting model, we used R code from Moen and Wiens (2017) to estimate standardized partial regression coefficients (SPRC). SPRC explain the contribution of each ecological variable to the best-fitting model, and show the impact of variables on evolutionary rates of body-size-related genes when others are held constant (Sokal and Rohlf, 1995).

We obtained the microhabitat data of 28 species from the IUCN database (IUCN, 2021) and previously published literatures (Bars-Closel et al., 2017; Meiri, 2018; Supplementary Table 4). Following the classification method of Bars-Closel et al. (2017), each species was assigned to one microhabitat type where they were active, including terrestrial, arboreal, fossorial, aquatic, semi-arboreal, semi-aquatic, and semi-fossorial. Pie et al. (2017) estimated phylogenetic principal components (PCs) for 19 climatic variables for squamate species. We chose PC1, PC2, PC3, ln-transformed Bio1 (annual mean temperature) and ln-transformed Bio12 (annual precipitation) as independent variables, and then extracted climatic niche data of 17 species studied (Pie et al., 2017; Supplementary Table 5). The number of ecoregions was defined as range size in the study. As Pyron and Wiens (2013) suggested, we first retrieved the geographic range of 28 species studied from the IUCN database, and assigned each species to one or more of 12 ecoregions (Li and Wiens, 2022; Supplementary Table 6). To improve normality, we transformed range size based on the equation log₁₀(x + 0.5), where x is the number of ecoregions and adding 0.5 avoids assigning zero to species found in only one ecoregion (Li and Wiens, 2022). As a measure of the rate of range expansion for each species, the total number of ecoregions was divided by the species’s stem age (ecoregions/age) (Li and Wiens, 2022). Thus, species that occur in many regions despite being relatively young will have a relatively rapid rate of inferred range expansion; conversely, species that occur in few regions despite being older will have a lower rate of range expansion.

Finally, we tested the relationship between ecological variables (microhabitat, climatic niche, range size, rate of range expansion) and phenotypes (body length, body mass) in 28 species of lizards and snakes. The PGLS regression analyses were performed in lizard-only, snake-only and squamate datasets. Specially, in the squamate dataset, we only examined the effect of ecological variables on body mass, because the measurement methods of body length in lizards and snakes were different. The specific analytical methods were the same as above.

The model 0 showed that the ω values of these 101 body-size-related genes were 0.121 ± 0.096 (mean ± SD) in the lizard-only dataset and 0.175 ± 0.144 (mean ± SD) in the snake-only dataset (Supplementary Table 7). The global ω values of snake were higher than those of lizard, indicating that body-size-related genes in snakes evolved faster than in lizards. In addition, the global ω values of lizards and snakes were much less than 1, suggesting that these genes were under strong purifying selection constraining their critical function in the body size evolution of squamates. The global ω values distribution of lizards and snakes were different, among which, the highest and lowest ω values (0.426 and 0.00225) in lizards were observed for NSD1 and GDF11, and the highest and lowest ω values (0.759 and 0.00491) in snakes were found for NPC2 and SOX4.

In site models, of 101 body-size-related genes in lizards, 9 (∼8.9%) supported the Model 8 (positively selected) (Supplementary Table 8), whereas in snakes, 21 (∼20.8%) were under positive selection (Supplementary Table 9). Fisher’s Extract Test demonstrated that there was a significant difference in the number of genes subject to positive selection between lizards and snakes (p = 0.028, two-tailed), which indicated that these clades might have disparate evolutionary patterns of body size. Notably, COL10A, GHR and NPC1 were co-positively selected in the lizard-only and snake-only datasets.

Furthermore, it was also observed that LCORL (0.348 and 0.409) and INSL5 (0.482 and 0.491) in the lizard-only and snake-only datasets had relatively higher global ω values (Supplementary Table 7). However, considering only the number of sites with ω values greater than 1 and posterior probabilities (PP) greater than 0.8, OBSL1 (11) of lizards and NPC1 (18), NPC2 (15) of snakes had the highest number of positively selected sites.

For branch model analyses, free-ratio models fitted the data better than one-ratio models at 53 genes in lizards and at 24 genes in snakes, indicating the heterogeneous selection among the phylogeny. Then, we marked the datasets based on body size, with the large-bodied species (S. merianae and V. komodoensis in lizards, P. bivittatus and O. hannah in snakes) as united foreground branches. After FDR correction, 10 body-size-related genes in lizards and 9 genes in snakes favored two-ratio models better than one-ratio models (p < 0.05, q < 0.05). Contrasting the two datasets, GALNS was the only gene with significantly different evolutionary rates between large-bodied and small-bodied lineages in both lizard-only and snake-only datasets. The positively selected genes identified by branch model were shown in Table 1.

Spearman correlation analyses showed that the correlation between lizards and snakes had the lowest rho (ρ) for large lizards and small snakes (Spearman’s ρ = 0.653; p < 0.01). The ω values were highly similar across taxonomic clades: large and small lizards (Spearman’s ρ = 0.857; p < 0.01), and large and small snakes (Spearman’s ρ = 0.700; p < 0.01) (Table 2). These results demonstrated that through relatively less sites influenced, they were sufficient to be detected under positive selection, yet were no able to cause significantly different evolutionary rates of body-size-related genes between large and small species in each clade. Specially, the correlation of ω values for large lizards and large snakes was also relatively high (Spearman’s ρ = 0.730; p < 0.01), suggesting that lizards and snakes were likely to evolve toward larger body size with similar patterns (Table 2).

In branch-site models, the terminal branches of large-bodied species were also set as united foreground branches. After FDR correction, CDKN2C, IGF1, IRS2 in lizards and FBN1, NSD1 in snakes were detected positive selection signatures in large-bodied lineages, and positively selected sites with posterior probabilities > 0.8 were shown in Table 3. It was thus clear that lizards and snakes did not have the same genes under positive selection after FDR correction. However, before correction, FBN1 was identified as a positively selected gene in both lizards and snakes, and had two identical sites under positive selection with PP > 0.8. CDKN2C in lizards and NSD1 in snakes had the highest number of positively selected sites, 24 and 56, respectively, indicating that they might be the key genes for regulating body differentiation in squamates.

To further identify body-size-related genes under disruptive selection in the different body size lineages, we analyzed the lizard-only and snake-only datasets using Clade model C. We found significant evidence for shifts in selection pressure in large-bodied species relative to other small species in 62 of 101 body-size-related genes from the lizard-only dataset (p < 0.05, q < 0.05). Among them, the foreground branches of 14 genes (CCNB3, CDKN2C, COL10A1, FOXO1, GPC3, GPC4, GRB10, IGF1, INSL5, LEPR, mTOR, NSD3, OBSL1, SOX21) had ω values higher than background and greater than 1, suggesting that these genes experienced the relaxation of selective pressure and reached the level of positive selection (Supplementary Table 10). In the snake-only dataset, CmC detected that the partitions of 35 body-size-related genes (from a total of 95 genes, and 6 genes were not identified in large-bodied snakes) were significantly better suited relative to M2a_rel models, indicating divergent selection applying to different body size lineages. The ω values of the foreground branches in 8 of 35 genes were higher than those of background, which were strong evidence for relaxation of constraint. Except CDKN2C, the retained seven genes (AKT1, BMP15, CDKN2C, IGFBP6, IRS4, NSD1, SOX11, SOX5) with ω > 1 were affected by positive selection (Supplementary Table 10). Interestingly, CDKN2C was under divergent selection in both lizards and snakes, but only in lizards reached the level of positive selection.

In addition to large-bodied species, there are also some small-bodied species in lizards and snakes. The terminal branches of these large-bodied and small-bodied species were also set as the foreground branches, and we used branch-site model to explore similar or different genes and sites under positive selection in lineages with different body size. Evidence of positive selection was determined along the terminal branch of H. yunnanensis, the smallest species in lizards, at seven genes (ACAN, FOXO6, GDF6, GDF7, GHR, IGFBP2, IRS2), and along the terminal branch of T. baileyi, the smallest species in snakes, at four genes (FBN1, IGF1, PRKG2, SOX21) (Supplementary Table 11). Moreover, positive selection was detected along the terminal branch of V. komodoensis (the largest species in lizards) at six genes (BMP10, FOXO4, IGF1, INS, LEPR, NSD1), as well as the lineage leading to S. merianae (another large-bodied species in lizards) at four genes (BMP15, CDKN2C, IRS2, LEPR) (Supplementary Table 11). We also found positive selection signals along the terminal branch of P. bivittatus (the largest species in snakes) at two genes (CDKN2C, FBN1) and along the terminal branch of O. hannah (another large-bodied species in snakes) at five genes (BMP15, IGFBP6, IGFBP7, NSD1, SOX11) (Supplementary Table 11). Comparing the results, we discovered some positively selected genes (e.g., FBN1, IGF1, IRS2) associated with large-bodied and small-bodied species, LEPR only relevant to large-bodied species in lizards, and some genes (e.g., BMP15, CDKN2C, NSD1) only related to large-bodied species in both lizards and snakes.

Comparing body size of 28 species studied, E. macularius is larger than the species of Gekkonidae and S. merianae is larger than the species of Lacertidae. Similarly, C. tigris is larger than P. mucrosquamatus and P. flavoviridis, and T. baileyi is smaller than T. sirtalis and T. elegans, and P. mucosa is larger than P. guttatus and P. obsoletus. Therefore, we set some ancestral branches as foreground and performed branch-site models in the lizard-only and snake-only datasets (Supplementary Table 3). Five genes (GDF6, GPC4, IGFBP2, SH2B3, SOX11) at the ancestral branch of Gekkonidae and six genes (CCNB3, CNP, IGFBP3, SOX9, IRS2, SMO) at the ancestral branch of Lacertidae displayed significant evidence of positive episodic selection (Supplementary Table 12). It is easy to see that both IGFBP2 and IGFBP3 belong to the IGFBP family, which may participate in the regulation of body size reduction. Furthermore, branch-site models showed that LCORL (296 G 0.604, 575 L 0.767, 964 V 0.595) was positively selected along the lineage to the last common ancestor (LCA) of P. mucrosquamatus and P. flavoviridis, but had no sites with PP > 0.8 (Supplementary Table 12). Positive selection was also identified along the LCA of T. sirtalis and T. elegans at three genes (NPC1, 553 Y 0.650; ACAN, 965 D 0.808; LCORL, 265 C 0.564, 1070 V 0.665) and along the LCA of P. guttatus and P. obsoletus at FBN1 (727 A 0.785) (Supplementary Table 12). Particularly, LCORL showed the evidence of episodic positive selection along two ancestral branches leading to relatively large-bodied species, suggesting its potential function to evolve larger. Nevertheless, we did not find common sites under positive selection, indicating that different sites in one gene could influence the same function.

To explicitly explore the association between evolutionary rate of body-size-related genes and morphological variables (body length, body mass), PGLS regressions were conducted in the lizard-only and snake-only datasets. The results in lizards revealed the evidence of negative association between log (root-to-tip ω) and log (body length)/log (body mass) at four genes (AKT1, BMP1, IGF2, SOX7), and positive association between log (root-to-tip ω) and log (body length)/log (body mass) at the gene SOX5 (body length: R² = 0.319, p = 0.033; body mass: R² = 0.269, p = 0.048) (Table 4). In the snake-only dataset, we detected significant negative association between log (root-to-tip ω) and log (body length) at eight genes (BMP5, BMP7, FGF4, IGFBP5, IGFBP6, SH2B3, SOX17, SOX30) and between log (root-to-tip ω) and log (body mass) at seven genes (BMP5, BMP7, CDKN2D, GPC1, SH2B1, SH2B3, SOX17) (Table 4). In addition, positive association also were found between log (root-to-tip ω) and log (body length)/log (body mass) at the gene GPC6 (body length: R² = 0.360, p = 0.018; body mass: R² = 0.330, p = 0.023). Among them, the evolutionary rate of BMP5, BMP7, GPC6, SH2B3, SOX17, and GPC6 were significantly correlated with body length and body mass. In short, genes associated with phenotypic data were different in both lizards and snakes, but involved in two gene families (BMP and SOX family).

The growth and development of squamates living in different habitats are affected by many ecological factors, and genes usually control the cell proliferation, apoptosis and other processes of organisms. Hence, we performed phylogenetically controlled regressions to probe into the link between evolutionary rate of body-size-related genes and ecological variables (microhabitat, climatic niche, range size, rate of range expansion). We discovered significant correlation between log (root-to-tip ω) and logit (microhabitat) at seven genes (BMP4, BMP6, CDKN1A, GALNS, IGFBP3, IGFBP6, SOX8) in lizards, and at five genes (CCNB3, FGF4, SH2B1, SOX6, SOX7) in snakes (Supplementary Table 13). It was obvious that no common microhabitat-related genes were found.

Moreover, AKT1, BMP1, IGF2, SOX5 in lizards and ADIPOR2, FOXO1, GPC6, KCNJ2, SOX30 in snakes were strongly associated with range size (number of ecoregions) (Supplementary Table 14). We also found AKT1, BMP6, BMP7, FOXO1, GALNS, GDF9, IGFBP2, IGFBP6, IRS1, NSD2, PLAGL2 in lizards and ADIPOR2, KCNJ2, SOX30 in snakes that showed the evidence of association between log (root-to-tip ω) and rate of range expansion in species (ecoregions/age) (Supplementary Table 14). It was clearly stated that AKT1 in lizards, and ADIPOR2, KCNJ2 and SOX30 in snakes were correlated with range size and rate of range expansion (p < 0.05).

Furthermore, we detected that CDKN1B, NPC1 and SOX11 were significantly associated with climatic-niche rate for PC1; AKT1, BMP15, EIF4EBP1, GDF9, KCNJ2, SOX18, SOX5, and SOX6 were significantly correlated with climatic-niche rate for PC2; CCNB1, FOXO6, GDF7, LEPR, SOX11, SOX6, and SOX8 were significantly linked to climatic-niche rate for PC3; AIP, BMP15, CCNB1, FOXO1, NOG, and SOX14 were significantly related with annual mean temperature (Bio1); CDKN1B, NPC1, PLAGL2 and SOX11 were significantly concerned with annual precipitation (Bio12) in lizards (Supplementary Table 15). In the snake-only dataset, the evidence of significant association were discovered between log (root-to-tip ω) and PC1 at four genes (CCNB1, CDKN1A, FBN1, SOX30), and PC2 at nine genes (ADAMTS10, AKT1, FBN1, FOXO6, GPC6, IGF1, IGFBP6, KCNJ2, SOX30), and PC3 at eight genes (ADAMTS10, AIP, BMP2, CDKN1A, GPC1, LEPR, OBSL1, SOX7), and ln (Bio1) at nine genes (AKT1, FBN1, FOXO1, GPC1, IGFBP4, LEPR, SH2B1, SOX30, SOX7), and ln (Bio12) at six genes (ADIPOR2, BMP6, CCNB1, FBN1, LCORL, SOX30) (Supplementary Table 15).

Finally, combining the results of the above eight ecological variables on the evolutionary rate of body-size-related genes, it was found that AKT1, BMP15, BMP6, CCNB1, CDKN1B, FOXO1, GALNS, GDF9, IGFBP6, NPC1, PLAGL2, SOX5, SOX6, SOX8, SOX11 in lizards and ADAMTS10, ADIPOR2, AKT1, CCNB1, CDKN1A, FBN1, FOXO1, GPC1, GPC6, KCNJ2, LEPR, SH2B3, SOX30, SOX7 in snakes were significantly correlated with multiple variables (p < 0.05). Thus, based on these genes and variables, phylogenetic multiple regression analyses were performed. In the lizard-only dataset (Table 5), NPC1 and SOX11 were significantly associated with PC1 and Bio12, where PC1 explained 64–71%, whereas Bio12 accounted for 17–27% of the variation in the evolutionary rate; AKT1 and GDF9 were strongly concerned with PC2 and ecoregions/age, which separately 41–64% and 24–48% of the variation; GALNS and IGFBP6 were significantly linked to microhabitat and ecoregions/age, which respectively accounted for 23–24% and 19–30% of the variation. In subsequent analyses of snakes (Table 6), ADAMTS10 was significantly connected with PC2 and PC3, which separately explained ∼28% and ∼19% of the variation in the evolutionary rate; ADIPOR2 was strongly associated with number of ecoregions and ecoregions/age, which respectively accounted for ∼23% and ∼11% of the variation; CCNB1 was correlated with PC1 (∼18%) and Bio12 (∼28%); FBN1 was associated with PC2 (∼23%) and Bio12 (∼40%); FOXO1 was strongly linked to number of ecoregions (∼0%) and Bio1 (∼43%); KCNJ2 was correlated with PC2 (∼36%) and number of ecoregions (∼26%); SOX7 was correlated with microhabitat (∼10%) and Bio1 (∼38%); SOX30 was significantly concerned with number of ecoregions (∼48%) and Bio1 (∼28%). Notably, we discovered the evidence of strong link between the evolutionary rate of body-size-related genes and five ecological variables (PC1, PC2, PC3, Bio1, Bio12), and relatively weak correlation between the evolutionary rate and microhabitat.

Due to natural selection, body size of squamates might be potentially related to ecological factors. Eventually, we explored the effects of ecological variables on body size. As were shown in Table 7, in the lizard-only dataset, only number of ecoregions had a significant impact on body length and body mass (p < 0.05); whereas in the snake-only dataset, no ecological factors were significantly correlated with body size. Specially, we also found that there was a weak association between microhabitat and body mass in the squamate dataset (R² = 0.104, p = 0.056) (Table 7 and Supplementary Table 16).

Site model found the number of positively selected genes in snakes was higher than that in lizards. According to data on body size for global squamates species (Feldman et al., 2016), the maximum body length and body mass of snakes have a wider ranger than that of lizards, and are obviously larger than lizards. Therefore, more genes may be required to be under positive selection to participate in the body size differentiation of snakes and regulate the process of growth and development. In addition, the global ω values of body-size-related genes in snakes was generally higher than that in lizards, which was another important evidence for the above inference.

Most notably, although body-size-related genes in snakes evolved at a faster rate overall, lizards and snakes also shared the same positively selected genes, COL10A, GHR, and NPC1. It has been reported that the increased and decreased expression of these three genes and mutations in genes all could cause the changes (overgrowth or dwarfism) in body size. Similarly, positive selection signals of COL10A, GHR, and NPC1 in our study provided the molecular evidence for their roles of regulating the adaptive evolution of body size in squamates. On the one hand, COL10A1 is overexpressed in a variety of tumors and plays a pivotal role in cell proliferation, migration, invasion, and tumor vasculature (Chapman et al., 2012; Li et al., 2018); on the other hand, the amino acid substitution G590R in the alpha1 (X) chain of type X collagen, which is encoded by COL10A1, resulted in dwarfism and growth plate abnormality in Sus scrofa (Nielsen et al., 2000). Moreover, studies in several vertebrate species have clearly demonstrated that GH/IGF-1 pathway is a major controller for growth rate and body size (Jia et al., 2018). GH binds to GHR to exert biological functions, thereby activating a series of intracellular pathways; GHR deficiency leads to growth and metabolic disturbances, whereas its increased expression is associated with various diseases, such as cancer (Dehkhoda et al., 2018). In addition, NPC1 protein plays an important role in essential metabolic activities during development, such as subcellular lipid transport, the balance of body weight, the function and formation of platelet, and its mutation is related to normal weight gain (Meyre et al., 2009), and NPC1-null mice showed delayed weight loss and decreased food intake (Xie et al., 1999).

The body size of squamates species varies significantly. This study respectively divided the representative species of lizards and snakes into large and small bodies, and then explored the evolutionary patterns of lizards and snakes by PMAL v4.9. Branch models displayed that GALNS was subjected to significantly different selective pressure in lizards and snakes with different body size. GALNS is a lysosomal enzyme that degrades keratan sulfate (KS) and chondroitin-6-sulfate (C6S), and its deficiency leads to dwarfism with abnormal bone and connective tissue (Yasuda et al., 2013). In our study, the evolutionary rate of GALNS in small-bodied lizards (ω = 0.118) was higher than that in large-bodied lizards (ω = 0.070), and similarly, the evolutionary rate of GALNS in small-bodied snakes (ω = 0.261) was higher than that in large-bodied snakes (ω = 0.097), which indicated that the activity of GALNS enzymes in small-bodied lizards and snakes might be reduced, and the accumulation of KS and C6S provided the possibility for the body size reduction of squamates. This result is in contrast to what was found in cetaceans, where GALNS of large-bodied cetaceans had a higher evolutionary rate than that of small-bodied cetaceans (Sun, 2017). These means that the change of evolutionary rate in GALNS could affect the activity of enzymes, and thus regulate the evolution of body size in vertebrates.

Spearman correlation analyses showed the evidence of the high correlation of ω values for large lizards and large snakes (Spearman’s ρ = 0.730; p < 0.01), suggesting the similar evolutionary patterns in the process of evolving into larger body size. This also provided a molecular basis for the partition in further selective pressure analyses, where large-bodied species were set as the foreground branches. In the lizard-only and snake-only datasets, CDKN2C was under divergent selection between the lineages of large-bodied species and small-bodied species, and exhibited the relaxation of purifying selection (ω_foreground > ω_background) in the lineages of large-bodied species, whereas CDKN2C was only positively selected in lizards (ω_foreground > 1). It has been reported that cell proliferation is tightly regulated by the dosage balance between cyclin-dependent kinases (CDK) and CDK inhibitors (CKI), which are critical cell cycle regulators) (Hubbi and Semenza, 2015). CDKN2 family specifically inhibits CDK4 and CDK6, resulting in their inability to bind and activate cyclin D, and then arresting in G1-phase (Sherr and Roberts, 1999). Expression data analyses of the pleiotropic genes obtained from genome-wide association studies (GWAS) revealed that CDKN2C was significantly correlated with body mass index (BMI) and insulin resistance (p = 1.99E-12 and p = 6.27E-11), suggesting its roles in type 2 diabetes (T2D) and obesity (Zeng et al., 2021). In addition, genomic analysis of Xenopus found that CDKN2C was differentially expressed between chromosomes L and S, and the conserved non-coding sequence (CNS) of CDKN2C genomic loci on X. laevis and X. tropicalis ChrS had enhanced activity in regulating expression levels (Tanaka et al., 2017), affecting the development and health of animals. In a word, the expression of CDKN2C was related to body mass, and the variability and fragility leading to the generation of new functions and expression patterns were discovered in the relative species of reptiles, which all improves the possibility of adjusting the body size of reptiles.

We examined the positive signal of FBN1 in the lineages of small-bodied and large-bodied snakes. FBN1 encodes fibrillin-1 with 350 kDa in length that has positive effect on the integrity and functions of all connective tissues (Sakai et al., 2016). The mutations of FBN1 are diverse, up to 2700 (Caputi et al., 2002), which could cause abnormal structure of microfibrils and even diseases, including Marfan syndrome (MFS), Weill-Marchesani syndrome (WMS) and so on. Mutations at different loci of FBN1 may lead to different phenotypes. For example, the major feature of Marfan syndrome is usually tall stature. However, some studies have found that only mutations in the TGFβ binding-protein-like domain 5 (TB5) of FBN1 triggered acromelic dysplasia with severe short stature phenotypes opposite of MFS, such as geleophysic dysplasia (GD), acromicric dysplasia (AD), WMS (Le Goff et al., 2011). In the present study, we detected a positively selected site (PP > 0.8) in the terminal branch of the small-bodied hot-spring snakes T. baileyi and thirteen positively selected sites (PP > 0.8) in the united branch of large-bodied snakes O. hannah and P. bivittatus. There were no overlapping positive selection sites between these two branches. Consistent with the descriptions in the literature, FBN1 was also found to play an important role in the increase and decrease of body size in cetaceans (Sun, 2017), and in our study this gene was positively selected in both large-bodied and small-bodied species of snakes, and sites under positive selection were completely different.

Our study also found that LCORL has higher global ω values in both lizards and snakes, and was positively selected along the lineages to the LCA of larger-bodied species (the LCA of P. mucrosquamatus and P. flavoviridis, the LCA of P. guttatus and P. obsoletus). LCORL contains a characteristic motif of transcription factors, and is considered to play a significant role during spermatogenesis in the testes (Lindholm-Perry et al., 2013); analyses of tissues from mice indicated that LCORL are capable of activating transcription (Metzger et al., 2013). In order to realize various functions with the same DNA-binding activity, LCORL may interact with other functional proteins to alter or specialize their functions or target genes (Kunieda et al., 2003). In addition, GWAS analyses showed that the polymorphisms of LCORL was significantly associated with the length of human trunk and hip axis, and the rate of height growth in children (Soranzo et al., 2009). The functions of LCORL, particularly similar to a transcription factor, explain its roles in regulating body size evolution in squamates. However, further association analyses with squamate phenotypic data are required, and we also need to explore how this gene works by a series of function experiments.

Genes (AKT1, BMP1, IGF1, SOX5, SOX7) whose evolutionary rates were significantly correlated with body length and body mass were exactly same, whereas only the evolutionary rates of BMP5, BMP7, GPC6, SH2B3 and SOX17 were strongly associated with both body length and body mass in snakes. Based on the results of selection pressure analyses, it was found that BMP1 in lizards and SH2B3 in snakes were not only genes that were significantly positively selected along the branches of large-bodied species, but also genes that were significantly negatively correlated with body length and body mass. BMPs are important growth regulators of embryogenesis and tissue homeostasis in the adult organism (Correns et al., 2021). Osteogenesis imperfecta (OI) is an inherited bone disease caused by gene mutations characterized by increased bone fragility and recurrent fractures, and more than 20 causative genes including BMP1 have been identified so far (Correns et al., 2021). BMP1/Tolloid (TLD) is the prototype of the metalloproteinases family, whose functions include the proteolytic removal of the carboxy-terminal propeptide in procollagens I, II, and III and the amino-terminal propeptide in procollagens V and XI and the regulation of embryonic development in different species (Xi et al., 2021). In a study of OI patients lacking BMP1/mTLD, delayed cleavage of type I collagen C-propeptide, disassembly of type I/V collagen fibrils, and impaired processing of small leucine-rich proteoglycans (SLRP) were observed (Pollitt et al., 2016). Hence, the evolutionary rate of BMP1 in lizards was significantly negatively correlated with body length (R² = 0.595, p = 0.002) and body mass (R² = 0.571, p = 0.003), and this gene is subject to positive selection, inferring that lizards might strengthen the function of BMP1 protein to regulate the bone development, because bone mass is also a measure of body size.

As a member of the SH2B adaptor protein family, SH2B3 has multiple functional domains, containing the SH2 domain necessary for binding to and repressing target proteins, and the PH domain that distinguishes and translocates phosphoinositide to the cell membrane, and participates in the signal transduction of regulating various cytokine signaling cascade (Blass et al., 2016). On the one hand, SH2B3 dysregulation is linked to aging, cardiovascular disease, myeloproliferative cancers and so on; on the other hand, SH2B3 binds and regulates several key signaling pathways negatively, such as cell proliferation and differentiation (Kuo et al., 2020). In the snake-only dataset, the evolutionary rate of SH2B3 was significantly negatively associated with body length (R² = 0.560, p = 0.001) and body mass (R² = 0.667, p < 0.001), and was positively selected overall. The evidence of positive selection signals was examined at the united branch of large-bodied species in snakes and at the ancestral branch of Gekkonidae in lizards, suggesting that SH2B3 might guide body size to increase and decrease by activating or inhibiting growth and development related signaling pathways in squamates. Nevertheless, in terms of the genes and sites under positive selection detected in our study, functional experiments are still needed to verify their specific mechanism.

Different ecological variables have different effects on lizards and snakes. We found that the influence of climate on body size evolution in lizards and snakes was not significant, which was consistent with previous study (Slavenko et al., 2019). Slavenko et al. (2019) found that body size evolution in global squamate reptiles was independent of multiple climatic components. In addition, no correlation between microhabitat and body length/body mass in the lizard-only and snake-only datasets was discovered (p > 0.05), but there was a weak correlation between microhabitat and body mass in the squamate dataset, which were generally different from previous studies. Aquatic lizards and snakes are usually larger than terrestrial lizards and snakes, and arboreal snakes are lighter than terrestrial snakes (Meiri, 2008; Feldman and Meiri, 2013). Actually, in squamates, diversification rate was more closely related to microhabitat than climate, which explained 37% of the variation in diversification rate between clades (Bars-Closel et al., 2017); similarly, strong effect of microhabitat on diversification rate was found in frogs (Moen and Wiens, 2017). Therefore, we think that the divergence between our results and previous studies is related to the limited sequenced squamate genomes and poor quality of genome assembly. In the future, for one thing, with the subsequent detection of reptilian genomes, we could add more species into selective analyses and association analyses, such as crocodiles, turtles and new squamates; for another thing, more ecological factors could be included in to further explore the general ecological laws.

In our study, the effects of ecological variables on body size of squamates were not obvious, but it was found that ecological variables had a greater impact on the evolutionary rate of body-size-related genes. In terms of the number of significantly correlated genes, the genes correlated to climatic-niche rate for PC2 were the most, followed by climatic-niche rate for PC3 and Bio1. From the best fitting models of genes significantly affected by multiple ecological variables, most of the models contained five climate-related variables, range size (number of ecoregions), and rate of range expansion in species (ecoregions/age), and the contribution of climate-related variables to the variation in evolutionary rate of body-size-related genes were relatively higher. Although studies have found that climate is not the major determinant of the body size evolution in squamates (Slavenko et al., 2019), many literatures show that body size of reptiles is closely related to climate. For instance, reptiles prefer to distribute in hot and arid regions (Roll et al., 2017); contrary to Bergman’s rule, turtles at lower latitudes generally exhibit larger body size (Angielczyk et al., 2015). Therefore, these also explain why our study discovered significant effects of climate-related variables on the evolutionary rate of body-size-related genes at the molecular level. Interestingly, SOX30 in the snake-only dataset was the gene affected by the largest number of ecological variables, and its best-fitting model contained two variables: number of ecoregions and Bio1 (annual mean temperature), which respectively explained ∼48% and ∼28% of the variation in the evolutionary rate. Furthermore, SOX30 also was detected under positive selection along the united branch of large-bodied snakes, and its evolutionary rate was significantly negatively correlated with body length in snakes, which were consistent with its important function to inhibit cell proliferation (Hao et al., 2018; Liu et al., 2020). As a master switch of desmosomal genes, SOX30 suppressed lung adenocarcinoma cell proliferation, migration and invasion by activating the transcription of desmosomal genes (Hao et al., 2018). Overexpression of SOX30 significantly reduces the expression of β-catenin, resulting in the inactivation of the Wnt/β-catenin pathway (Liu et al., 2020), which is a critical pathway for regulating the animal body size. All in all, for squamate reptiles, ecological factors affect the evolutionary rate of body-size-related genes at the molecular level and body size at the physiological level. Due to the organism’s self-regulation mechanism, small changes in ecological factors may not be reflected in phenotype, but may cause the variation in the evolutionary rate of body-size-related genes, thereby controlling the synthesis of related proteins and changing the activity of related enzymes in response to environmental changes.