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.