Data analysis
We used R 4.0.2. for data analysis (R Core Team 2017). First, we reduced the dimensionality of the climatic variables using principal component analysis (PCA). Then we performed a phylogenetically controlled analysis to examine the relationship between size, climatic factors (PC1 and PC2; see results), and the mean butterfly reflectance over VIS-NIR range. The phylogenetic relationships of the analysed species were inferred from a published time-calibrated molecular phylogenetic tree (Wiemers et al. 2020). To account for the uncertainties in topology and node age, we generated 1,000 trees randomly sampled from the posterior distribution and performed the analysis using all 1000 trees.
For each tree, we fitted multivariate phylogenetic regressions implemented in the ‘mvMORPH’ package (Clavel et al. 2015). This enabled us to include multiple response variables in a single model. We set size, climatic principal components (PCs; up to PC2), and the interaction between size and climatic PCs as predictor variables. Our response variables consisted of the mean reflectance of the six body regions (DT, DB, DE, VT, VB, DE). We compared the goodness of fit of the models assuming either Brownian motion, Ornstein-Uhlenbeck, Early Burst, or Pagel’s λ models of trait evolution and used the models with the lowest Generalized Information Criterion (GIC) (Konishi & Kitagawa 1996; Hernández et al. 2013). Ornstein-Uhlenbeck models showed slightly better fit than other models in multivariate phylogenetic regressions, but Pagel’s λ models outperformed others in allpost-hoc PGLS models. Thus, we consistently assumed Pagel’s λ models of trait evolution. Nevertheless, the results were robust regardless of which models we assumed. We also used GIC to select the best model among all candidate models.
Using the predictors that remained significant in the above model, we performed post-hoc phylogenetic generalised least squares (PGLS) on the mean reflectance of each body region to further examine which body regions were specifically explained by each predictor. We implemented a ‘gls’ function to run PGLS (Revell 2012). We iteratively performed PGLS on 1,000 trees for all body regions and estimated the 95% confidence intervals of statistics. The 95% confidence limits of statistics from all analyses using 1000 trees varied only slightly and showed essentially the same results as the MCC tree results (all parameters estimated were within ± 0.02 range) because of low phylogenetic uncertainty. Thus, we report the results from the MCC tree only in the results. Here, we used Akaike Information Criteria (AIC) to select the best model and controlled for the false-discovery rate by adjusting P-values (Benjamini & Hochberg 1995). The strength of phylogenetic signal λ was estimated using the maximum clade credibility (MCC) tree.
We additionally examined the effects of our predictors on VIS reflectance and NIR reflectance separately using the same modelling approach and model structure. To examine whether NIR reflectance showed adaptive features even after accounting for its high correlation with VIS reflectance, we fitted linear regressions to each NIR reflectance (for each body region) using log form of VIS reflectance as a predictor (average r2 among all body regions ≈ 0.77) and extracted residuals of the fitted models. These residuals indicate the degree of NIR reflectance independent of VIS reflectance (i.e. after accounting for the correlation between VIS and NIR reflectance). Though the use of residuals has been criticised due to potential for biased parameter estimation (Freckleton 2009), residuals are appropriate when the effect of one predictor on a response variable should be controlled ‘before and independent’ of the other predictors such as in our case; the effects of VIS reflectance on NIR reflectance should be considered before estimating the effects of other climatic and size variables and should not affect the parameters of other predictors in the model (see Supplementary materials for the statistical justification of using residuals). We fitted multivariate phylogenetic regressions using size, climatic PCs, and the interaction between those two as predictors and the residuals for each body region as response variables. Then we performed a post-hoc PGLS analysis for each body region using the retained significant terms as predictors as above.
In the main results, we analysed all specimen images regardless of the sex and the presence of sexual dimorphism. However, because sexually dimorphic male reflectance is likely to be heavily shaped by sexual selection (Lande 1980; van der Bijl et al. 2020), we repeated the analysis without sexually dimorphic male specimens and compared the results.
For the analysis of ventral-dorsal reflectance differences, we first compared the reflectance between dorsal and ventral parts using paired t-tests. Next, to examine whether climate variables predict these differences, we calculated the reflectance differences between ventral and dorsal parts by subtracting dorsal reflectance from ventral reflectance for each body region. Using these ventral-dorsal differences in the three body regions as independent variables, we fitted phylogenetic multivariate multiple regressions with climatic PCs, size, and the interaction between them as dependent variables. We then performed post-hoc PGLS analyses for each body region using the retained significant terms as predictors.