Colour variation in the Italian wall lizards was assessed in a three-step analysis (Sacchi et al., 2021). First, we generated the individuals’ frequency distributions of hue, saturation, and brightness by using values computed on the whole samples of pixels selected for each individual (HSV colour spectrum). Second, three principal component analyses (PCA), one each for H, S, and V, respectively, were used on the colour spectra, and the first components, explaining 30.7%, 49.0% and 67.2% of the total variance, were used as a proxy to summarise the inter-individual variability of colouration. The hue PC score accounted for the opposite variation of the orange-yellow and green interval in the hue colour spectrum; negative scores were associated with a higher peak in the orange-yellow interval whereas positive scores were associated with a higher peak in the green interval. The saturation PC score accounted for the increase of colour saturation with an increased score, whereas the value PC score accounted for the increase of brightness with higher scores (Fig. 6). Third, the PC scores were analysed through random intercept linear mixed models (LMM) including a single-component cosinor function to model the effect of the season (Refinetti et al., 2007; Cornelissen, 2014). Cosinor models were originally developed to model circadian rhythm in physiological processes (Halberg et al., 1967), but they can be also used to model periodic variations in ecological variables (Mangiacotti et al., 2019; Sacchi et al., 2020). In Cosinor models the response variable (Y) is assumed to depend on time (t) following a regular cycle, which is incorporated in a linear model through a cosine function:where M is the MESOR (Midline Statistic Of Rhythm, i.e. the time-corrected mean of the response), A is the amplitude (maximum absolute deviation from MESOR), τ the period of the cycle, the acrophase (i.e. the timing of highest values), and e(t) the error term (Cornelissen, 2014). The model can be linearized by rewriting the formula: Y(t)=M+βx+γz+e(t); being x=cos(2πt/τ) and z=sin (2πt/τ) the cosinor terms, and β=Acos and γ=− Asin the cosinor coefficients (Cornelissen, 2014). In our model, the two cosinor terms entered the LMM as fixed effects, with time expressed as Julian date (1= 1 January) and τ=365 to account for circannual rhythms around the time-corrected mean of H, S, V values. Additional fixed effects were sex and body size (i.e. standardized SVL). The two-way interactions sex×size, sex×x, and sex×z were also added to account for possible differential effect of sex on size and season on colour expression. Population entered the model as a random effect. The dependent variables were the PC scores of H, S and V values, and three independent models were run. LMMs were fit in a Bayesian analytical framework available in the package JAGS 4.3.0 (http://mcmc-jags.sourceforge.net/), using flat priors for coefficients and intercept (μ=0 and σ=0.001), and uninformative half-Cauchy priors (x0=0, γ=25) for both σ²error and σ²population. For all models, Markov Chain Monte Carlo parameters were set as follows: number of independent chains=three; number of iterations=34,000; burning=4000; thinning=three. Convergence was checked and results from the posterior distribution are reported as the half sample mode (HSM, Bickel and Frühwirth, 2006) plus 95% (or 50%) highest density intervals (HDI95; HDI50; Kruschke, 2010). In Bayesian statistics, the HSM is a commonly used estimator of the central tendency of posterior probability distribution robust to outliers, whereas the HDI95 defines the interval that includes the parameter with 95% probability. Parameter values in the centre of the HDI have higher credibility than parameter values at the limits of the interval. Therefore, when the HDIs of two groups do not overlap, there is a credible evidence for different group means. By contrast, to the extent of the two groups’ HDIs overlap there is no credible evidence of difference between the means. When comparing two groups (e.g. males and females), we therefore reported in Table 1 (bold) the posterior probability of their difference being different (i.e. higher or lower) from zero. Data preparation, model settings, call to JAGS, and posterior elaborations were done in R 4.0 using the package R2jags (Su and Yajima, 2015), modeest (Poncet, 2012), and HDInterval (Meredith and Kruschke, 2018).