Juvenile Pacific cod were tagged, acclimated to laboratory conditions, then gradually acclimated to four experimental temperatures (0, 5, 9, 16). After 6 weeks in treatments, fish were sacrificed, liver and fin tissues were collected for RNASeq and lcWGS, respectively, and measurements taken of fish length, wet weight, and liver weight.

On 11/21/22 160 tagged fish entered experimental tanks for acclimation, 40 per treatment (n=10 fish / tank, 4 tanks/treatment). On 12/28/22 temperatures began to slowly increase until all target experimental temperatures (0C, 5C, 9C, 16C) were reached on 1/8/23. Of the 40 fish per treatment, all survived the 0C, 5C, and 9C treatments. Four fish died in the 16degC treatment on 1/22/23 (tank 1), 1/24/23 (tank 2), 1/29/23 (tank 4), and 2/1/23 (tank 2). All survivors were euthanized/sampled on 2/8 & 2/9. Standard length and wet weight were 1) collected before tank acclimation, 2) before experimental temperature, and 3) at treatment termination so that growth and condition could be assessed prior to and during experimental treatments. Liver wet weight and tissues for sequencing were collected at treatment termination.

In this notebook I look at effects of temperature on distributions of growth rate (length, weight) body condition index (wet weight / length), liver condition (i.e. hepato-somatic index, liver wet weight / whole body wet weight), and lipids.

Load libraries and source scripts

list.of.packages <- c("tidyverse", "readxl", "janitor", "purrr", "ggpubr", "googlesheets4", "plotly", "lubridate", "scales", "knitr", "cowplot")

# Load all required libraries
all.packages <- c(list.of.packages)
lapply(all.packages, FUN = function(X) {
  do.call("require", list(X))
})

`%!in%` = Negate(`%in%`)

Load data and calculate per-fish metrics

Specific growth rates (SGR) is calculated for each fish based on wet weights and standard length using the below equation, where M1 and M2 are measurements (standard length, wet weight) taken before (t1) and after (t2) the acclimation period to assess baseline growth rates, and before and after experimental temperature treatments to assess effects of temperature (Table 1).

SGR(standard length, wet weight) = 100 * (e^g-1) where g = (ln(M2) - ln(M1)) / t2-t1

Fulton’s condition index (Kwet) was approximated as K(wet) = 100 * wet total weight (g) / (standard length (mm)/10)^3

The hepato-somatic index (HSI) was calculated as HSI =100 x wet liver weight (g)wet total weight (g).

load(file = "../rnaseq/aligned/counts.ts")

phenotypes <- read_excel("../data/Pcod Temp Growth experiment 2022-23 DATA.xlsx", sheet = "AllData") %>% clean_names() %>%
  mutate(acc_t1=as.Date("2022-11-21", format="%Y-%m-%d"), exp_t1=as.Date("2022-12-27", format= "%Y-%m-%d"), dissection_date=as.Date(dissection_date, format= "%Y-%m-%d")) %>%
  mutate_at(c("tank", "temperature", "microchip_id", "genetic_sampling_count"), factor) %>%
  dplyr::rename("sl_final"="sl_mm", "wwt_final"="whole_body_ww_g") %>% 
  mutate(SGR.sl.accl=round(100*(exp((log(sl_12272022)-log(sl_11212022))/(as.numeric(exp_t1-acc_t1)))-1), digits=5),
         SGR.sl.trt=round(100*(exp((log(sl_final)-log(sl_12272022))/(as.numeric(dissection_date-exp_t1)))-1), digits=5),
         SGR.ww.accl=round(100*(exp((log(wwt_12272022)-log(wwt_11212022))/(as.numeric(exp_t1-acc_t1)))-1), digits=5),
         SGR.ww.trt=round(100*(exp((log(wwt_final)-log(wwt_12272022))/(as.numeric(dissection_date-exp_t1)))-1), digits=5),
         Kwet=100*(wwt_final/(sl_final/10)^3),
         hsi=100*total_liver_ww_mg/(wwt_final),
         mort=as.factor(case_when(is.na(mort_date) ~ "No", TRUE ~ "Yes")),
         rna=case_when(genetic_sampling_count %in% gsub("sample_", "", rownames(counts.ts)) ~ "Yes", TRUE ~ "No"),
         measure_date=case_when(is.na(mort_date) ~ ymd(dissection_date), TRUE ~ mort_date)) %>%
  mutate(days_growth = as.numeric(difftime(ymd(measure_date), ymd("2022-12-27"), units = "days")))

phenotypes %>% head()
## # A tibble: 6 × 34
##   microchip_id sl_11212022 wwt_11212022 tank  temperature sl_12272022
##   <fct>              <dbl>        <dbl> <fct> <fct>             <dbl>
## 1 9397                  86         7.38 12    0                    93
## 2 9467                  92         7.98 12    0                    97
## 3 5059                 102        11.6  12    0                   113
## 4 9461                  82         6.15 12    0                    88
## 5 9560                  95         9.59 12    0                   101
## 6 7386                  95         9.23 12    0                   108
## # ℹ 28 more variables: wwt_12272022 <dbl>, mort_date <dttm>,
## #   dissection_date <date>, sl_final <dbl>, wwt_final <dbl>,
## #   total_liver_ww_mg <dbl>, liverfor_lipids_ww_mg <dbl>,
## #   muscle_w_wfor_lipids_mg <dbl>, genetic_sampling_count <fct>,
## #   dissection_comments <chr>, molecular_notes <chr>, fin <chr>, liver <chr>,
## #   blood <chr>, spleen <chr>, gill <chr>, acc_t1 <date>, exp_t1 <date>,
## #   SGR.sl.accl <dbl>, SGR.sl.trt <dbl>, SGR.ww.accl <dbl>, SGR.ww.trt <dbl>, …
# Wet weight over time
phenotypes %>%
  group_by(temperature) %>%
  dplyr::summarise(sl_mean.1=mean(sl_11212022),sl_sd.1=sd(sl_11212022), 
            sl_mean.2=mean(sl_12272022),sl_sd.2=sd(sl_12272022),
            sl_mean.3=mean(sl_final),sl_sd.3=sd(sl_final),
            
            wwt_mean.1=mean(wwt_11212022),wwt_sd.1=sd(wwt_11212022), 
            wwt_mean.2=mean(wwt_12272022),wwt_sd.2=sd(wwt_12272022),
            wwt_mean.3=mean(wwt_final),wwt_sd.3=sd(wwt_final)) %>%
  pivot_longer(cols = -temperature) %>%
  separate(name, sep = "\\.", into = c("metric", "time")) %>%
  mutate_at(c("metric"), factor) %>% mutate(time=as.numeric(time)) %>%
  filter(metric=="wwt_mean") %>%
  ggplot() + geom_line(aes(x=time, y=value, color=temperature), cex=1.5) + theme_minimal() +
  scale_color_manual(name="Temperature", 
                       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) + 
  ggtitle("Wet weight over time") + xlab(NULL) + ylab("Wet weight") +
  scale_x_continuous(breaks = c(1, 2, 3), labels = c("Pre-Acclimation", "Pre-Treatment", "Treatment\ntermination"))

### Standard length by temperature over time - mean and SD
phenotypes %>%
  group_by(temperature) %>%
  dplyr::summarise(sl_mean.1=mean(sl_11212022),sl_sd.1=sd(sl_11212022), 
            sl_mean.2=mean(sl_12272022),sl_sd.2=sd(sl_12272022),
            sl_mean.3=mean(sl_final),sl_sd.3=sd(sl_final),
            
            wwt_mean.1=mean(wwt_11212022),wwt_sd.1=sd(wwt_11212022), 
            wwt_mean.2=mean(wwt_12272022),wwt_sd.2=sd(wwt_12272022),
            wwt_mean.3=mean(wwt_final),wwt_sd.3=sd(wwt_final)) %>%
  pivot_longer(cols = -temperature) %>%
  separate(name, sep = "\\.", into = c("metric", "time")) %>%
  mutate_at(c("metric"), factor) %>% mutate(time=as.numeric(time)) %>%
  filter(metric=="sl_mean") %>%
  ggplot() + geom_line(aes(x=time, y=value, color=temperature), cex=1.5) + theme_minimal() +
  scale_color_manual(name="Temperature", 
                       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) + 
  ggtitle("Standard length over time") + xlab(NULL) + ylab("Standard length") +
  scale_x_continuous(breaks = c(1, 2, 3), labels = c("Pre-Acclimation", "Pre-Treatment", "Treatment\ntermination"))

Specific growth while acclimating, Standard Length

Statistics - did it differ among treatments? No. 

phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(SGR.sl.accl), sd=sd(SGR.sl.accl))
## # A tibble: 4 × 3
##   temperature  mean     sd
##   <fct>       <dbl>  <dbl>
## 1 0           0.256 0.0815
## 2 5           0.255 0.0956
## 3 9           0.240 0.0664
## 4 16          0.282 0.0730
summary(aov(SGR.sl.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq  Mean Sq F value Pr(>F)
## temperature   3 0.0341 0.011364   1.773  0.155
## Residuals   152 0.9741 0.006409
TukeyHSD(aov(SGR.sl.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = SGR.sl.accl ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##             diff          lwr        upr     p adj
## 5-0  -0.00075475 -0.047254817 0.04574532 0.9999730
## 9-0  -0.01565400 -0.062154067 0.03084607 0.8180483
## 16-0  0.02624806 -0.021526222 0.07402233 0.4843983
## 9-5  -0.01489925 -0.061399317 0.03160082 0.8390225
## 16-5  0.02700281 -0.020771472 0.07477708 0.4592131
## 16-9  0.04190206 -0.005872222 0.08967633 0.1076556
phenotypes %>% 
    ggplot(aes(y=SGR.sl.accl, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin(trim = F) + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Specific growth rate while acclimating, standard length (~ %/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_x_discrete(drop=T) + #Do drop empty factors for temperature contrasts
  ylab("mm / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) + ylab(NULL)

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=SGR.sl.accl/days_growth, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Growth rate while acclimating, standard length (mm/day)") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) + 
#   xlab(NULL) + ylab(NULL)

Specific growth in treatments, Standard Length

Did SGR(sl) differ among treatments? Yes! Among all pairwise comparisons.

Note: should see if this data fits the model/estimates using this formula from Laurel et al. 2015, https://doi.org/10.1007/s00300-015-1761-5 G indiv = y o + aT + bT 2 + cT 3

summary(aov(SGR.sl.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## temperature   3 1.1303  0.3768   90.09 <2e-16 ***
## Residuals   152 0.6357  0.0042                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(SGR.sl.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = SGR.sl.trt ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##             diff         lwr          upr     p adj
## 5-0   0.11158550  0.07402061  0.149150392 0.0000000
## 9-0   0.22290525  0.18534036  0.260470142 0.0000000
## 16-0  0.18173008  0.14313583  0.220324341 0.0000000
## 9-5   0.11131975  0.07375486  0.148884642 0.0000000
## 16-5  0.07014458  0.03155033  0.108738841 0.0000311
## 16-9 -0.04117517 -0.07976942 -0.002580909 0.0315640
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(SGR.sl.trt), sd=sd(SGR.sl.trt))
## # A tibble: 4 × 3
##   temperature   mean     sd
##   <fct>        <dbl>  <dbl>
## 1 0           0.0934 0.0396
## 2 5           0.205  0.0588
## 3 9           0.316  0.0695
## 4 16          0.275  0.0848
phenotypes %>% 
    ggplot(aes(y=SGR.sl.trt, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Specific growth rate in treatments, standard length (~%day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  ylab("mm / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) + ylab(NULL)

SGR.sl.model <- lm(SGR.sl.trt ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes)
summary(SGR.sl.model)
## 
## Call:
## lm(formula = SGR.sl.trt ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.252484 -0.031726  0.007516  0.046239  0.213046 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                      0.216789   0.005785  37.474
## poly(as.numeric(as.character(temperature)), 3)1  0.760624   0.073177  10.394
## poly(as.numeric(as.character(temperature)), 3)2 -0.648036   0.073177  -8.856
## poly(as.numeric(as.character(temperature)), 3)3 -0.251380   0.073177  -3.435
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1  < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)2 1.74e-15 ***
## poly(as.numeric(as.character(temperature)), 3)3 0.000759 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07318 on 156 degrees of freedom
## Multiple R-squared:  0.5597, Adjusted R-squared:  0.5512 
## F-statistic: 66.09 on 3 and 156 DF,  p-value: < 2.2e-16
SGR.sl.predict <- data.frame(temperature=seq(0,16,length.out=100))
SGR.sl.predict$predicted <- predict(SGR.sl.model, SGR.sl.predict)

phenotypes %>% 
    ggplot(aes(y=SGR.sl.trt, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Specific growth rate in treatments, standard length (~%day)") +
  geom_line(data=SGR.sl.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("% / day")

round(SGR.sl.model$coefficients, digits=3) %>% as.vector() %>% as.data.frame(row.names = c("y","a","b","c")) %>% set_names(NULL) %>% kable(caption="Model coefficients")
Model coefficients
y 0.217
a 0.761
b -0.648
c -0.251 
# What is the maximum standard length growth rate and associated predicted temperature as per this model?  0.35%/day at 12C.
SGR.sl.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted growth rate & associated temperature, standard length, %/day")
Maximum predicted growth rate & associated temperature, standard length, %/day
temperature predicted
11.96 0.35

Specific growth while acclimating, wet weight

No difference among temperatures

summary(aov(SGR.ww.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)
## temperature   3  0.055 0.01821   0.859  0.464
## Residuals   152  3.221 0.02119
TukeyHSD(aov(SGR.ww.accl ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = SGR.ww.accl ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##              diff         lwr        upr     p adj
## 5-0   0.007924250 -0.07663468 0.09248318 0.9949093
## 9-0  -0.039128500 -0.12368743 0.04543043 0.6265002
## 16-0  0.001401917 -0.08547412 0.08827795 0.9999734
## 9-5  -0.047052750 -0.13161168 0.03750618 0.4731401
## 16-5 -0.006522333 -0.09339837 0.08035370 0.9973613
## 16-9  0.040530417 -0.04634562 0.12740645 0.6202693
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(SGR.ww.accl), sd=sd(SGR.ww.accl))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0           0.700 0.168
## 2 5           0.708 0.151
## 3 9           0.661 0.134
## 4 16          0.702 0.122
phenotypes %>% 
    ggplot(aes(y=SGR.ww.accl, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Growth rate while acclimating, wet weight (g/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  ylab("g / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL)

#### Specific growth in treatments, Wet Weight

Did SGR(wet weight) differ among treatments? Yes! Among all pairwise comparisons.

Note: should see if this data fits the model/estimates using this formula from Laurel et al. 2015, https://doi.org/10.1007/s00300-015-1761-5 G indiv = y o + aT + bT 2 + cT 3

summary(aov(SGR.ww.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## temperature   3 10.121   3.374   116.4 <2e-16 ***
## Residuals   152  4.405   0.029                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(SGR.ww.trt ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = SGR.ww.trt ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##            diff        lwr         upr     p adj
## 5-0   0.3083317  0.2094470  0.40721653 0.0000000
## 9-0   0.6631442  0.5642595  0.76202903 0.0000000
## 16-0  0.5414311  0.4398367  0.64302557 0.0000000
## 9-5   0.3548125  0.2559277  0.45369728 0.0000000
## 16-5  0.2330994  0.1315049  0.33469382 0.0000001
## 16-9 -0.1217131 -0.2233076 -0.02011868 0.0117878
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(SGR.ww.trt), sd=sd(SGR.ww.trt))
## # A tibble: 4 × 3
##   temperature   mean    sd
##   <fct>        <dbl> <dbl>
## 1 0           0.0109 0.108
## 2 5           0.319  0.168
## 3 9           0.674  0.190
## 4 16          0.552  0.202
phenotypes %>% 
    ggplot(aes(y=SGR.ww.trt, x=temperature, fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_point(position = position_jitterdodge(jitter.width = 1)) + theme_minimal() +
    ggtitle("Specific growth rate in treatments, wet weight (~%/day)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  ylab("mm / day") +
  scale_fill_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) + ylab(NULL)

SGR.ww.model <- lm(SGR.ww.trt ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes)
summary(SGR.ww.model)
## 
## Call:
## lm(formula = SGR.ww.trt ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.88482 -0.09333  0.00724  0.11338  0.49895 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                      0.36926    0.01652  22.358
## poly(as.numeric(as.character(temperature)), 3)1  2.24342    0.20891  10.739
## poly(as.numeric(as.character(temperature)), 3)2 -1.89754    0.20891  -9.083
## poly(as.numeric(as.character(temperature)), 3)3 -0.86586    0.20891  -4.145
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1  < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)2 4.46e-16 ***
## poly(as.numeric(as.character(temperature)), 3)3 5.57e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2089 on 156 degrees of freedom
## Multiple R-squared:  0.5795, Adjusted R-squared:  0.5714 
## F-statistic: 71.67 on 3 and 156 DF,  p-value: < 2.2e-16
SGR.ww.predict <- data.frame(temperature=seq(0,16,length.out=100))
SGR.ww.predict$predicted <- predict(SGR.ww.model, SGR.ww.predict)

phenotypes %>% 
    ggplot(aes(y=SGR.ww.trt, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Specific growth rate in treatments, wet weight (~%/day)") +
  geom_line(data=SGR.ww.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("% / day")

round(SGR.ww.model$coefficients, digits=3) %>% as.vector() %>% as.data.frame(row.names = c("y","a","b","c")) %>% set_names(NULL) %>% kable(caption="Model coefficients")
Model coefficients
y 0.369
a 2.243
b -1.898
c -0.866 
# What is the maximum standard length growth rate and associated predicted temperature as per this model?  0.35%/day at 12C.
SGR.ww.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted growth rate & associated temperature, wet weight, %/day")
Maximum predicted growth rate & associated temperature, wet weight, %/day
temperature predicted
11.96 0.79

Combined growth rate plot

phenotypes %>% 
  group_by(temperature) %>% summarise(
    SGR.sl.mean=mean(SGR.sl.trt), SGR.sl.sd=sd(SGR.sl.trt),
    SGR.ww.mean=mean(SGR.ww.trt), SGR.ww.sd=sd(SGR.ww.trt)) %>%  
  ggplot() + 
  geom_point(aes(x=as.numeric(as.character(temperature)), y=SGR.sl.mean), size=3) + 
  geom_point(aes(x=as.numeric(as.character(temperature)), y=SGR.ww.mean), size=3) + 
  geom_errorbar(aes(x = as.numeric(as.character(temperature)), ymin = SGR.sl.mean - SGR.sl.sd, ymax = SGR.sl.mean + SGR.sl.sd), width = 0.2) +
  geom_errorbar(aes(x = as.numeric(as.character(temperature)), ymin = SGR.ww.mean - SGR.ww.sd, ymax = SGR.ww.mean + SGR.ww.sd), width = 0.2) +
  theme_minimal() +
  ggtitle("Specific growth rate in treatments, ~% per day") +
  geom_line(data=SGR.sl.predict, aes(x=temperature, y=predicted), color="gray25", size=1) +
  geom_line(data=SGR.ww.predict, aes(x=temperature, y=predicted), color="gray25", size=1) +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_x_continuous(breaks = c(seq(0,16,2))) +
  xlab("Temperature") +   ylab("% / day") + 
  geom_text(aes(x = 13.6, y = 0.83, label="Standard Length"), size=4, color="gray25") +
  geom_text(aes(x = 13.2, y = 0.39, label="Wet Weight"), size=4, color="gray25") +
  geom_vline(xintercept = SGR.sl.predict[which.max(SGR.sl.predict$predicted),]$temperature, 
             linetype = "dashed", color = "gray25", size = 1) +
  geom_vline(xintercept = SGR.ww.predict[which.max(SGR.ww.predict$predicted),]$temperature, 
             linetype = "dashed", color = "gray25", size = 1)

Body condition index at termination (Kwet)

Not significantly affected by temperature (but marginal, p=0.07)

summary(aov(Kwet ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq  Mean Sq F value Pr(>F)  
## temperature   3 0.0199 0.006631   2.384 0.0715 .
## Residuals   152 0.4227 0.002781                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(Kwet ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Kwet ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##               diff         lwr         upr     p adj
## 5-0   0.0006042697 -0.03002831 0.031236850 0.9999515
## 9-0  -0.0023184309 -0.03295101 0.028314149 0.9972972
## 16-0 -0.0272442230 -0.05871621 0.004227762 0.1149886
## 9-5  -0.0029227006 -0.03355528 0.027709880 0.9946319
## 16-5 -0.0278484927 -0.05932048 0.003623492 0.1028730
## 16-9 -0.0249257921 -0.05639778 0.006546192 0.1719787
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(Kwet), sd=sd(Kwet))
## # A tibble: 4 × 3
##   temperature  mean     sd
##   <fct>       <dbl>  <dbl>
## 1 0           0.971 0.0523
## 2 5           0.972 0.0492
## 3 9           0.969 0.0629
## 4 16          0.944 0.0438
phenotypes %>% 
    ggplot(aes(y=Kwet, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Body condition index\n(wet weight / standard length") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Wet weight / standard length")

Hepato-somatic index

Very different among temperature treatments except for 0 vs. 5C.

summary(aov(hsi ~ temperature, phenotypes %>% filter(mort=="No")))
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## temperature   3  49.45  16.482   31.45 7.14e-16 ***
## Residuals   152  79.66   0.524                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(hsi ~ temperature, phenotypes %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hsi ~ temperature, data = phenotypes %>% filter(mort == "No"))
## 
## $temperature
##            diff        lwr        upr     p adj
## 5-0  -0.2335143 -0.6540247  0.1869961 0.4749567
## 9-0  -0.8358065 -1.2563169 -0.4152961 0.0000045
## 16-0 -1.4737284 -1.9057617 -1.0416950 0.0000000
## 9-5  -0.6022922 -1.0228026 -0.1817818 0.0015764
## 16-5 -1.2402141 -1.6722475 -0.8081807 0.0000000
## 16-9 -0.6379219 -1.0699553 -0.2058885 0.0010420
phenotypes %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(hsi), sd=sd(hsi))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            3.61 0.674
## 2 5            3.37 0.877
## 3 9            2.77 0.704
## 4 16           2.13 0.600
phenotypes %>%
    ggplot(aes(y=hsi, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Hepato-somatic index\n(100 * liver wet weight / total wet weight)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
       values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("100 * Liver wet weight / Total wet weight")

# phenotypes %>% filter(mort=="No") %>%
#     ggplot(aes(y=hsi, x=temperature, fill=temperature)) + 
#   geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
#     ggtitle("Appr. hepato-somatic index\n(liver wet weight / total wet weight)") +
#   scale_fill_manual(name="Temperature", 
#                        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c"))

Lipid data!

lipids <- read_excel("../data/Lipid class liver data_091324.xlsx") %>% clean_names() %>%
  mutate(sample_id=as.factor(sample_id))
## New names:
## • `TAG` -> `TAG...20`
## • `FFA` -> `FFA...21`
## • `ST` -> `ST...22`
## • `Polar Lipids` -> `Polar Lipids...23`
## • `Total` -> `Total...24`
## • `TAG` -> `TAG...28`
## • `FFA` -> `FFA...29`
## • `ST` -> `ST...30`
## • `Polar Lipids` -> `Polar Lipids...31`
## • `Total` -> `Total...32`

Add lipid data to phenotype data, and calculate composite performance index

Composite Performance Index calculated as SGR(length) * Kwet * HSI *

phenotypes2 <- phenotypes %>% 
  left_join(lipids %>% dplyr::select(sample_id, genetics_sample_number, tag_20:total_32), 
            by=c("microchip_id"="sample_id")) %>%
    mutate(G.sl.scaled=rescale(SGR.sl.trt, to=c(0,1)),
           G.ww.scaled=rescale(SGR.ww.trt, to=c(0,1)),
           Kwet.scaled=rescale(Kwet, to=c(0,1)),
           hsi.scaled=rescale(hsi, to=c(0,1)),
           total.lipid.scaled=rescale(total_32, to=c(0,1))) %>% 
    mutate(pi=100*G.sl.scaled*G.ww.scaled*Kwet.scaled*hsi.scaled*total.lipid.scaled)

head(phenotypes2)
## # A tibble: 6 × 54
##   microchip_id sl_11212022 wwt_11212022 tank  temperature sl_12272022
##   <fct>              <dbl>        <dbl> <fct> <fct>             <dbl>
## 1 9397                  86         7.38 12    0                    93
## 2 9467                  92         7.98 12    0                    97
## 3 5059                 102        11.6  12    0                   113
## 4 9461                  82         6.15 12    0                    88
## 5 9560                  95         9.59 12    0                   101
## 6 7386                  95         9.23 12    0                   108
## # ℹ 48 more variables: wwt_12272022 <dbl>, mort_date <dttm>,
## #   dissection_date <date>, sl_final <dbl>, wwt_final <dbl>,
## #   total_liver_ww_mg <dbl>, liverfor_lipids_ww_mg <dbl>,
## #   muscle_w_wfor_lipids_mg <dbl>, genetic_sampling_count <fct>,
## #   dissection_comments <chr>, molecular_notes <chr>, fin <chr>, liver <chr>,
## #   blood <chr>, spleen <chr>, gill <chr>, acc_t1 <date>, exp_t1 <date>,
## #   SGR.sl.accl <dbl>, SGR.sl.trt <dbl>, SGR.ww.accl <dbl>, SGR.ww.trt <dbl>, …

Liver lipids consitutents

Total Lipids

summary(aov(total_32 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## temperature  3 237633   79211   7.526 0.000144 ***
## Residuals   94 989295   10524                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 58 observations deleted due to missingness
TukeyHSD(aov(total_32 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = total_32 ~ temperature, data = phenotypes2 %>% filter(mort == "No"))
## 
## $temperature
##            diff         lwr       upr     p adj
## 5-0    43.23684  -33.444705 119.91839 0.4567136
## 9-0    67.56489   -9.116656 144.24644 0.1042676
## 16-0  -62.35665 -139.816708  15.10342 0.1588418
## 9-5    24.32805  -51.567003 100.22310 0.8360115
## 16-5 -105.59349 -182.275040 -28.91194 0.0028107
## 16-9 -129.92154 -206.603089 -53.23999 0.0001469
phenotypes2 %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(total_32, na.rm=T), sd=sd(total_32, na.rm=T))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            318.  87.8
## 2 5            361. 104. 
## 3 9            385. 105. 
## 4 16           255. 112.
phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=total_32, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Total lipid content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Total lipid Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

total.model <- lm(total_32 ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes2)
summary(total.model)
## 
## Call:
## lm(formula = total_32 ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -276.974  -50.619    3.187   68.455  255.134 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                       327.39      10.52  31.107
## poly(as.numeric(as.character(temperature)), 3)1  -332.86     133.57  -2.492
## poly(as.numeric(as.character(temperature)), 3)2  -565.26     133.24  -4.243
## poly(as.numeric(as.character(temperature)), 3)3  -131.43     132.58  -0.991
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1   0.0144 *  
## poly(as.numeric(as.character(temperature)), 3)2 5.14e-05 ***
## poly(as.numeric(as.character(temperature)), 3)3   0.3240    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 104.7 on 95 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.2117, Adjusted R-squared:  0.1868 
## F-statistic: 8.503 on 3 and 95 DF,  p-value: 4.637e-05
total.predict <- data.frame(temperature=seq(0,16,length.out=100))
total.predict$predicted <- predict(total.model, total.predict)

# What is the maximum total lipid and associated predicted temperature as per this model?  0.35%/day at 12C.
total.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted total lipid & associated temperature")
Maximum predicted total lipid & associated temperature
temperature predicted
9.05 385.19 
phenotypes2 %>% 
    ggplot(aes(y=total_32, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Total lipid content") +
  geom_line(data=total.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("% / day")
## Warning: Removed 61 rows containing missing values (`geom_point()`).

Triglycerides

summary(aov(tag_28 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## temperature  3 224741   74914   7.062 0.000248 ***
## Residuals   94 997145   10608                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 58 observations deleted due to missingness
TukeyHSD(aov(tag_28 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = tag_28 ~ temperature, data = phenotypes2 %>% filter(mort == "No"))
## 
## $temperature
##            diff        lwr       upr     p adj
## 5-0    37.71970  -39.26549 114.70490 0.5769094
## 9-0    63.94261  -13.04258 140.92780 0.1385091
## 16-0  -63.71042 -141.47721  14.05636 0.1472629
## 9-5    26.22291  -49.97267 102.41849 0.8047497
## 16-5 -101.43013 -178.41532 -24.44494 0.0046533
## 16-9 -127.65304 -204.63823 -50.66784 0.0002099
phenotypes2 %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(tag_28, na.rm=T), sd=sd(tag_28, na.rm=T))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            293.  87.1
## 2 5            331. 104. 
## 3 9            357. 105. 
## 4 16           230. 113.
phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=tag_28, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Triglyceride content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("TAG Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

tag.model <- lm(tag_28 ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes2)
summary(tag.model)
## 
## Call:
## lm(formula = tag_28 ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -278.176  -52.689    4.626   70.764  256.929 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                       300.64      10.55  28.492
## poly(as.numeric(as.character(temperature)), 3)1  -333.53     133.91  -2.491
## poly(as.numeric(as.character(temperature)), 3)2  -540.77     133.58  -4.048
## poly(as.numeric(as.character(temperature)), 3)3  -143.19     132.92  -1.077
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1 0.014485 *  
## poly(as.numeric(as.character(temperature)), 3)2 0.000105 ***
## poly(as.numeric(as.character(temperature)), 3)3 0.284066    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 105 on 95 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.202,  Adjusted R-squared:  0.1768 
## F-statistic: 8.018 on 3 and 95 DF,  p-value: 8.084e-05
tag.predict <- data.frame(temperature=seq(0,16,length.out=100))
tag.predict$predicted <- predict(tag.model, tag.predict)

# What is the maximum triglyceride and associated predicted temperature as per this model?  0.35%/day at 12C.
tag.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted triglyceride & associated temperature")
Maximum predicted triglyceride & associated temperature
temperature predicted
9.21 357.49 
phenotypes2 %>% 
    ggplot(aes(y=tag_28, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Triglyceride content") +
  geom_line(data=tag.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("TAG Content")
## Warning: Removed 61 rows containing missing values (`geom_point()`).

Free Fatty Acids

summary(aov(ffa_29 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## temperature  3  109.5   36.49   7.448 0.000158 ***
## Residuals   94  460.5    4.90                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 58 observations deleted due to missingness
TukeyHSD(aov(Kwet ~ temperature, phenotypes2 %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Kwet ~ temperature, data = phenotypes2 %>% filter(mort == "No"))
## 
## $temperature
##               diff         lwr         upr     p adj
## 5-0   0.0006042697 -0.03002831 0.031236850 0.9999515
## 9-0  -0.0023184309 -0.03295101 0.028314149 0.9972972
## 16-0 -0.0272442230 -0.05871621 0.004227762 0.1149886
## 9-5  -0.0029227006 -0.03355528 0.027709880 0.9946319
## 16-5 -0.0278484927 -0.05932048 0.003623492 0.1028730
## 16-9 -0.0249257921 -0.05639778 0.006546192 0.1719787
phenotypes2 %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(total_32, na.rm=T), sd=sd(total_32, na.rm=T))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            318.  87.8
## 2 5            361. 104. 
## 3 9            385. 105. 
## 4 16           255. 112.
phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=ffa_29, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Free fatty acid content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("FFA Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

ffa.model <- lm(ffa_29 ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes2)
summary(ffa.model)
## 
## Call:
## lm(formula = ffa_29 ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8355 -1.4248 -0.3153  1.4970  5.7127 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                       4.4480     0.2219  20.044
## poly(as.numeric(as.character(temperature)), 3)1  -7.7234     2.8162  -2.742
## poly(as.numeric(as.character(temperature)), 3)2 -11.0858     2.8093  -3.946
## poly(as.numeric(as.character(temperature)), 3)3  -2.4583     2.7953  -0.879
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1 0.007289 ** 
## poly(as.numeric(as.character(temperature)), 3)2 0.000152 ***
## poly(as.numeric(as.character(temperature)), 3)3 0.381380    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.208 on 95 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.2029, Adjusted R-squared:  0.1777 
## F-statistic: 8.061 on 3 and 95 DF,  p-value: 7.692e-05
ffa.predict <- data.frame(temperature=seq(0,16,length.out=100))
ffa.predict$predicted <- predict(ffa.model, ffa.predict)

# What is the maximum FFA and associated predicted temperature as per this model?  0.35%/day at 12C.
ffa.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted free fatty acid & associated temperature")
Maximum predicted free fatty acid & associated temperature
temperature predicted
8.73 5.55 
phenotypes2 %>% 
    ggplot(aes(y=ffa_29, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Free fatty acid content") +
  geom_line(data=ffa.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("FFA Content")
## Warning: Removed 61 rows containing missing values (`geom_point()`).

Sterols

summary(aov(st_30 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## temperature  3  14.78   4.928   10.58 4.67e-06 ***
## Residuals   94  43.80   0.466                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 58 observations deleted due to missingness
TukeyHSD(aov(st_30 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = st_30 ~ temperature, data = phenotypes2 %>% filter(mort == "No"))
## 
## $temperature
##             diff         lwr         upr     p adj
## 5-0   1.09446739  0.58425391  1.60468086 0.0000012
## 9-0   0.55701862  0.04680514  1.06723209 0.0266961
## 16-0  0.63830828  0.12291484  1.15370172 0.0088463
## 9-5  -0.53744877 -1.04242915 -0.03246839 0.0323788
## 16-5 -0.45615911 -0.96637258  0.05405437 0.0966813
## 16-9  0.08128966 -0.42892381  0.59150314 0.9754894
phenotypes2 %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(st_30, na.rm=T), sd=sd(st_30, na.rm=T))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            2.23 0.471
## 2 5            3.32 0.923
## 3 9            2.78 0.695
## 4 16           2.87 0.536
phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=st_30, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Sterol content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("sterol Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

st.model <- lm(st_30 ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes2)
summary(st.model)
## 
## Call:
## lm(formula = st_30 ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.26046 -0.32453  0.07564  0.35719  2.14094 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                       2.7946     0.0685  40.798
## poly(as.numeric(as.character(temperature)), 3)1   1.7988     0.8693   2.069
## poly(as.numeric(as.character(temperature)), 3)2  -3.1500     0.8672  -3.633
## poly(as.numeric(as.character(temperature)), 3)3   3.3056     0.8629   3.831
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1 0.041240 *  
## poly(as.numeric(as.character(temperature)), 3)2 0.000455 ***
## poly(as.numeric(as.character(temperature)), 3)3 0.000229 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6815 on 95 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.2503, Adjusted R-squared:  0.2266 
## F-statistic: 10.57 on 3 and 95 DF,  p-value: 4.607e-06
st.predict <- data.frame(temperature=seq(0,16,length.out=100))
st.predict$predicted <- predict(st.model, st.predict)

# What is the maximum sterol and associated predicted temperature as per this model?  0.35%/day at 12C.
st.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted sterol & associated temperature")
Maximum predicted sterol & associated temperature
temperature predicted
4.53 3.33 
phenotypes2 %>% 
    ggplot(aes(y=st_30, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Sterol content") +
  geom_line(data=st.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("Sterol Content")
## Warning: Removed 61 rows containing missing values (`geom_point()`).

Polar Lipids

summary(aov(polar_lipids_31 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##             Df Sum Sq Mean Sq F value Pr(>F)  
## temperature  3  165.5   55.18   2.319 0.0804 .
## Residuals   94 2236.8   23.80                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 58 observations deleted due to missingness
TukeyHSD(aov(polar_lipids_31 ~ temperature, phenotypes2 %>% filter(mort=="No")))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = polar_lipids_31 ~ temperature, data = phenotypes2 %>% filter(mort == "No"))
## 
## $temperature
##           diff          lwr      upr     p adj
## 5-0   3.637550 -0.008653012 7.283752 0.0507837
## 9-0   1.893293 -1.752909255 5.539496 0.5286287
## 16-0  2.308196 -1.375024451 5.991417 0.3619231
## 9-5  -1.744256 -5.353060907 1.864548 0.5878696
## 16-5 -1.329353 -4.975555807 2.316849 0.7759266
## 16-9  0.414903 -3.231299565 4.061106 0.9907754
phenotypes2 %>% filter(mort=="No") %>% group_by(temperature) %>% dplyr::summarise(mean=mean(total_32, na.rm=T), sd=sd(total_32, na.rm=T))
## # A tibble: 4 × 3
##   temperature  mean    sd
##   <fct>       <dbl> <dbl>
## 1 0            318.  87.8
## 2 5            361. 104. 
## 3 9            385. 105. 
## 4 16           255. 112.
phenotypes2 %>% 
  filter(!is.na(total_32)) %>%
    ggplot(aes(y=polar_lipids_31, x=interaction(mort, temperature), fill=temperature, shape=mort, alpha=mort, label=genetic_sampling_count)) + 
  geom_violin() + geom_jitter(width=0.25) + theme_minimal() +
    ggtitle("Polar lipid content (ug/mg)") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  xlab(NULL) + ylab("Polar lipid Content")
## Warning: Groups with fewer than two data points have been dropped.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

Composite Performance Index

phenotypes2 %>% ggplot(aes(x=temperature, y=pi, fill=temperature)) + 
  geom_violin(alpha=0.6) + geom_jitter(width = 0.2) + 
    scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme_minimal() + ylab("Composite Performance Index\nGSR(sl) * GSR(ww) * Kwet * HSI * Total Lipid")

Model composite performance index by temperature

pi.model <- lm(pi ~ poly(as.numeric(as.character(temperature)), 3), data = phenotypes2)
summary(pi.model)
## 
## Call:
## lm(formula = pi ~ poly(as.numeric(as.character(temperature)), 
##     3), data = phenotypes2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1959  -2.7916  -0.3744   1.6485  21.6455 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                       4.8652     0.4882   9.965
## poly(as.numeric(as.character(temperature)), 3)1  13.5395     6.1960   2.185
## poly(as.numeric(as.character(temperature)), 3)2 -35.0231     6.1806  -5.667
## poly(as.numeric(as.character(temperature)), 3)3 -18.6743     6.1500  -3.036
##                                                 Pr(>|t|)    
## (Intercept)                                      < 2e-16 ***
## poly(as.numeric(as.character(temperature)), 3)1  0.03133 *  
## poly(as.numeric(as.character(temperature)), 3)2 1.55e-07 ***
## poly(as.numeric(as.character(temperature)), 3)3  0.00309 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.857 on 95 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.3263, Adjusted R-squared:  0.305 
## F-statistic: 15.34 on 3 and 95 DF,  p-value: 3.256e-08
pi.predict <- data.frame(temperature=seq(0,16,length.out=100))
pi.predict$predicted <- predict(pi.model, pi.predict)

# What is the maximum standard length growth rate and associated predicted temperature as per this model?  0.35%/day at 12C.
pi.predict %>% arrange(desc(predicted)) %>% head(n=1) %>% round(digits = 2) %>% kable(caption = "Maximum predicted composite performance index & associated temperature, wet weight, %/day")
Maximum predicted composite performance index & associated temperature, wet weight, %/day
temperature predicted
11.31 11.58 
phenotypes2 %>% 
    ggplot(aes(y=pi, x=as.numeric(as.character(temperature)))) + 
  geom_point(aes(color=temperature)) + theme_minimal() +
  ggtitle("Composite Performance Index") +
  geom_line(data=pi.predict, aes(x=temperature, y=predicted), color="gray45", size=1) +
  #  stat_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 3),colour="gray50") +
  scale_alpha_manual(values=c(0.75,0.5)) +
  scale_color_manual(name="Temperature", values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  xlab(NULL) +   ylab("% / day")

phenotypes2 %>% dplyr::select(microchip_id, genetic_sampling_count, temperature, G.sl.scaled:pi) %>% 
  pivot_longer(cols = c(G.sl.scaled:total.lipid.scaled), names_to = "biometric", values_to = "scaled_value") %>% 
  ggplot(aes(y=pi, x=scaled_value, color=temperature)) + geom_point() + 
  facet_wrap(~biometric) + 
  theme_minimal_grid() + xlab("Scaled value") + ylab("Composite Performancde Index") + 
    scale_color_manual(name="Temperature", 
        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  geom_smooth(method = "glm", se = F) + theme(legend.position = "bottomright")

phenotypes3 <- phenotypes2 %>% select(microchip_id, genetic_sampling_count, temperature, SGR.sl.trt, SGR.ww.trt, Kwet, hsi, mort, tag_28, ffa_29, st_30, polar_lipids_31, total_32, pi) %>%
  filter(mort=="No")
save(phenotypes3, file="../data/phenotypes3")

All biometrics in one plot

phenotypes3 %>% filter(mort=="No") %>%  dplyr::select(-mort) %>% 
  dplyr::rename("K(wet)"="Kwet", "HSI"="hsi", "Triglycerides"="tag_28", "Free Fatty Acids"="ffa_29", 
         "Sterols"="st_30", "Polar Lipids"="polar_lipids_31", "Total Lipids"="total_32",
         "Performance Index"="pi", "Growth Rate, Length"="SGR.sl.trt",
         "Growth Rate, Weight"="SGR.ww.trt") %>%  
  select(-`Performance Index`) %>% 
  pivot_longer(`Growth Rate, Length`:HSI, names_to = "metric", values_to = "value") %>% 
  mutate(metric=as.factor(metric)) %>%
 
  ggplot(aes(x=temperature, y=value, fill=temperature)) + 
  #geom_jitter(width=0.25, size=1.5) + 
  geom_boxplot(alpha=0.5) + 
  theme_minimal() +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~metric, scales = "free", nrow=3) + xlab(NULL) + ylab(NULL)

### All lipid metrics in one plot

phenotypes3 %>% filter(mort=="No") %>%  dplyr::select(-mort) %>% 
  dplyr::rename("K(wet)"="Kwet", "HSI"="hsi", "Triglycerides"="tag_28", "Free Fatty Acids"="ffa_29", 
         "Sterols"="st_30", "Polar Lipids"="polar_lipids_31", "Total Lipids"="total_32",
         "Performance Index"="pi", "Growth Rate, Length"="SGR.sl.trt",
         "Growth Rate, Weight"="SGR.ww.trt") %>%  
  select(-`Performance Index`) %>% 
  pivot_longer(`Triglycerides`:`Total Lipids`, names_to = "metric", values_to = "value") %>% 
  mutate(metric=as.factor(metric)) %>%
 
  ggplot(aes(x=temperature, y=value, fill=temperature)) + 
  #geom_jitter(width=0.25, size=1.5) + 
  geom_boxplot(alpha=0.5) + 
  theme_minimal() +
  scale_fill_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~metric, scales = "free", nrow=3) + xlab(NULL) + ylab(NULL)

### Correlations among metrics

require(corrplot)
corr.metrics <- cor(phenotypes2 %>% select(SGR.sl.trt, SGR.ww.trt, Kwet, hsi, mort, days_growth, tag_28, ffa_29, st_30, polar_lipids_31, total_32, pi) %>%
  filter(mort=="No") %>% mutate(growthrate.sl=SGR.sl.trt/days_growth, growthrate.wwt=SGR.ww.trt/days_growth) %>% select(-SGR.sl.trt, -SGR.ww.trt, -mort, -days_growth), use="complete.obs")

corrplot::corrplot(corr.metrics, method = "number", type="lower", diag=T, tl.cex=.8, number.cex = 0.8, sig.level = 0.05, tl.pos="tl")
corrplot::corrplot(corr.metrics, method = "color", type="upper", diag=F, sig.level = 0.05, add=T, tl.pos="n", )

Residuals, all metrics ~ temperature

a <- lm(total_32~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
b <- lm(polar_lipids_31~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
c <- lm(st_30~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
d <- lm(ffa_29~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
e <- lm(tag_28~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
f <- lm(Kwet~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
g <- lm(hsi~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
h <- lm(`SGR.sl.trt`~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))
i <- lm(`SGR.ww.trt`~as.numeric(as.character(temperature)), data=phenotypes3 %>% filter(!is.na(total_32)))

phenotypes3.resids <- cbind(
  sample=as.character((phenotypes3 %>% filter(!is.na(total_32)))$microchip_id),
  temperature=as.character((phenotypes3 %>% filter(!is.na(total_32)))$temperature),        
        resid_total.lipids=resid(a),
         resid_polar=resid(b),
         resid_sterol=resid(c),
         resid_ffa=resid(d),
         resid_trigl=resid(e),
         resid_Kwet=resid(f),
         resid_HSI=resid(g),
         resid_sl.growthrate=resid(h),
         resid_wwt.growthrate=resid(i)) %>% as.data.frame()
phenotypes3.resids %>% 
  mutate(temperature=factor(temperature, ordered=T, levels=c(0,5,9,16))) %>%
  pivot_longer(starts_with("resid"), names_to = "metric", values_to = "resid") %>% 
  mutate(temperature=as.factor(temperature), metric=as.factor(metric), resid=as.numeric(resid)) %>% 
 
  ggplot(aes(x=temperature, y=resid, color=temperature, text=sample)) + 
  geom_jitter(width=0.25, size=1.5) + theme_minimal() +
  scale_color_manual(name="Temperature", 
     values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  facet_wrap(~metric, scales = "free")

Correlations, residuals ~ temperature

corr.metrics.resids <- cor(phenotypes3.resids %>% select(-sample, -temperature) %>% mutate_all(as.numeric), use="complete.obs")

corrplot::corrplot(corr.metrics.resids, method = "number", type="lower", diag=T, tl.cex=.8, number.cex = 0.8, sig.level = 0.05, tl.pos="tl")
corrplot::corrplot(corr.metrics.resids, method = "color", type="upper", diag=F, sig.level = 0.05, add = TRUE, tl.pos = "n")

# anova(a) #total lipid
# anova(b) #polar lipids <-- only metric not affected by temperature 
# anova(c) #sterols 
# anova(d) #free fatty acids 
# anova(e) #triglycerols 
# anova(f) #condition index
# anova(g) #HSI
# anova(h) #growth rate - length
# anova(i) #growth rate - weight

Look at total lipids ~ condition index

#ggplotly(
  phenotypes3 %>% 
    ggplot(aes(y=total_32, x=Kwet, color=temperature)) + 
  geom_point() + theme_minimal() +
    ggtitle("Total Lipids / Condition Index") +
  scale_color_manual(name="Temperature", 
        values=c("0"="#2c7bb6", "5"="#abd9e9", "9"="#fdae61", "16"="#d7191c")) +
  geom_smooth(method = "lm") +   facet_wrap(~temperature) +
  xlab("Condition Index") + ylab("Total Lipids")#)