Multivariate analysis of eye colour in Flatback (Natator depressus) hatchlings in relation to nest temperature
Ella-Monique Mason
2025-10-23
Libraries
library(tidyverse)
library(vegan)
library(viridisLite)Import + process eye data
eye <- read.csv("eyeRBG.csv", stringsAsFactors=TRUE)
summary(eye)## Nest_id Colour_channel MEAN MODE STDV
## N6H7 : 4 blue :90 Min. :26.03 Min. :19.00 Min. : 3.86
## R4H5 : 4 green:90 1st Qu.:44.77 1st Qu.:36.00 1st Qu.:12.26
## N1H1 : 3 red :90 Median :49.92 Median :42.50 Median :15.84
## N1H10 : 3 Mean :50.68 Mean :43.26 Mean :16.59
## N1H2 : 3 3rd Qu.:56.62 3rd Qu.:50.00 3rd Qu.:20.27
## N1H3 : 3 Max. :72.58 Max. :74.00 Max. :31.71
## (Other):250head(eye)## Nest_id Colour_channel MEAN MODE STDV
## 1 N1H1 red 31.97 28 8.01
## 2 N1H1 green 39.63 37 9.63
## 3 N1H1 blue 41.04 37 10.22
## 4 N1H2 red 44.17 41 12.09
## 5 N1H2 green 53.19 52 12.54
## 6 N1H2 blue 49.14 47 12.63eye_rgb <- eye %>%
group_by(Nest_id, Colour_channel) %>%
summarise(mean_value = mean(MEAN, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = Colour_channel, values_from = mean_value) %>%
rename(R = red, G = green, B = blue) %>%
mutate(Nest_code = gsub("H.*", "", Nest_id))
summary(eye_rgb)## Nest_id B G R
## N1H1 : 1 Min. :29.05 Min. :31.83 Min. :26.03
## N1H10 : 1 1st Qu.:46.95 1st Qu.:48.80 1st Qu.:40.78
## N1H2 : 1 Median :51.38 Median :53.91 Median :45.59
## N1H3 : 1 Mean :51.91 Mean :54.39 Mean :45.69
## N1H4 : 1 3rd Qu.:56.98 3rd Qu.:60.94 3rd Qu.:50.82
## N1H5 : 1 Max. :69.50 Max. :72.58 Max. :63.07
## (Other):84 NA's :2
## Nest_code
## Length:90
## Class :character
## Mode :character
##
##
##
## head(eye_rgb)## # A tibble: 6 × 5
## Nest_id B G R Nest_code
## <fct> <dbl> <dbl> <dbl> <chr>
## 1 N1H1 41.0 39.6 32.0 N1
## 2 N1H10 54.0 57.4 47.8 N1
## 3 N1H2 49.1 53.2 44.2 N1
## 4 N1H3 33.1 31.8 26.0 N1
## 5 N1H4 51.1 49.4 40.8 N1
## 6 N1H5 61.0 62.7 51.5 N1Add temp
temp <- read_csv("fb_temp.csv")
bottom_vars <- grep("_B$", names(temp), value = TRUE)
nest_medians <- sapply(bottom_vars, function(nest) median(temp[[nest]], na.rm = TRUE))
temp_df <- data.frame(
Nest_code = gsub("_B", "", names(nest_medians)),
Bottom_temp = round(nest_medians, 2)
)
eye_rgb_temp <- left_join(eye_rgb, temp_df, by = "Nest_code") %>%
drop_na()
summary(eye_rgb_temp)## Nest_id B G R
## N1H1 : 1 Min. :29.05 Min. :31.83 Min. :26.03
## N1H10 : 1 1st Qu.:46.95 1st Qu.:48.52 1st Qu.:40.72
## N1H2 : 1 Median :50.86 Median :53.64 Median :45.53
## N1H3 : 1 Mean :51.65 Mean :54.02 Mean :45.55
## N1H4 : 1 3rd Qu.:56.77 3rd Qu.:60.66 3rd Qu.:50.46
## N1H5 : 1 Max. :69.50 Max. :72.58 Max. :63.07
## (Other):80
## Nest_code Bottom_temp
## Length:86 Min. :29.46
## Class :character 1st Qu.:29.79
## Mode :character Median :30.02
## Mean :30.10
## 3rd Qu.:30.67
## Max. :30.82
## head(eye_rgb_temp)## # A tibble: 6 × 6
## Nest_id B G R Nest_code Bottom_temp
## <fct> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 N1H1 41.0 39.6 32.0 N1 30.8
## 2 N1H10 54.0 57.4 47.8 N1 30.8
## 3 N1H2 49.1 53.2 44.2 N1 30.8
## 4 N1H3 33.1 31.8 26.0 N1 30.8
## 5 N1H4 51.1 49.4 40.8 N1 30.8
## 6 N1H5 61.0 62.7 51.5 N1 30.8hatchling level proportions
eye_rgb_prop <- eye_rgb_temp %>%
mutate(
total_rgb = R + G + B,
prop_R = R / total_rgb,
prop_G = G / total_rgb,
prop_B = B / total_rgb
)
summary(eye_rgb_prop)## Nest_id B G R
## N1H1 : 1 Min. :29.05 Min. :31.83 Min. :26.03
## N1H10 : 1 1st Qu.:46.95 1st Qu.:48.52 1st Qu.:40.72
## N1H2 : 1 Median :50.86 Median :53.64 Median :45.53
## N1H3 : 1 Mean :51.65 Mean :54.02 Mean :45.55
## N1H4 : 1 3rd Qu.:56.77 3rd Qu.:60.66 3rd Qu.:50.46
## N1H5 : 1 Max. :69.50 Max. :72.58 Max. :63.07
## (Other):80
## Nest_code Bottom_temp total_rgb prop_R
## Length:86 Min. :29.46 Min. : 90.95 Min. :0.2679
## Class :character 1st Qu.:29.79 1st Qu.:136.63 1st Qu.:0.2920
## Mode :character Median :30.02 Median :149.93 Median :0.3036
## Mean :30.10 Mean :151.22 Mean :0.3014
## 3rd Qu.:30.67 3rd Qu.:168.22 3rd Qu.:0.3111
## Max. :30.82 Max. :203.22 Max. :0.3292
##
## prop_G prop_B
## Min. :0.3414 Min. :0.3070
## 1st Qu.:0.3528 1st Qu.:0.3327
## Median :0.3573 Median :0.3377
## Mean :0.3570 Mean :0.3415
## 3rd Qu.:0.3609 3rd Qu.:0.3499
## Max. :0.3724 Max. :0.3813
## head(eye_rgb_prop)## # A tibble: 6 × 10
## Nest_id B G R Nest_code Bottom_temp total_rgb prop_R prop_G prop_B
## <fct> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 N1H1 41.0 39.6 32.0 N1 30.8 113. 0.284 0.352 0.364
## 2 N1H10 54.0 57.4 47.8 N1 30.8 159. 0.300 0.360 0.339
## 3 N1H2 49.1 53.2 44.2 N1 30.8 146. 0.302 0.363 0.335
## 4 N1H3 33.1 31.8 26.0 N1 30.8 91.0 0.286 0.350 0.364
## 5 N1H4 51.1 49.4 40.8 N1 30.8 141. 0.289 0.349 0.362
## 6 N1H5 61.0 62.7 51.5 N1 30.8 175. 0.294 0.358 0.348Unconstrained ordination PCA
With a vector for bottom temp
# set colour gradient for e.g. propB
b_vals <- eye_rgb_prop$prop_B
cols_B <- colorRampPalette(c("#deebf7", "#3182bd"))(100)[cut(b_vals, 100, include.lowest = TRUE)]
# matrix of proportional RGB values
rgb_mat <- eye_rgb_prop %>%
select(prop_R, prop_G, prop_B) %>%
as.matrix()
# create PCA - also check Env3/511 multivariate stats notes
pca_rgb <- rda(rgb_mat, scale = TRUE)
# extract sample ("sites") scores (PCA coordinates)
sc <- scores(pca_rgb, display = "sites", scaling = 2)
# draw the plot
plot(pca_rgb, type = "n", scaling = 2) # empty axes
points(sc[,1], sc[,2], pch = 19, cex = 1, col = cols_B)
# legend bar
brks <- pretty(range(b_vals, na.rm = TRUE), 5)
legend_cols <- colorRampPalette(c("#deebf7", "#3182bd"))(length(brks)-1)
par(xpd = NA)
legend("topright", inset = 0.02, title = "prop_B",
legend = paste(head(brks, -1), "–", tail(brks, -1)),
fill = legend_cols, border = NA, cex = 0.8)
par(xpd = FALSE)
# add temp gradient
#with(eye_rgb_prop, ordisurf(pca_rgb, Bottom_temp, add = TRUE, col = "red"))
# add temp vector
fit <- envfit(pca_rgb, eye_rgb_prop$Bottom_temp)
plot(fit, col = "red", add = TRUE) # arrow points toward increasing temperature
# add RGB vectors
sp <- scores(pca_rgb, display = "species", scaling = 2)
mul <- ordiArrowMul(sp, fill = 0.99) # scales arrows
# draw arrows
arrows(0, 0, sp[,1] * mul, sp[,2] * mul, length = 0.06, angle = 25, lwd = 1.1)
# put labels slightly inside each arrow tip
lab_pos <- sp * (mul * 1.05)
text(lab_pos[,1], lab_pos[,2], labels = rownames(sp), font = 1)
RDA Constrained ordination
RDA - How much does bottom temp explain the RGB proportions
Bottom temp explained 3.7% of the variance in the hatchling eye RGB proportions (permutation test, P=0.056)
# constrain by temperature
rda_rgb <- rda(rgb_mat ~ Bottom_temp, data = eye_rgb_prop)
RsquareAdj(rda_rgb)## $r.squared
## [1] 0.03727381
##
## $adj.r.squared
## [1] 0.02581278anova(rda_rgb, permutations = 999) # overall## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = rgb_mat ~ Bottom_temp, data = eye_rgb_prop)
## Df Variance F Pr(>F)
## Model 1 0.00001511 3.2522 0.057 .
## Residual 84 0.00039021
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(rda_rgb, by = "axis", permutations = 999) # tests RDA1 specifically## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = rgb_mat ~ Bottom_temp, data = eye_rgb_prop)
## Df Variance F Pr(>F)
## RDA1 1 0.00001511 3.2522 0.055 .
## Residual 84 0.00039021
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(rda_rgb)##
## Call:
## rda(formula = rgb_mat ~ Bottom_temp, data = eye_rgb_prop)
##
## Partitioning of variance:
## Inertia Proportion
## Total 4.053e-04 1.00000
## Constrained 1.511e-05 0.03727
## Unconstrained 3.902e-04 0.96273
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 PC1 PC2
## Eigenvalue 1.511e-05 0.0003347 5.549e-05
## Proportion Explained 3.727e-02 0.8258251 1.369e-01
## Cumulative Proportion 3.727e-02 0.8630989 1.000e+00
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1
## Eigenvalue 1.511e-05
## Proportion Explained 1.000e+00
## Cumulative Proportion 1.000e+00plot(rda_rgb, display = "sites")