---
title: "Social familiarity improves fast-start escape performance in schooling fish"
author: "Lauren E Nadler"
date: "June 10, 2021"
output: html_document
---
Using social groups (i.e. schools) of the tropical damselfish Chromis viridis, we tested how familiarity through repeated social interactions influences fast-start responses, the primary defensive behaviour in a range of taxa, including fish, sharks, and larval amphibians. We focused on reactivity through response latency and kinematic performance (i.e. agility and propulsion) following a simulated predator attack, while distinguishing between first and subsequent responders (direct response to stimulation versus response triggered by integrated direct and social stimulation, respectively). Schools composed of eight familiar or unfamiliar C. viridis individuals (n = 12 schools per treatment, referred to as familiar and unfamiliar schools, respectively, hereafter) were tested in a laminar flow tank, using a current speed (3.2 cm/s) that mimicked their natural habitat conditions on a calm weather day.
```{r, message = F}
#LIBRARIES
library(Matrix)
library(lme4)
library(carData)
library(car)
library(MuMIn)
library(emmeans)
library(mvtnorm)
library(survival)
library(TH.data)
library(MASS)
library(multcomp)
library(partR2)
library(patchwork)
library(msm)
library(polycor)
library(ltm)
library(ggplot2)
library(zoo)
library(lmtest)
library(glmmTMB)
library(ARTool)
library(optimx)
```
# LATENCY
Individuals were assigned a responder number based on the sequential order of latencies in their respective group. Responder number was designated using standard competition ranking (e.g. if two individuals tied for first place, the responder numbers for that group would be 1, 1, 3, 4, 5, 6, 7, 8; sample sizes for each responder number by treatment are detailed in Table S1; ties occurred in 13 out of 24 schools tested; number of ties: 2-4 individuals per school; number of schools exhibiting ties by treatment: n=9 familiar schools, n=4 unfamiliar schools). While the responses for non-first responders (i.e. responder numbers 2 through 8) could have been initiated either directly by the stimulus or indirectly by their group-mates´ responses, first responders´ responses could only have been initiated directly by the stimulus.
```{r, message = F, warning=FALSE}
#dotplot of latency versus responder number
escape <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape.csv", header=T)
escape$Responder.number = as.factor(escape$Responder.number)
ggplot(escape, aes(x=Responder.number, y=Latency)) + geom_dotplot(aes(color=School.Type, fill=School.Type), binaxis = "y", stackdir = "center", dotsize=0.35) + xlab(bquote('Responder number')) + scale_y_log10(name ="Latency (ms)") + theme_bw() + scale_colour_manual(values = c("#386cb0","#fdb462"), name ="School.Type") + scale_fill_manual(values = c("#386cb0","#fdb462"), name ="School.Type") + theme(legend.title = element_blank()) + theme(legend.position = c(0.2, 0.9)) + theme(legend.box = "horizontal") + theme(legend.text=element_text(size=25)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 25))
#latency - first responders model
ties <- read.csv("/Users/laurennadler/Desktop/MyRData/Familiarity_First responder ties-3.csv", header=T)
summary(ties)
first.glm = glm(Latency~School.Type*StimDist, data=ties, weights=Frequency)
summary(first.glm)
anova(first.glm, test = "F")
plot(first.glm)
shapiro.test(resid(first.glm))
bartlett.test(resid(first.glm), ties$School.Type)
r.squaredGLMM(first.glm)
biserial.cor(ties$Latency, ties$School.Type)
#latency - subsequent responders model
subsequent <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape-3.csv", header=T)
summary(subsequent)
latsub.lmer = lmer(Latency ~ School.Type*Responder.number*StimDist + (1|School), data=subsequent, na.action = "na.omit")
plot(latsub.lmer)
qqnorm(resid(latsub.lmer))
qqline(resid(latsub.lmer))
shapiro.test(resid(latsub.lmer))
bartlett.test(resid(latsub.lmer), subsequent$School.Type) #doesn't meet assumptions so using boxcox transformation
latsub.lmer.bc <- powerTransform(latsub.lmer, family="bcPower")
summary(latsub.lmer.bc)
latsub.bc.lmer = lmer(bcPower(Latency, latsub.lmer.bc$roundlam) ~ School.Type*Responder.number*StimDist + (1|School), data=subsequent, na.action = "na.omit")
plot(latsub.bc.lmer)
qqnorm(resid(latsub.bc.lmer))
qqline(resid(latsub.bc.lmer))
shapiro.test(resid(latsub.bc.lmer))
bartlett.test(resid(latsub.bc.lmer), subsequent$School.Type)
summary(latsub.bc.lmer)
Anova(latsub.bc.lmer, test="F")
r.squaredGLMM(latsub.bc.lmer)
sub <- partR2(latsub.bc.lmer, partvars = c("School.Type", "Responder.number", "School.Type:Responder.number"), nboot=1000)
print(sub)
summary(sub, round_to = 2)
forestplot(sub, type = "SC")
```
The fast-start response is commonly initiated by a pair of higher order command neurons called Mauthner cells (M-cells). If M-cell functionality is lost through experimental ablation, rapid fast-start responses (e.g. latency <16ms from the stimulus onset, as measured in Eaton et al. 1981, J Comp Physiol A) are absent, though longer latency, slower responses can be initiated by other reticulospinal neurons in the brainstem escape network. Thus, we also examined if there was a significant difference in the proportion of fast versus slow responses with familiarity.
```{r, message = F, warning=FALSE}
#fast-slow model
escape <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape.csv", header=T)
summary(escape)
latency.lmer = glmer(Fast.Slow.numeric ~ School.Type + (1|School), data=escape, family = "binomial")
summary(latency.lmer)
Anova(latency.lmer)
r.squaredGLMM(latency.lmer)
fastslow <- partR2(latency.lmer, max_level = 1, nboot=1000)
print(fastslow)
summary(fastslow, round_to = 2)
forestplot(fastslow, type = "SC")
```
# AVERAGE TURNING RATE & DISTANCE COVERED
We quantified fast-start characteristics associated with the initial unilateral body bend following stimulation (i.e., stage 1 of the fast-start) and the subsequent contralateral body bend (i.e., stage 2 of the fast-start). Individual fast-start kinematic performance was characterized through average turning rate (i.e. the ratio of the turning angle achieved during stage 1 to the stage 1 duration) and distance covered (i.e. distance moved by the centre of mass during the first 42 ms of the reaction, the mean time to achieve stages 1 and 2 in this species according to published data in Nadler et al. 2018, Biol Open).
```{r, message = F, warning=FALSE}
#dotplot of average turning rate versus responder number
escape <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape.csv", header=T)
escape$Responder.number = as.factor(escape$Responder.number)
turning <- ggplot(escape, aes(x=Responder.number, y=Turning)) + geom_dotplot(aes(color=School.Type, fill=School.Type), binaxis = "y", stackdir = "center", dotsize=0.35) + xlab(bquote('Responder number')) + ylab(bquote("Average turning rate (°/s)")) + theme_bw()
turning + scale_colour_manual(values = c("#386cb0","#fdb462"), name ="School.Type") + scale_fill_manual(values = c("#386cb0","#fdb462"), name ="School.Type") + theme(legend.title = element_blank()) + theme(legend.position = "none") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 25))
#turning - first responders model
first <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape-4.csv", header=T)
summary(first)
turn1st.lmer = lmer(Turning ~ Treatment*Stim.Dist + (1|School), data=first, na.action = "na.omit")
plot(turn1st.lmer)
qqnorm(resid(turn1st.lmer))
qqline(resid(turn1st.lmer))
shapiro.test(resid(turn1st.lmer))
bartlett.test(resid(turn1st.lmer), first$Treatment)
summary(turn1st.lmer)
Anova(turn1st.lmer, test = "F")
r.squaredGLMM(turn1st.lmer)
turn1 <- partR2(turn1st.lmer, nboot=1000)
print(turn1)
summary(turn1, round_to = 2)
forestplot(turn1, type = "SC")
#turning - subsequent responders models
subsequent <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape-3.csv", header=T)
summary(subsequent)
turnsub.lmer = lmer(Turning ~ School.Type*Responder.number*StimDist + (1|School), data=subsequent, na.action = "na.omit")
plot(turnsub.lmer)
qqnorm(resid(turnsub.lmer))
qqline(resid(turnsub.lmer))
shapiro.test(resid(turnsub.lmer))
bartlett.test(resid(turnsub.lmer), subsequent$School.Type)
summary(turnsub.lmer)
Anova(turnsub.lmer, test = "F")
r.squaredGLMM(turnsub.lmer)
turnsub <- partR2(turnsub.lmer, partvars = c("School.Type", "Responder.number", "StimDist"), nboot=1000)
print(turnsub)
summary(turnsub, round_to = 2)
forestplot(turnsub, type = "SC")
```
Note that for both distance covered datasets I had to remove the missing values to calculate the structure coefficient effect sizes (Stoffel et al. 2021 BioRxiv specifies that partR2 can´t handle missing values). For distance covered, individuals that were less than 1 body length from a wall of the experimental arena were excluded from this analysis due to the constraining effects of the wall on their ability to accelerate with their response.
```{r, message = F, warning=FALSE}
#dotplot of distance covered versus responder number
escape <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape.csv", header=T)
escape$Responder.number = as.factor(escape$Responder.number)
ggplot(escape, aes(x=Responder.number, y=DistCovered)) + geom_dotplot(aes(color=School.Type, fill=School.Type), binaxis = "y", stackdir = "center", dotsize=0.35) + xlab(bquote('Responder number')) + ylab(bquote("Distance covered (mm)")) + theme_bw() + scale_colour_manual(values = c("#386cb0","#fdb462"), name ="School.Type") + scale_fill_manual(values = c("#386cb0","#fdb462"), name ="School.Type") + theme(legend.title = element_blank()) + theme(legend.position = "none") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 25))
#distance - first responders model
first.dist.noNA <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape-6.csv", header=T)
summary(first.dist.noNA)
dist1st.lmer = lmer(DistCovered ~ Treatment*Stim.Dist + (1|School), data=first.dist.noNA, na.action = "na.omit")
summary(dist1st.lmer)
plot(dist1st.lmer)
qqnorm(resid(dist1st.lmer))
qqline(resid(dist1st.lmer))
shapiro.test(resid(dist1st.lmer))
bartlett.test(resid(dist1st.lmer), first.dist.noNA$Treatment)
Anova(dist1st.lmer, test="F")
r.squaredGLMM(dist1st.lmer)
dist1 <- partR2(dist1st.lmer, partvars = c("Treatment", "Stim.Dist", "Treatment:Stim.Dist"), nboot=1000)
print(dist1)
summary(dist1, round_to = 2)
forestplot(dist1, type = "SC")
#distance - subsequent responders model
subs.dist.noNA <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscape-5.csv", header=T)
summary(subs.dist.noNA)
distsub.lmer = lmer(DistCovered ~ School.Type*Responder.number*StimDist + (1|School), data=subs.dist.noNA, na.action = "na.omit")
summary(distsub.lmer)
plot(distsub.lmer)
qqnorm(resid(distsub.lmer))
qqline(resid(distsub.lmer))
shapiro.test(resid(distsub.lmer))
bartlett.test(resid(distsub.lmer), subs.dist.noNA$School.Type)
Anova(distsub.lmer, test="F")
r.squaredGLMM(distsub.lmer)
distsub <- partR2(distsub.lmer, partvars = c("School.Type", "Responder.number", "School.Type:Responder.number"), nboot=1000)
print(distsub)
summary(distsub, round_to = 2)
forestplot(distsub, type = "SC")
```
# SCHOOL BEHAVIOUR AFTER STIMULATION
School cohesion and coordination were measured during the 100 ms time period following the simulated predator attack. This short time frame following a predator attack is crucial for survival. School cohesion was assessed using nearest neighbour distance (NND; distance to the closest neighbour) and school area (horizontal spread of the school), while school alignment (length of the mean circular vector, a measure of the variation in fish orientations) was used as a proxy for spatial coordination of individuals within the school. These school characteristics were measured at set time points throughout the escape response (Time = 0, ~20, and ~100 ms after stimulus onset), representative of before, during, and after the escape response on average in this species.
```{r message=FALSE, warning=FALSE}
#nearest neighbor distance model
NND <- read.csv("/Users/laurennadler/Desktop/MyRData/FamEscapeNND.csv", header=T)
NND$Time = as.factor(NND$Time)
summary(NND)
nnd.lmer = lmer(NND ~ School.Type*Time + (1|School), data=NND)
plot(nnd.lmer)
qqnorm(resid(nnd.lmer))
qqline(resid(nnd.lmer))
shapiro.test(resid(nnd.lmer))
bartlett.test(resid(nnd.lmer), NND$Time) #doesn't meet assumptions so using boxcox transformation
nnd.lmer.bc <- powerTransform(nnd.lmer, family="bcPower")
summary(nnd.lmer.bc)
nnd.bc.lmer = lmer(bcPower(NND, nnd.lmer.bc$roundlam) ~ School.Type*Time + (1|School), data=NND)
plot(nnd.bc.lmer)
qqnorm(resid(nnd.bc.lmer))
qqline(resid(nnd.bc.lmer))
shapiro.test(resid(nnd.bc.lmer))
bartlett.test(resid(nnd.bc.lmer), NND$School.Time)
summary(nnd.bc.lmer)
Anova(nnd.bc.lmer, test = "F")
r.squaredGLMM(nnd.bc.lmer)
nndpartR2 <- partR2(nnd.bc.lmer, nboot=1000)
print(nndpartR2)
summary(nndpartR2, round_to = 2)
forestplot(nndpartR2, type = "SC")
# *NND tukeys test for figure lettering
nnd.bc2.lmer = lmer(bcPower(NND, nnd.lmer.bc$roundlam) ~ School.Time + (1|School), data=NND)
summary(glht(nnd.bc2.lmer, linfct=mcp(School.Time="Tukey")))
#school alignment model
align <- read.csv("/Users/laurennadler/Desktop/MyRData/EscapeFamAlignment.csv", header=T)
align$Time = as.factor(align$Time)
summary(align)
align.lmer = lmer(Alignment ~ Treatment*Time + (1|School), data=align)
plot(align.lmer)
qqnorm(resid(align.lmer))
qqline(resid(align.lmer))
shapiro.test(resid(align.lmer))
bartlett.test(resid(align.lmer), align$Treatment.Time)
summary(align.lmer)
Anova(align.lmer, test="F")
r.squaredGLMM(align.lmer)
alignpartR2 <- partR2(align.lmer, nboot=1000)
print(alignpartR2)
summary(alignpartR2, round_to = 2)
forestplot(alignpartR2, type = "SC")
# *Alignment tukeys test for figure lettering
align2.lmer = lmer(Alignment ~ Treatment.Time + (1|School), data=align)
summary(glht(align2.lmer, linfct=mcp(Treatment.Time="Tukey")))
#school area model
area <- read.csv("/Users/laurennadler/Desktop/MyRData/EscapeFamArea.csv", header=T)
area$Time = as.factor(area$Time)
summary(area)
area.lmer = lmer(Area ~ Treatment*Time + (1|School), data=area)
plot(area.lmer)
qqnorm(resid(area.lmer))
qqline(resid(area.lmer))
shapiro.test(resid(area.lmer))
bartlett.test(resid(area.lmer), area$Treatment.Time) #doesn't meet assumptions so trying box cox transformation
area.lmer.bc <- powerTransform(area.lmer, family="bcPower")
summary(area.lmer.bc) #indicates that log transformation appropriate (LRT test indicates transformation parameter - lambda - = 0)
logarea.lmer = lmer(log(Area) ~ Treatment*Time + (1|School), data=area)
Anova(logarea.lmer, test = "F")
r.squaredGLMM(logarea.lmer)
areapartR2 <- partR2(logarea.lmer, nboot=1000)
print(areapartR2)
summary(areapartR2, round_to = 2)
forestplot(areapartR2, type = "SC")
# *Area tukeys test for figure lettering
area2.lmer = lmer(log(Area) ~ Treatment.Time + (1|School), data=area)
summary(glht(area2.lmer, linfct=mcp(Treatment.Time="Tukey")))
```
# TIME TO FAMILIARITY
Prior to escape response testing, we confirmed that our study species achieves familiarity in a comparable time frame to another tropical fish species, the guppy Poecilia reticulata (Griffiths & Magurran 1997, Anim Behav), using a choice test methodology (supplementary methods and results related to this preliminary study included in the supplementary information of the paper).
```{r message=FALSE, warning=FALSE}
#time to familiarity model
time <- read.csv("/Users/laurennadler/Desktop/MyRData/Familiarity_time to fam.csv", header=T)
summary(time)
time.glm = glm(Pref.fam. ~ Day + School, family = binomial, data=time)
summary(time.glm)
anova(time.glm, test = "Chisq")
```