Statistical analyses

For each species, we quantified SSD with the Gibbons and Lovich index [59], SSD (Sexual Size Dimorphism) = (SVL of longest sex / SVL of shortest sex) – 1. The resulting value is then made negative if males are the larger sex and positive if females are the larger sex [3]. To evaluate the direction and magnitude of sexual dimorphism, a permutations test was performed on SSD, with p-value < 0.05, to determine if the sexual dimorphism index of each species was statistically significant. Those species that equaled or exceeded 15% of the SSD Index of their own family were also considered dimorphic, despite not having a significant p-value in the permutation test. We categorized the species into 3 SSD groups: male-biased SSD (when sexual dimorphism was deviated towards males); female-biased SSD (when sexual dimorphism was deviated towards females) and monomorphic (when they had no marked sexual dimorphism).

We calculated the mean values of the morphological variables per species per sex and we used Kolmogorov–Smirnov (KS) tests to ensure normality. The variables without normal distribution were logarithmically (log₁₀) transformed priori to analysis. As an estimate of intrasexual selection of reproduction-related phenotypic traits, we calculated an Index sexual dimorphism to HW, HH, LH and TP_male, as the ratio of HW in males to HW in females (IHW=HW_males /HW_females) (see [17, 60, 61]), and the same for the other variables. To estimate the fecundity selection of reproduction-related phenotypic traits, we calculated an Index sexual dimorphism to TL, AW and TP_female, as the ratio of TL in females to TL in males (ITL = TL_females /TL_males), and the same for the other variables.

When examining data from phylogenetically related species, data points cannot be considered as statistically independent due to shared evolutionary history [58, 62]. So we performed the phylogenetic size-correction analysis [63] by using phylo.resid (a module of Phytools for R developed by Revell [64]), over SSD, IHW, IHH, ILH, ITP_male, ITL, IAW, ITP_female and clutch/litter size. The resultant residuals from the phylogenetic size-correction were then used in the subsequent analyses.

To analyze the relationship between SSD and body size of the species, we ran Phylogenetic Generalized Least Squares (PGLS) using a model with SVLLog₁₀ as predictor variable and SSD as dependent variable. To analyze the effect of sexual dimorphism on reproduction-related phenotypic traits, PGLS was run using models with residuals of IHW, IHH, ILH, ITP_male and of ITL, IAW, ITP_female as dependent variables and the residual SSD and SSD groups (males-biased SSD, females-biased SSD and monomorphic) as predictor variables.

To analyze the effect of sexual dimorphism and female reproduction-related phenotypic traits (only those significantly related to residual SSD according to the previous PGLS) on fecundity, we ran PGLS using models with residual clutch/litter size as the dependent variable and the residual SSD and SSD groups (males-biased SSD, females-biased SSD) as predictor variables. We also ran PGLS using models with residual clutch/litter size as the dependent variable and the residual ITL and SSD groups (males-biased SSD, females-biased SSD) as predictor variables. In these analyses, we eliminated the monomorphic species group because there was not enough data.

To analyze the effect of sexual dimorphism on female reproduction-related phenotypic traits (only those significantly related to residual SSD according to the previous PGLS) considering the species’ reproductive mode, we ran PGLS using models with the residual ITL as dependent variable and the residual of SSD and the reproductive mode (oviparous and viviparous) as predictor variables. Also, to analyze the effect of female reproduction-related phenotypic traits of (only those significantly related to residual SSD according to the previous PGLS) on fecundity, we ran PGLS using models with the residual clutch/litter size as dependent variable and the residual of ITL and the reproductive mode (oviparous and viviparous) as predictor variables.

We estimated Pagel’s phylogenetic signal (λ) from the residual errors simultaneously on the regression parameters of phylogenetic generalized least squares models (PGLS) analyses. Analyses were made in ‘caper’ [65] and ‘ape’ [66] packages, both developed in R [67].