The aim

  • To test whether the modelling results were sensitive to the photoreceptor ratios used, we repeated these analyses using published ratios that represent the known variation in insects. Specifically, we compared the changes in contrast when shifting the LWS peak from 580 nm to 660 nm between tetrachomat visual systems with different photoreceptor ratios.

  • We applied the same modelling method described in the Main models but using ratios of different insects, including the jewel beetle, Agrilus planipennis (same as the ratio used in the main text; 1.14: 1: 1.26: 1.38; UVS: SWS: MWS: LWS), and two butterflies, Papilio Xuthus (1: 1: 4.08: 2.92) and Heliconius sp. type III (1.00: 1.44: 2.22: 11.11).

  • The model parameters and the statistical methods remain the same as described in main text of the paper.

library(pavo)
library(dplyr)
library(stringr)
library(tidyr) #for gather() function
library(ggplot2)
library(lme4)
library(car)
library(multcomp)
library(boot) #for mean() function
library(pander) #for creating tidy tables
library(ggpubr) #for ggarrange() function
wl <- specsensbuprest.model2[,1]
peaks <- gather(specsensbuprest.model2[,2:9], peak, value) %>% 
  cbind(wl)

#order the peaks in the legend
peak.order <- c("UVS.355.A1.", "SWS.445.A1.", "MWS.530.A1.", "LWS.570.A1..filter580.", "LWS.570.A1..filter600.", "LWS.570.A1..filter620.", "LWS.570.A1..filter640.", "LWS.570.A1..filter660.") 

ggplot(peaks,
       aes(x = wl, 
                 y = value, 
                 col = peak))+
  geom_line()+
  guides(color = guide_legend(title = "peak sensitivity"))+
  scale_color_manual(
    values = c("darkorchid4", "dodgerblue3", "olivedrab4", "orange1", "orange3","darkorange3", "orangered1", "red2"),
    labels = c("355 nm", "445 nm","530 nm", "580 nm", "600 nm", "620 nm", "640 nm", "660 nm" ),
    breaks = peak.order)+
  xlab("Wavelength (nm)")+ 
  ylab("Relative spectral sensitivity")+
  theme_classic()
Figure caption: Sensitivity curves of the photoreceptors with different peak wavelengths.

Figure caption: Sensitivity curves of the photoreceptors with different peak wavelengths.



Data description

We used the same spectral data and daylight illumination described in the Main models section.



Run the visual models

  • We created 3 different visual systems using different receptor ratios from 3 insect species (UVS, SWS, MWS, LWS):

  • Agrilus planipennis: 1.14, 1, 1.26, 1.38

  • Papilio xuthus: 1, 1, 4.08, 2.92

  • Heliconius sp. type III : 1, 1.44, 2.22, 11.11


  • We ran the visual models with the same parameter settings and steps as in the Main models.

  • We first calculated the quantum catches, then calculated the contrasts.

get.d65.vismodel <- function(i){
  
  vs.i <- vismodel(dataset[1:501,], 
                   visual = i, #this need to change according to the visual system
                   bkg = aveleaf$aveleaf, 
                   illum = irradiance.d65[1:501,2], 
                   qcatch = 'fi', 
                   relative = FALSE,
                   vonkries = TRUE)
  return(vs.i)
}
#shared between 3 species
sens.list <- list(VS580 <- specsensbuprest.model2[, 1:5],
                  VS600 <- specsensbuprest.model2[, c(1,2,3,4,6)],
                  VS620 <- specsensbuprest.model2[, c(1,2,3,4,7)],
                  VS640 <- specsensbuprest.model2[, c(1,2,3,4,8)],
                  VS660 <- specsensbuprest.model2[, c(1,2,3,4,9)])
d65.vismodel.output <- list(VS580 = NA, VS600 = NA, VS620 = NA, VS640 = NA, VS660 = NA)

# Set a loop for calculating quantum catch for each visual system
for (i in 1:length(sens.list)){
  
  vs.result.d65 <- get.d65.vismodel(sens.list[[i]]) # for D65
  d65.vismodel.output[[i]] <- vs.result.d65
}

buprest580 <- d65.vismodel.output$VS580
buprest600 <- d65.vismodel.output$VS600
buprest620 <- d65.vismodel.output$VS620
buprest640 <- d65.vismodel.output$VS640
buprest660 <- d65.vismodel.output$VS660
recep.density.list <- list(Agrilus <- c(1.14, 1, 1.26, 1.38),
                                Papilio <- c(1, 1, 4.08, 2.92),
                                Heliconius <- c(1, 1.44, 2.22, 11.11))

get.sens.coldist <- function(i, sp){
  
  con <- coldist(modeldata = i, # put output of vismodel()
                   noise = "neural", 
                   achro = FALSE, 
                   n = recep.density.list[[sp]],
                   weber = 0.12,
                   weber.ref = 4)
  
  return(con)
}

vismodel.list <- list(buprest580, 
                      buprest600, 
                      buprest620, 
                      buprest640, 
                      buprest660)
#for Agrilus planipennis
Agrilus.contrast.output <- list(VS580 = NA, VS600 = NA, VS620 = NA, VS640 = NA, VS660 = NA)

for (i in 1:length(vismodel.list)) {
  
  contrast.result.i <- get.sens.coldist(vismodel.list[[i]], 1)
  Agrilus.contrast.output[[i]] <- contrast.result.i
}

Cbuprest580 <- Agrilus.contrast.output$VS580
Cbuprest600 <- Agrilus.contrast.output$VS600
Cbuprest620 <- Agrilus.contrast.output$VS620
Cbuprest640 <- Agrilus.contrast.output$VS640
Cbuprest660 <- Agrilus.contrast.output$VS660


#for Papilio xuthus
Papilio.contrast.output <- list(VS580 = NA, VS600 = NA, VS620 = NA, VS640 = NA, VS660 = NA)

for (i in 1:length(vismodel.list)) {
  
  contrast.result.i <- get.sens.coldist(vismodel.list[[i]], 2)
  Papilio.contrast.output[[i]] <- contrast.result.i
}

Cpapilio580 <- Papilio.contrast.output$VS580
Cpapilio600 <- Papilio.contrast.output$VS600
Cpapilio620 <- Papilio.contrast.output$VS620
Cpapilio640 <- Papilio.contrast.output$VS640
Cpapilio660 <- Papilio.contrast.output$VS660

#for Heliconius sp
Heliconius.contrast.output <- list(VS580 = NA, VS600 = NA, VS620 = NA, VS640 = NA, VS660 = NA)

for (i in 1:length(vismodel.list)) {
  
  contrast.result.i <- get.sens.coldist(vismodel.list[[i]], 3)
  Heliconius.contrast.output[[i]] <- contrast.result.i
}

Cheliconius580 <- Heliconius.contrast.output$VS580
Cheliconius600 <- Heliconius.contrast.output$VS600
Cheliconius620 <- Heliconius.contrast.output$VS620
Cheliconius640 <- Heliconius.contrast.output$VS640
Cheliconius660 <- Heliconius.contrast.output$VS660
#Organize data before GLMM
fl.vissys <- list("VS580.fl", "VS600.fl", "VS620.fl", "VS640.fl", "VS660.fl") 
bl.vissys <- list("VS580.bl", "VS600.bl", "VS620.bl", "VS640.bl", "VS660.bl")
bf.vissys <- list("VS580.bf", "VS600.bf", "VS620.bf", "VS640.bf", "VS660.bf") 


# for Agrilus planipennis
Agrilus.vissys <- list(Cbuprest580,
                       Cbuprest600, 
                       Cbuprest620, 
                       Cbuprest640, 
                       Cbuprest660)
##flower vs leaf
allvis.fl.Agrilus <- data_frame()
for (i in 1:length(Agrilus.vissys)) {
  
  temp.i <- Agrilus.vissys[[i]] %>% 
    filter(str_detect(patch1,"flower")) %>% 
    filter(str_detect(patch2,"leaves")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(fl.vissys[[i]],1)
  
  allvis.fl.Agrilus <- temp.i %>% rbind(allvis.fl.Agrilus)  
}

### beetle vs leaf
allvis.bl.Agrilus <- data_frame()
for (i in 1:length(Agrilus.vissys)) {
  
  temp.i <- Agrilus.vissys[[i]] %>% 
    filter(str_detect(patch1,"leaves")) %>% 
    filter(str_detect(patch2,"beetle")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(bl.vissys[[i]],1)
  
  allvis.bl.Agrilus <- temp.i %>% rbind(allvis.bl.Agrilus)  
}

##beetle vs flower
allvis.bf.Agrilus <- data_frame()
for (i in 1:length(Agrilus.vissys)) {
  
  temp.i <- Agrilus.vissys[[i]] %>% 
    filter(str_detect(patch1,"flower")) %>% 
    filter(str_detect(patch2,"beetle")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(bf.vissys[[i]],1)
  
  allvis.bf.Agrilus <- temp.i %>% rbind(allvis.bf.Agrilus)  
}

#for Papilio xuthus
Papilio.vissys <- list(Cpapilio580,
                       Cpapilio600, 
                       Cpapilio620, 
                       Cpapilio640, 
                       Cpapilio660)
##flower vs leaf
allvis.fl.Papilio <- data_frame()
for (i in 1:length(Papilio.vissys)) {
  
  temp.i <- Papilio.vissys[[i]] %>% 
    filter(str_detect(patch1,"flower")) %>% 
    filter(str_detect(patch2,"leaves")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(fl.vissys[[i]],1)
  
  allvis.fl.Papilio <- temp.i %>% rbind(allvis.fl.Papilio)  
}

##beetle vs leaf
allvis.bl.Papilio <- data_frame()
for (i in 1:length(Papilio.vissys)) {
  
  temp.i <- Papilio.vissys[[i]] %>% 
    filter(str_detect(patch1,"leaves")) %>% 
    filter(str_detect(patch2,"beetle")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(bl.vissys[[i]],1)
  
  allvis.bl.Papilio <- temp.i %>% rbind(allvis.bl.Papilio)  
}

##beetle vs flower
allvis.bf.Papilio <- data_frame()
for (i in 1:length(Papilio.vissys)) {
  
  temp.i <- Papilio.vissys[[i]] %>% 
    filter(str_detect(patch1,"flower")) %>% 
    filter(str_detect(patch2,"beetle")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(bf.vissys[[i]],1)
  
  allvis.bf.Papilio <- temp.i %>% rbind(allvis.bf.Papilio)  
}

#for Heliconius sp
Heliconius.vissys <- list(Cheliconius580,
                       Cheliconius600, 
                       Cheliconius620, 
                       Cheliconius640, 
                       Cheliconius660)
##flower vs leaf
allvis.fl.Heliconius <- data_frame()
for (i in 1:length(Heliconius.vissys)) {
  
  temp.i <- Heliconius.vissys[[i]] %>% 
    filter(str_detect(patch1,"flower")) %>% 
    filter(str_detect(patch2,"leaves")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(fl.vissys[[i]],1)
  
  allvis.fl.Heliconius <- temp.i %>% rbind(allvis.fl.Heliconius)  
}

##beetle vs leaf
allvis.bl.Heliconius <- data_frame()
for (i in 1:length(Heliconius.vissys)) {
  
  temp.i <- Heliconius.vissys[[i]] %>% 
    filter(str_detect(patch1,"leaves")) %>% 
    filter(str_detect(patch2,"beetle")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(bl.vissys[[i]],1)
  
  allvis.bl.Heliconius <- temp.i %>% rbind(allvis.bl.Heliconius)  
}

##beetle vs flower
allvis.bf.Heliconius <- data_frame()
for (i in 1:length(Heliconius.vissys)) {
  
  temp.i <- Heliconius.vissys[[i]] %>% 
    filter(str_detect(patch1,"flower")) %>% 
    filter(str_detect(patch2,"beetle")) %>% 
    dplyr::select(1:3)
  
  temp.i$vissys <- strrep(bf.vissys[[i]],1)
  
  allvis.bf.Heliconius <- temp.i %>% rbind(allvis.bf.Heliconius)  
}


Compare contrast



Click the tabs to see the results for each comparison group ( Flower vs. Leaf / Beetle vs. Leaf / Beetle vs. Flower) in 3 different speceis (Agrilus planipennis / Papilio xuthus / Heliconius sp).

get.lmer.sens <- function(inputdat){
  
  lmer(dS ~ vissys + (1|patch2) + (1|patch1), 
       data = inputdat, REML = F)
  #REML=F, in order to fit the model using the likelihood ratio test. Otherwise, the lmer default will fit the model using the REML (REstricted Maximum Likelihood) criterion.
}

get.posthocsum <- function(modelx){
  
  summary(glht(modelx, 
               linfct = mcp(vissys = "Tukey")), 
          test = adjusted("bonferroni"))
}
# for Agrilus planipennis
## Flower vs Leaf
m.flower.vs.leaf.bup <- get.lmer.sens(allvis.fl.Agrilus)
sum.flower.vs.leaf.bup <- get.posthocsum(m.flower.vs.leaf.bup)

##Beetle vs Leaf
m.beetle.vs.leaf.bup <- get.lmer.sens(allvis.bl.Agrilus)
sum.beetle.vs.leaf.bup <- get.posthocsum(m.beetle.vs.leaf.bup)

##Beetle vs Flower
m.beetle.vs.flower.bup <- get.lmer.sens(allvis.bf.Agrilus)
sum.beetle.vs.flower.bup <- get.posthocsum(m.beetle.vs.flower.bup)

#for Papilio xuthus
##Flower vs Leaf
m.flower.vs.leaf.pap <- get.lmer.sens(allvis.fl.Papilio)
sum.flower.vs.leaf.pap <- get.posthocsum(m.flower.vs.leaf.pap)

##Beetle vs Leaf
m.beetle.vs.leaf.pap <- get.lmer.sens(allvis.bl.Papilio)
sum.beetle.vs.leaf.pap <- get.posthocsum(m.beetle.vs.leaf.pap)

##Beetle vs Flower
m.beetle.vs.flower.pap <- get.lmer.sens(allvis.bf.Papilio)
sum.beetle.vs.flower.pap <- get.posthocsum(m.beetle.vs.flower.pap)

#for Heliconius sp
##Flower vs Leaf
m.flower.vs.leaf.heli <- get.lmer.sens(allvis.fl.Heliconius)
sum.flower.vs.leaf.heli <- get.posthocsum(m.flower.vs.leaf.heli)

##Beetle vs Leaf
m.beetle.vs.leaf.heli <- get.lmer.sens(allvis.bl.Heliconius)
sum.beetle.vs.leaf.heli <- get.posthocsum(m.beetle.vs.leaf.heli)

##Beetle vs Flower
m.beetle.vs.flower.heli <- get.lmer.sens(allvis.bf.Heliconius)
sum.beetle.vs.flower.heli <- get.posthocsum(m.beetle.vs.flower.heli)
vislist.heatmap <- list(
  c("VS 600", "VS 620","VS 640","VS 660","VS 620","VS 640","VS 660","VS 640","VS 660","VS 660"),
  c("VS 580","VS 580","VS 580","VS 580","VS 600","VS 600","VS 600","VS 620","VS 620","VS 640"))

# Set up the function
get.heatmap.datframe <- function(i,z){
  
  heatmapdat.i <- data.frame(as.numeric(str_extract(i[[1]],"([0-9]+).*$"))) %>%
    cbind(as.numeric(str_extract(i[[2]],"([0-9]+).*$"))) %>% 
    cbind(as.numeric(str_extract(i[[3]],"([0-9]+).*$"))) %>%
    dplyr::rename(flower.vs.leaf = 1, beetle.vs.leaf = 2, beetle.vs.flower = 3) %>% 
    cbind(z[[1]]) %>% 
    cbind(z[[2]]) 
  
  colnames(heatmapdat.i)[4:5] <- c("VislistA", "VislistB")
  
  return(heatmapdat.i)
}
# Create pvalue dataset
## for Agrilus planipennis
pdata.Agrilus <- list(sum.flower.vs.leaf.bup[["test"]][["pvalues"]],
                      sum.beetle.vs.leaf.bup[["test"]][["pvalues"]],
                      sum.beetle.vs.flower.bup[["test"]][["pvalues"]])

heat.Agrilus <- get.heatmap.datframe(pdata.Agrilus, vislist.heatmap)

## for Papilio xuthus
pdata.Papilio <- list(sum.flower.vs.leaf.pap[["test"]][["pvalues"]],
                      sum.beetle.vs.leaf.pap[["test"]][["pvalues"]],
                      sum.beetle.vs.flower.pap[["test"]][["pvalues"]])

heat.Papilio <- get.heatmap.datframe(pdata.Papilio, vislist.heatmap)

## for Heliconius sp
pdata.Heliconius <- list(sum.flower.vs.leaf.heli[["test"]][["pvalues"]],
                      sum.beetle.vs.leaf.heli[["test"]][["pvalues"]],
                      sum.beetle.vs.flower.heli[["test"]][["pvalues"]])

heat.Heliconius <- get.heatmap.datframe(pdata.Heliconius, vislist.heatmap)


# Assign asterisk signs
pdatalist <- list(heat.Agrilus, heat.Papilio, heat.Heliconius)

for (i in 1:length(pdatalist)) {
  
  pdatalist[[i]]$sig.flower.vs.leaf[pdatalist[[i]]$flower.vs.leaf > 0.05] <- ""
  pdatalist[[i]]$sig.flower.vs.leaf[pdatalist[[i]]$flower.vs.leaf < 0.05] <- "*"
  pdatalist[[i]]$sig.flower.vs.leaf[pdatalist[[i]]$flower.vs.leaf < 0.01] <- "**"
  pdatalist[[i]]$sig.flower.vs.leaf[pdatalist[[i]]$flower.vs.leaf < 0.0001] <- "***"
  pdatalist[[i]]$sig.beetle.vs.leaf[pdatalist[[i]]$beetle.vs.leaf > 0.05] <- ""
  pdatalist[[i]]$sig.beetle.vs.leaf[pdatalist[[i]]$beetle.vs.leaf < 0.05] <- "*"
  pdatalist[[i]]$sig.beetle.vs.leaf[pdatalist[[i]]$beetle.vs.leaf < 0.01] <- "**"
  pdatalist[[i]]$sig.beetle.vs.leaf[pdatalist[[i]]$beetle.vs.leaf < 0.0001] <- "***"
  pdatalist[[i]]$sig.beetle.vs.flower[pdatalist[[i]]$beetle.vs.flower > 0.05] <- ""
  pdatalist[[i]]$sig.beetle.vs.flower[pdatalist[[i]]$beetle.vs.flower < 0.05] <- "*"
  pdatalist[[i]]$sig.beetle.vs.flower[pdatalist[[i]]$beetle.vs.flower < 0.01] <- "**"
  pdatalist[[i]]$sig.beetle.vs.flower[pdatalist[[i]]$beetle.vs.flower < 0.0001] <- "***"
  
}

# Because we did the assignment on a data list, now we need to link the loop output back to the data respectively
heat.Agrilus <- pdatalist[[1]] 
heat.Papilio <- pdatalist[[2]]
heat.Heliconius <- pdatalist[[3]]
# a = data
# b = 1/2/3; 1 for flower.vs.leaf, 2 for beetle.vs.leaf, 3 for beetle.vs.flower
get.pheatmap <- function(a, b){ggplot(data = a, 
                                        aes(x = VislistA, 
                                            y = VislistB, 
                                            fill = a[,b])) + 
    geom_tile(colour = "white", size = 4)+
    geom_text(aes(VislistA, 
                  VislistB, 
                  label = paste(format(round(a[,b], 2), nsmall = 2), 
                                a[, b+5])))+
    scale_fill_continuous(high = "#132B43", low = "#56B1F7", limit = c(0,1))+ #delete if want to reverse the colour
    theme_bw()+
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          panel.grid = element_blank(), 
          panel.border = element_blank(), 
          axis.ticks = element_blank() )+
    labs(fill = "p vaalue")
  
}

Flower vs Leaf

Agrilus planipennis

Anova(m.flower.vs.leaf.bup) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 2268 4 0



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Agrilus, 1)



Original model output

sum.flower.vs.leaf.bup
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.fl - VS580.fl == 0  0.42097    0.04777   8.812   <2e-16 ***
## VS620.fl - VS580.fl == 0  1.17568    0.04777  24.609   <2e-16 ***
## VS640.fl - VS580.fl == 0  1.80399    0.04777  37.760   <2e-16 ***
## VS660.fl - VS580.fl == 0  1.76193    0.04777  36.880   <2e-16 ***
## VS620.fl - VS600.fl == 0  0.75471    0.04777  15.797   <2e-16 ***
## VS640.fl - VS600.fl == 0  1.38301    0.04777  28.949   <2e-16 ***
## VS660.fl - VS600.fl == 0  1.34096    0.04777  28.068   <2e-16 ***
## VS640.fl - VS620.fl == 0  0.62831    0.04777  13.151   <2e-16 ***
## VS660.fl - VS620.fl == 0  0.58625    0.04777  12.271   <2e-16 ***
## VS660.fl - VS640.fl == 0 -0.04206    0.04777  -0.880        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Papilio xuthus

Anova(m.flower.vs.leaf.pap) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 5425 4 0



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Papilio, 1)



Original model output

sum.flower.vs.leaf.pap
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.fl - VS580.fl == 0  0.57457    0.03751  15.319   <2e-16 ***
## VS620.fl - VS580.fl == 0  1.47910    0.03751  39.435   <2e-16 ***
## VS640.fl - VS580.fl == 0  2.19498    0.03751  58.521   <2e-16 ***
## VS660.fl - VS580.fl == 0  2.17608    0.03751  58.017   <2e-16 ***
## VS620.fl - VS600.fl == 0  0.90454    0.03751  24.116   <2e-16 ***
## VS640.fl - VS600.fl == 0  1.62041    0.03751  43.202   <2e-16 ***
## VS660.fl - VS600.fl == 0  1.60151    0.03751  42.698   <2e-16 ***
## VS640.fl - VS620.fl == 0  0.71587    0.03751  19.086   <2e-16 ***
## VS660.fl - VS620.fl == 0  0.69697    0.03751  18.582   <2e-16 ***
## VS660.fl - VS640.fl == 0 -0.01890    0.03751  -0.504        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Heliconius sp

Anova(m.flower.vs.leaf.heli) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 6518 4 0



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Heliconius, 1)



Original model output

sum.flower.vs.leaf.heli
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.fl - VS580.fl == 0  0.35075    0.02222  15.782   <2e-16 ***
## VS620.fl - VS580.fl == 0  0.94433    0.02222  42.491   <2e-16 ***
## VS640.fl - VS580.fl == 0  1.42640    0.02222  64.183   <2e-16 ***
## VS660.fl - VS580.fl == 0  1.39784    0.02222  62.898   <2e-16 ***
## VS620.fl - VS600.fl == 0  0.59358    0.02222  26.709   <2e-16 ***
## VS640.fl - VS600.fl == 0  1.07566    0.02222  48.400   <2e-16 ***
## VS660.fl - VS600.fl == 0  1.04710    0.02222  47.115   <2e-16 ***
## VS640.fl - VS620.fl == 0  0.48207    0.02222  21.691   <2e-16 ***
## VS660.fl - VS620.fl == 0  0.45351    0.02222  20.406   <2e-16 ***
## VS660.fl - VS640.fl == 0 -0.02856    0.02222  -1.285        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Beetle vs Leaf

Agrilus planipennis

Anova(m.beetle.vs.leaf.bup) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 2119 4 0



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Agrilus, 2)



Original model output

sum.beetle.vs.leaf.bup
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.bl - VS580.bl == 0  0.57356    0.04944  11.602  < 2e-16 ***
## VS620.bl - VS580.bl == 0  1.22449    0.04944  24.770  < 2e-16 ***
## VS640.bl - VS580.bl == 0  1.72490    0.04944  34.892  < 2e-16 ***
## VS660.bl - VS580.bl == 0  1.93477    0.04944  39.137  < 2e-16 ***
## VS620.bl - VS600.bl == 0  0.65093    0.04944  13.167  < 2e-16 ***
## VS640.bl - VS600.bl == 0  1.15135    0.04944  23.290  < 2e-16 ***
## VS660.bl - VS600.bl == 0  1.36121    0.04944  27.535  < 2e-16 ***
## VS640.bl - VS620.bl == 0  0.50041    0.04944  10.123  < 2e-16 ***
## VS660.bl - VS620.bl == 0  0.71028    0.04944  14.368  < 2e-16 ***
## VS660.bl - VS640.bl == 0  0.20986    0.04944   4.245 0.000218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Papilio xuthus

Anova(m.beetle.vs.leaf.pap) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 4046 4 0



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Papilio, 2)



Original model output

sum.beetle.vs.leaf.pap
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.bl - VS580.bl == 0  0.71260    0.04192  16.998  < 2e-16 ***
## VS620.bl - VS580.bl == 0  1.47367    0.04192  35.152  < 2e-16 ***
## VS640.bl - VS580.bl == 0  2.03704    0.04192  48.590  < 2e-16 ***
## VS660.bl - VS580.bl == 0  2.27321    0.04192  54.224  < 2e-16 ***
## VS620.bl - VS600.bl == 0  0.76107    0.04192  18.154  < 2e-16 ***
## VS640.bl - VS600.bl == 0  1.32444    0.04192  31.592  < 2e-16 ***
## VS660.bl - VS600.bl == 0  1.56060    0.04192  37.226  < 2e-16 ***
## VS640.bl - VS620.bl == 0  0.56337    0.04192  13.438  < 2e-16 ***
## VS660.bl - VS620.bl == 0  0.79954    0.04192  19.072  < 2e-16 ***
## VS660.bl - VS640.bl == 0  0.23617    0.04192   5.633 1.77e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Heliconius sp

Anova(m.beetle.vs.leaf.heli) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 3617 4 0



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Heliconius, 2)



Original model output

sum.beetle.vs.leaf.heli
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.bl - VS580.bl == 0  0.46256    0.02883  16.042  < 2e-16 ***
## VS620.bl - VS580.bl == 0  0.95329    0.02883  33.062  < 2e-16 ***
## VS640.bl - VS580.bl == 0  1.32086    0.02883  45.809  < 2e-16 ***
## VS660.bl - VS580.bl == 0  1.48159    0.02883  51.384  < 2e-16 ***
## VS620.bl - VS600.bl == 0  0.49073    0.02883  17.019  < 2e-16 ***
## VS640.bl - VS600.bl == 0  0.85829    0.02883  29.767  < 2e-16 ***
## VS660.bl - VS600.bl == 0  1.01903    0.02883  35.342  < 2e-16 ***
## VS640.bl - VS620.bl == 0  0.36756    0.02883  12.748  < 2e-16 ***
## VS660.bl - VS620.bl == 0  0.52830    0.02883  18.322  < 2e-16 ***
## VS660.bl - VS640.bl == 0  0.16074    0.02883   5.575 2.48e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Beetle vs Flower

Agrilus planipennis

Anova(m.beetle.vs.flower.bup) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 607.9 4 3.052e-130



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Agrilus, 3)



Original model output

sum.beetle.vs.flower.bup
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.bf - VS580.bf == 0  0.37988    0.07044   5.393 6.93e-07 ***
## VS620.bf - VS580.bf == 0  0.84419    0.07044  11.985  < 2e-16 ***
## VS640.bf - VS580.bf == 0  1.22628    0.07044  17.409  < 2e-16 ***
## VS660.bf - VS580.bf == 0  1.51287    0.07044  21.477  < 2e-16 ***
## VS620.bf - VS600.bf == 0  0.46431    0.07044   6.592 4.35e-10 ***
## VS640.bf - VS600.bf == 0  0.84640    0.07044  12.016  < 2e-16 ***
## VS660.bf - VS600.bf == 0  1.13299    0.07044  16.084  < 2e-16 ***
## VS640.bf - VS620.bf == 0  0.38209    0.07044   5.424 5.82e-07 ***
## VS660.bf - VS620.bf == 0  0.66868    0.07044   9.493  < 2e-16 ***
## VS660.bf - VS640.bf == 0  0.28659    0.07044   4.069 0.000473 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Papilio xuthus

Anova(m.beetle.vs.flower.pap) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 1264 4 2.33e-272



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Papilio, 3)



Original model output

sum.beetle.vs.flower.pap
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.bf - VS580.bf == 0  0.55397    0.06493   8.531  < 2e-16 ***
## VS620.bf - VS580.bf == 0  1.17015    0.06493  18.021  < 2e-16 ***
## VS640.bf - VS580.bf == 0  1.65504    0.06493  25.488  < 2e-16 ***
## VS660.bf - VS580.bf == 0  2.01940    0.06493  31.099  < 2e-16 ***
## VS620.bf - VS600.bf == 0  0.61618    0.06493   9.489  < 2e-16 ***
## VS640.bf - VS600.bf == 0  1.10107    0.06493  16.957  < 2e-16 ***
## VS660.bf - VS600.bf == 0  1.46543    0.06493  22.568  < 2e-16 ***
## VS640.bf - VS620.bf == 0  0.48489    0.06493   7.467 8.17e-13 ***
## VS660.bf - VS620.bf == 0  0.84925    0.06493  13.079  < 2e-16 ***
## VS660.bf - VS640.bf == 0  0.36436    0.06493   5.611 2.01e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Heliconius sp

Anova(m.beetle.vs.flower.heli) %>% pander()
Analysis of Deviance Table (Type II Wald chisquare tests)
  Chisq Df Pr(>Chisq)
vissys 1218 4 1.996e-262



Click the tabs to see the p-value summary plot or the original model output

Pair-wise p-values

get.pheatmap(heat.Heliconius, 3)



Original model output

sum.beetle.vs.flower.heli
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = dS ~ vissys + (1 | patch2) + (1 | patch1), data = inputdat, 
##     REML = F)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)    
## VS600.bf - VS580.bf == 0   0.3330     0.0405   8.221 2.22e-15 ***
## VS620.bf - VS580.bf == 0   0.7093     0.0405  17.514  < 2e-16 ***
## VS640.bf - VS580.bf == 0   1.0101     0.0405  24.940  < 2e-16 ***
## VS660.bf - VS580.bf == 0   1.2356     0.0405  30.510  < 2e-16 ***
## VS620.bf - VS600.bf == 0   0.3764     0.0405   9.293  < 2e-16 ***
## VS640.bf - VS600.bf == 0   0.6771     0.0405  16.719  < 2e-16 ***
## VS660.bf - VS600.bf == 0   0.9027     0.0405  22.289  < 2e-16 ***
## VS640.bf - VS620.bf == 0   0.3008     0.0405   7.426 1.12e-12 ***
## VS660.bf - VS620.bf == 0   0.5263     0.0405  12.996  < 2e-16 ***
## VS660.bf - VS640.bf == 0   0.2256     0.0405   5.570 2.55e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)



Plot the contrasts

  • Colours approximate human perception of either flower (left and right) or beetle (centre) colouration.
  • Lines of the same colors connect the same sample
#import the sample color code 
color.code <- read.csv("../data/color code list.csv", header=TRUE) 

# Set up a function to calculate means for the plot
get.plotmeans <- function(a, b, c){
  if (a == "aim1"){
    if(b == "bl"){
      plotmean.i <- c %>% 
        group_by(patch2, vissys) %>%
        summarize(mean.dS.sub = mean(dS)) %>% 
        ungroup() %>% 
        rename(beetleID = patch2) %>% 
        merge(beetle.colour.aim1[,c("species", "colour")], by.x = c("beetleID"), by.y = c("species")) %>% 
        distinct()
      return(plotmean.i)
    } else {
      plotmean.i <- c %>% 
        group_by(patch1, vissys) %>%
        summarize(mean.dS.sub = mean(dS)) %>% 
        ungroup() %>% 
        rename(flowerID = patch1) %>% 
        merge(flower.colour.aim1[,c("species", "colour")], by.x = c("flowerID"), by.y = c("species")) %>% 
        distinct()
      return(plotmean.i)
    }
  } else if (a == "aim2"){
    if (b == "bl"){
      plotmean.i <- c %>% 
        group_by(patch2, vissys) %>%
        summarize(mean.dS.sub = mean(dS)) %>% 
        ungroup() %>% 
        rename(beetleID = patch2) %>% 
        merge(beetle.colour.aim2[,c("species", "colour")], by.x = c("beetleID"), by.y = c("species")) %>% 
        distinct()
      return(plotmean.i)
    } else {
      plotmean.i <- c %>% 
        group_by(patch1, vissys) %>%
        summarize(mean.dS.sub = mean(dS)) %>% 
        ungroup() %>% 
        rename(flowerID = patch1) %>% 
        merge(flower.colour.aim2[,c("species", "colour")], by.x = c("flowerID"), by.y = c("species")) %>% 
        distinct()
      return(plotmean.i)
    }
  }
}

# What should be put in the function get.plotmeans()
# a = "aim1" or "aim2"
# b = comparison type: "fl", "bl", "bf"
# c = input data, e.g. allvis.fl_d65.aim1
#create a list of species name used in the spec data in oder to merge() with mean data set later
dataset.transpose <- gather(dataset[,2:131], 
                            key = "species", 
                            value = "reflectance", 
                            na.rm = FALSE,
                            convert = FALSE, factor_key = FALSE)
name.list <- unique(dataset.transpose$species) %>% 
  sort() %>% #order it alphabetically 
  data.frame() %>% 
  dplyr::rename(species = ".") #make it a data frame and name the column "species"

# create flower and beetle name list 
colourgrouplist <- c( "flower", "beetle")
tempcolordat <- vector("list", 4)

## for Aim 2
for(i in seq_along(colourgrouplist)){
  
  ouput.i <- name.list %>% 
    filter(str_detect(species, colourgrouplist[i])) %>% 
    cbind(color.code %>% 
            filter(str_detect(type, colourgrouplist[i])) %>% 
            arrange(name)) %>% 
    dplyr::select(-type) %>% 
    mutate(count = 5) %>% 
    uncount(count) 
  
  tempcolordat[[i+2]] <- ouput.i # here same in i+2 because 1-2 is for aim1
  
}
flower.colour.aim2 <- tempcolordat[[3]]
beetle.colour.aim2 <- tempcolordat[[4]]
allvis.list <- list(allvis.fl.Agrilus, allvis.bl.Agrilus, allvis.bf.Agrilus,
                    allvis.fl.Papilio, allvis.bl.Papilio, allvis.bf.Papilio,
                    allvis.fl.Heliconius, allvis.bl.Heliconius, allvis.bf.Heliconius)

tempcolordat <- vector("list", length(allvis.list))

for(i in seq_along(allvis.list)){
  
  temp.meandS.i <- allvis.list[[i]] %>% 
    group_by(vissys) %>%
    summarize(mean.dS = mean(dS))
  
  tempcolordat[[i]] <- temp.meandS.i
  
}
title.list <- list("Flower vs. Leaf", "Beetle vs. Leaf", "Beetle vs. Flower")

vis.comp.list <- list(fl <- c("VS580.fl", "VS600.fl", "VS620.fl", "VS640.fl", "VS660.fl"),
                      bl <- c("VS580.bl", "VS600.bl", "VS620.bl", "VS640.bl", "VS660.bl"),
                      bf <- c("VS580.bf", "VS600.bf", "VS620.bf", "VS640.bf", "VS660.bf"))
stats.list <- list(fl.Agrilus <- c("a","b","c","d","d"),
                   bl.Agrilus <- c("a","b","c","d","e"),
                   bf.Agrilus <- c("a","b","c","d","e"),
                   fl.Papilio <- c("a","b","c","d","d"),
                   bl.Papilio <- c("a","b","c","d","e"),
                   bf.Papilio <- c("a","b","c","d","e"),
                   fl.Heliconius <- c("a","b","c","d","d"),
                   bl.Heliconius <- c("a","b","c","d","e"),
                   bf.Heliconius <- c("a","b","c","d","e"))
backgorund.list <- c("grey80", NA) # choose NA for beetle.vs.leaf

get.contrastplot.sens <- function(dat, title, visomplist, compgroupnstats, background){
  ggplot(dat, aes(x = vissys, 
                y = mean.dS.sub,
                group = factor(dat[,1])))+
    geom_point(col = dat$colour,
               size = 1, alpha = 0.7) + 
    geom_line(col = dat$colour, 
              size = 0.5, alpha = 0.7)+
    xlab("Visual system") +     
    ylab("Chromatic contrast (JND)") +
    ylim(0, 21)+
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = backgorund.list[background], size = NA), 
          axis.title.x = element_text(size = 12),  
          axis.text.x  = element_text(size = 8, colour = "black", angle = 90),
          axis.title.y = element_text(size = 12, vjust = 1),
          axis.text.y = element_text(size = 8, colour = "black"),
          axis.line.x = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
          axis.line.y = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
          legend.justification = c(1,0), 
          legend.position = c(1,0.45),
          legend.key = element_blank(),
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 9))+
    scale_x_discrete(limits = vis.comp.list[[visomplist]],
                     labels = c("VS 580", "VS 600","VS 620","VS 640","VS 660"))+
    geom_point(data = dat, 
               aes(x = vis.comp.list[[visomplist]][1], 
                   y = tempcolordat[[compgroupnstats]]$mean.dS[1]),
               col = "black", size = 3)+
    geom_point(data = dat, 
               aes(x = vis.comp.list[[visomplist]][2], 
                   y = tempcolordat[[compgroupnstats]]$mean.dS[2]), 
               col = "black", size = 3)+
    geom_point(data = dat, 
               aes(x = vis.comp.list[[visomplist]][3], 
                   y = tempcolordat[[compgroupnstats]]$mean.dS[3]), 
               col = "black", size = 3)+
    geom_point(data = dat, 
               aes(x = vis.comp.list[[visomplist]][4], 
                   y = tempcolordat[[compgroupnstats]]$mean.dS[4]), 
               col = "black", size = 3)+
    geom_point(data = dat, 
               aes(x = vis.comp.list[[visomplist]][5], 
                   y = tempcolordat[[compgroupnstats]]$mean.dS[5]), 
               col = "black", size = 3)+
    ggtitle(title.list[[title]])+
    annotate("text", x = vis.comp.list[[visomplist]], y = 21, 
             label = stats.list[[compgroupnstats]])
    
}

Agrilus planipennis

# create a data set for means
mean.fl.Agrilus <- get.plotmeans("aim2", "fl", allvis.fl.Agrilus)
mean.bl.Agrilus <- get.plotmeans("aim2", "bl", allvis.bl.Agrilus)
mean.bf.Agrilus <- get.plotmeans("aim2", "bf", allvis.bf.Agrilus)

# plot individual panel
fl.Agrilus <- get.contrastplot.sens(mean.fl.Agrilus, 1, 1, 1, 1)
bl.Agrilus <- get.contrastplot.sens(mean.bl.Agrilus, 2, 2, 2, 2)
bf.Agrilus <- get.contrastplot.sens(mean.bf.Agrilus, 3, 3, 3, 1)

# bind the 3 panels
figure.bup <- ggarrange(fl.Agrilus,
                             bl.Agrilus,
                             bf.Agrilus, 
                             ncol = 3, nrow = 1)
figure.bup
Figure caption: Comparison of chromatic contrast between visual systems with the LWS photoreceptor peaking at different wavelengths in the jewel beetle, _Agrilus planipenni_. Black dots show the means of the representative contrasts in visual systems. Each coloured dot represents the average contrast of each flower pattern to all leaves (left panel), each beetle colour to all leaves (middle panel) or each flower colour to all beetle colours (right panel). Colours of the dots correspond to the human-visible colour of the flower (left and right panels) or beetle (middle panel) with the lines of the same colour connecting the results between different visual systems. This is for graphical representation only; statistical tests are based on all pairwise combinations of spectra and not averages. Letters on the top of each panel show the significant difference in contrast between visual systems. Dotted horizontal line indicates JND=1. Two contrasts > 12 JND are from flowers that have high UV - blue chroma compared to beetles and leaves.

Figure caption: Comparison of chromatic contrast between visual systems with the LWS photoreceptor peaking at different wavelengths in the jewel beetle, Agrilus planipenni. Black dots show the means of the representative contrasts in visual systems. Each coloured dot represents the average contrast of each flower pattern to all leaves (left panel), each beetle colour to all leaves (middle panel) or each flower colour to all beetle colours (right panel). Colours of the dots correspond to the human-visible colour of the flower (left and right panels) or beetle (middle panel) with the lines of the same colour connecting the results between different visual systems. This is for graphical representation only; statistical tests are based on all pairwise combinations of spectra and not averages. Letters on the top of each panel show the significant difference in contrast between visual systems. Dotted horizontal line indicates JND=1. Two contrasts > 12 JND are from flowers that have high UV - blue chroma compared to beetles and leaves.



Papilio xuthus

# create a data set for means
mean.fl.Papilio <- get.plotmeans("aim2", "fl", allvis.fl.Papilio)
mean.bl.Papilio <- get.plotmeans("aim2", "bl", allvis.bl.Papilio)
mean.bf.Papilio <- get.plotmeans("aim2", "bf", allvis.bf.Papilio)

# plot individual panel
fl.Papilio <- get.contrastplot.sens(mean.fl.Papilio, 1, 1, 4, 1)
bl.Papilio <- get.contrastplot.sens(mean.bl.Papilio, 2, 2, 5, 2)
bf.Papilio <- get.contrastplot.sens(mean.bf.Papilio, 3, 3, 6, 1)

# bind the 3 panels
figure.pap <- ggarrange(fl.Papilio,
                        bl.Papilio,
                        bf.Papilio, 
                        ncol = 3, nrow = 1)
figure.pap
Figure caption: Comparison of chromatic contrast between visual systems with the LWS photoreceptor peaking at different wavelengths in the butterfly, _Papilio xuthus_. Black dots show the means of the representative contrasts in visual systems. Each coloured dot represents the average contrast of each flower pattern to all leaves (left panel), each beetle colour to all leaves (middle panel) or each flower colour to all beetle colours (right panel). Colours of the dots correspond to the human-visible colour of the flower (left and right panels) or beetle (middle panel) with the lines of the same colour connecting the results between different visual systems. This is for graphical representation only; statistical tests are based on all pairwise combinations of spectra and not averages. Letters on the top of each panel show the significant difference in contrast between visual systems. Dotted horizontal line indicates JND=1.

Figure caption: Comparison of chromatic contrast between visual systems with the LWS photoreceptor peaking at different wavelengths in the butterfly, Papilio xuthus. Black dots show the means of the representative contrasts in visual systems. Each coloured dot represents the average contrast of each flower pattern to all leaves (left panel), each beetle colour to all leaves (middle panel) or each flower colour to all beetle colours (right panel). Colours of the dots correspond to the human-visible colour of the flower (left and right panels) or beetle (middle panel) with the lines of the same colour connecting the results between different visual systems. This is for graphical representation only; statistical tests are based on all pairwise combinations of spectra and not averages. Letters on the top of each panel show the significant difference in contrast between visual systems. Dotted horizontal line indicates JND=1.



Heliconius sp

# create a data set for means
mean.fl.Heliconius <- get.plotmeans("aim2", "fl", allvis.fl.Heliconius)
mean.bl.Heliconius <- get.plotmeans("aim2", "bl", allvis.bl.Heliconius)
mean.bf.Heliconius <- get.plotmeans("aim2", "bf", allvis.bf.Heliconius)

# plot individual panel
fl.Heliconius <- get.contrastplot.sens(mean.fl.Heliconius, 1, 1, 7, 1)
bl.Heliconius <- get.contrastplot.sens(mean.bl.Heliconius, 2, 2, 8, 2)
bf.Heliconius <- get.contrastplot.sens(mean.bf.Heliconius, 3, 3, 9, 1)

# bind the 3 panels
figure.heli <- ggarrange(fl.Heliconius,
                             bl.Heliconius,
                             bf.Heliconius, 
                             ncol = 3, nrow = 1)
figure.heli
Figure caption: Comparison of chromatic contrast between visual systems with the LWS photoreceptor peaking at different wavelengths in the butterfly, _Heliconius_ sp. Black dots show the means of the representative contrasts in visual systems. Each coloured dot represents the average contrast of each flower pattern to all leaves (left panel), each beetle colour to all leaves (middle panel) or each flower colour to all beetle colours (right panel). Colours of the dots correspond to the human-visible colour of the flower (left and right panels) or beetle (middle panel) with the lines of the same colour connecting the results between different visual systems. This is for graphical representation only; statistical tests are based on all pairwise combinations of spectra and not averages. Letters on the top of each panel show the significant difference in contrast between visual systems. Dotted horizontal line indicates JND=1.

Figure caption: Comparison of chromatic contrast between visual systems with the LWS photoreceptor peaking at different wavelengths in the butterfly, Heliconius sp. Black dots show the means of the representative contrasts in visual systems. Each coloured dot represents the average contrast of each flower pattern to all leaves (left panel), each beetle colour to all leaves (middle panel) or each flower colour to all beetle colours (right panel). Colours of the dots correspond to the human-visible colour of the flower (left and right panels) or beetle (middle panel) with the lines of the same colour connecting the results between different visual systems. This is for graphical representation only; statistical tests are based on all pairwise combinations of spectra and not averages. Letters on the top of each panel show the significant difference in contrast between visual systems. Dotted horizontal line indicates JND=1.



Summary

The results were qualitatively the same for different photoreceptor ratios among three comparison groups – contrast increased for LWS photoreceptors with peak sensitivity from 580 to at least 640 nm. This indicates that our results were not biased by photoreceptor ratios in our hypothetical visual systems.