Physiological state matching in a pair bonded poison frog

Data analyses and R script

General aim and experimental design

Here, we tested for evidence of emotional contagion in a pair bonding amphibian, the mimetic poison frog (Ranitomeya imitator). In this species, males and females form prolonged partnerships characterized by affiliative interactions, the mutual defense of a shared territory, and joint care for offspring. By subjecting male-female partner dyads to an “empathy assay” (see Figure 1) similar to one developed for rodents, we tested the hypothesis that male observers would display emotional contagion and ingroup bias towards female partners (demonstrators) that were subjected to a stressor. Specifically, we predicted that males would state match female partners hormonally and behaviorally despite never experiencing or observing the stressor themselves, and that this response would be biased towards partners relative to familiar non-partner females. Furthermore, we predicted that this response would increase with the longevity and lifetime reproductive output of partnerships. To better interpret the results from the stress experiment and establish the semi-natural (baseline) corticosterone state of pairs, we also examined corticosterone levels and hormone matching of experimentally naïve pairs during cohabitation.

Figure 1: Design for testing partner-selective emotional contagion in pair bonded male R. imitator poison frogs. Pair bonded males were assayed for state matching the behavior and corticosterone level of partner females that underwent a stress treatment that males did not observe or experience themselves. To examine partner-specificity, each male was assayed with their female partner and a familiar female non-partner of similar reproductive salience. To better interpret these experimental results and establish the semi-natural (baseline) corticosterone state of pairs, we also examined corticosterone levels and matching of experimentally naïve pairs within their housing terraria. Picture is of a pair bond within its housing terrarium, taken by Daniel Shaykevich.

Libraries

The following packages were used to perform the analyses (display code):

library(ggplot2)
library(emmeans)
library(lmerTest)
library(lme4)
library(png)
library(grid)
library(lmtest)
library(MuMIn)
library(ggpubr)
library(nlme)
library(GGally)
library(PerformanceAnalytics)
library(psych)
library(knitr)
library(performance)
library(see)
library(viridis)
library(car)
library(rstatix)
library(FSA)
library(repmis)
library(kableExtra)
library(webshot)
library(Matrix)
library(tidyverse)
library(DHARMa)
library(sjPlot)
library(dlookr)
library(coin)

Graphical parameters setting

We set some general graphical parameters that were used for most of the figures:

my_theme <- theme_bw() + theme(panel.grid.major = element_blank(), 
                  panel.grid.minor = element_blank(),
                  axis.text.x = NULL, legend.position = "none") 

my_geom_smooth <- geom_smooth(method="lm", col="black", lwd=0.5)

my_geom_point <- geom_point(shape=21, bg="grey", size=3.5)

Data

We used the following data sets for all the analyses (https://doi.org/10.5061/dryad.05qfttfbn):

1- ‘Empathy’, which has the concentration of corticosterone of every individual and information related to their experimental condition (i.e. partner/non-partner, stress/non-stress, baseline/treatment/housing tank).

2- ‘CORTcorr’, which its a ‘widened’ version of the ‘Empathy’ data set, where corticosterone concentrations are paired based on the experimental dyads. Because some of the hormone samples did not meet the quality controls (i.e. intra-assay CV > 20%) we kept only pairs that had pre and post-experimental measurements in all four experimental conditions (partner + non-stress, partner + stress, non-partner + non-stress, non-partner + stress).

3- ‘bhs’, which has values of 9 scored individual and joint (male-female dyad) behaviors of every individual frog in every experimental condition (Table 1).

cod.beh <- data.frame("Behavior" = c("Gaze - bouts", "Gaze - duration", 
                                   "Approach - bouts", "Approach - duration",
                                   "Freeze - bouts", "Freeze - duration",
                                   "Refuge use - duration", "Activity - bouts",
                                   "Activity - speed", "Activity - distance",
                                   "Space use", "Male-female proximity",
                                   "Male-female joint refuge use", 
                                   "Coordinated freezing"),
                      "Description" = c("Total number of times animal is facing towards other conspecific", "Total duration animal is facing towards other conspecific",
                                      "Total number of times animal moves towards other conspecific",
                                      "Total duration of times animal moves towards other conspecific", 
                                      "Total number of times animal is motionless for at least 1 min.  and completely visible",
                                      "Total duration of time animal is motionless for at least 1 min.  and completely visible",
                                      "Animal is attempting to be concealed within moss or in cup",
                                      "Total number of times animal is in motion, either stationary or through space",
                                      "Average speed animal has moved through space",
                                      "Total distance animal has moved",
                                      "Total area of space that the animal used (cm2)",
                                      "The average distance between animals",
                                      "Both animals are refuging together in moss or cup",
                                      "Both animals are freezing at same time"),
                      "Scoring software" = c(rep("BORIS",8), rep("Annolid",4), rep("BORIS",2)))

cod.beh %>% 
  kable(digits = 3, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, bold = T) %>% 
  pack_rows("Individual",1,11) %>% 
  pack_rows("Male-female dyad",12,14)

Behavior	Description	Scoring.software
Individual
Gaze - bouts	Total number of times animal is facing towards other conspecific	BORIS
Gaze - duration	Total duration animal is facing towards other conspecific	BORIS
Approach - bouts	Total number of times animal moves towards other conspecific	BORIS
Approach - duration	Total duration of times animal moves towards other conspecific	BORIS
Freeze - bouts	Total number of times animal is motionless for at least 1 min. and completely visible	BORIS
Freeze - duration	Total duration of time animal is motionless for at least 1 min. and completely visible	BORIS
Refuge use - duration	Animal is attempting to be concealed within moss or in cup	BORIS
Activity - bouts	Total number of times animal is in motion, either stationary or through space	BORIS
Activity - speed	Average speed animal has moved through space	Annolid
Activity - distance	Total distance animal has moved	Annolid
Space use	Total area of space that the animal used (cm2)	Annolid
Male-female dyad
Male-female proximity	The average distance between animals	Annolid
Male-female joint refuge use	Both animals are refuging together in moss or cup	BORIS
Coordinated freezing	Both animals are freezing at same time	BORIS

Table 1: Description of coded behaviors

Analyses

Males selectively state match their partners for corticosterone irrespective of the endurance or lifetime reproductive output of partnerships

Prior to analysis, corticosterone concentration was log-transformed to achieve normal distribution when necessary. For all models, we conducted standard model diagnostics and assessed for outliers and deviations across quantiles using the ‘DHARMa’ package Hartig (2022). First, we tested corticosterone state matching of male-female dyads within and across conditions. We ran analyses separately for each condition (housing tank, experimental pre-treatment baseline, and post- stress/non-stress treatment) and dyad type (partner and non-partner). Since corticosterone concentration for the housing condition was normally distributed, we performed a simple linear correlation using Pearson’s correlation coefficient. For the experimental conditions, we performed a linear mixed model LMM using the package ‘lme4’, with male corticosterone log-concentration as the response variable, female corticosterone log-concentration and post-treatment condition (post-stress or post-non-stress) as fixed predictors, and frog identity and trial order as random factors. In brief, male corticosterone levels positively correlated with those of female partners within the housing baseline (Pearson’s correlation test: t = 2.69, r(8) = 0.69, P = 0.02; Figure 2 A) and experimental pre-treatment baseline conditions (LMM: β = 0.5, t = 5.2, R2 = 0.642, P < 0.001; Figure 2 B). There was a similar trend after controlling for the post-treatment conditions (stress/non-stress), although to a statistically non-significant extent (β = 0.001, t = 2.33, R2 = 0.83, P = 0.06; Figure 2 C), likely due to small sample size (n = 5 and 6, respectively) and increased response variation (see Supplementary Table 2 and Table 3 for details). Similarly, complementary research in voles has also discovered large effect sized, despite limited sample size (Burkett et al. 2016). By contrast, male corticosterone levels did not correlate with those of female non-partners across any condition (pre-treatment: LMM: β = 0.2, t = 0.72, R2 = 0.041, P = 0.48; post stress/non-stress: LMM: β = -0.0007, t = -0.28, R2 = 0.74, P = 0.78; Fig. 2D, E; see Supplementary Table 2 and Table 3 for details).

##################
#### PARTNERS ####
##################

#Correlation test in CORT levels between males and females in housing condition
House_corr <- Empathy %>% 
  filter(Condition2 == "housing_tank") %>% 
  pivot_wider(names_from = "Sex", values_from = "CORT_conc", id_cols = "CoupleID") %>% 
  do(tidy(cor.test(.$Male,.$Female)))

Figure2A <- Empathy %>% 
  filter(Condition2 == "housing_tank") %>% 
  pivot_wider(names_from = "Sex", values_from = "CORT_conc", id_cols = "CoupleID") %>% 
  ggplot(aes(Female, Male)) + my_geom_point + 
  my_geom_smooth + my_theme + 
  xlab("Females' CORT levels (arbitrary units)") +
  ylab("Males' CORT levels (arbitrary units)") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank())

#LMM to test relationship in CORT levels between males and females partners in baseline condition
PARTcorr_Bl <- CORTcorr %>%
  filter(MFamiliarity == "partner", MTrial == "Baseline") %>% 
  mutate(order.trial = c(1,2,1,2,1,2,1,1,2,1,1,2,1,2,1,2,1))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

#Model including the order of the trial
BlPARTNERSMod.ord <- lmer(log(MCORT_conc)~log(FCORT_conc) + (1|order.trial), 
                      data=PARTcorr_Bl) 
R2_Bl_partners <- r.squaredGLMM(BlPARTNERSMod.ord)
#ID as random factor led to singular fit. Since we are not including conditions (stress/non-stress), we exclude ID

#Check model fit
disp.test.1 <- function() {testDispersion(BlPARTNERSMod.ord)}
simulation.1 <- simulateResiduals(fittedModel = BlPARTNERSMod.ord, plot = F)
test.plots.1 <- function() {plot(simulation.1)}

Figure2B <- ggplot(PARTcorr_Bl, aes(log(FCORT_conc), log(MCORT_conc))) + 
  my_geom_point + my_geom_smooth + my_theme +
  xlab("Females' CORT levels (arbitrary units)") +
  ylab("Males' CORT levels (arbitrary units)")

#LMM to test relationship in CORT levels between males and females partners in treatment condition (stress/non-stress)
PARTcorr_Tr <- CORTcorr %>% 
  filter(MFamiliarity == "partner", MTrial == "Treatment") %>% 
  mutate(order.trial = c(1,2,1,1,2,1,2,1,2,1,1))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

#Model including the order of the trial
TrPARTNERSMod.ord <- lmer(log(MCORT_conc) ~ FCORT_conc + MCondition + 
                          (1|order.trial) + (1|MCoupleID), data = PARTcorr_Tr)

R2_Tr_partners <- r.squaredGLMM(TrPARTNERSMod.ord)

#Check model fit
disp.test.2 <- function() {testDispersion(TrPARTNERSMod.ord)}
simulation.2 <- simulateResiduals(fittedModel = TrPARTNERSMod.ord, plot = F)
test.plots.2 <- function() {plot(simulation.2)}

Figure2C <- ggplot(PARTcorr_Tr, aes(log(FCORT_conc),log(MCORT_conc))) + 
  geom_point(aes(fill=MCondition), shape=21, size=3.5) + 
  scale_fill_manual(values = c("black", "white")) +
  my_theme + 
  geom_smooth(aes(group = MCondition), method="lm", col="black", lwd=0.5, se=F) +
  xlab("Females' CORT levels (arbitrary units)") +
  ylab("Males' CORT levels (arbitrary units)")

######################
#### NON-PARTNERS ####
######################

#LMM to test relationship in CORT levels between males and females non-partners in baseline condition
NoPARTcorr_Bl <- CORTcorr %>% 
  filter(MFamiliarity == "non_partner", MTrial == "Baseline") %>% 
  mutate(order.trial = c(1,2,1,2,1,2,1,2,1,2,1,2,1))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

#Model including order of the test
BlNONPARTNERSMod.ord <- lmer(log(MCORT_conc)~log(FCORT_conc) + (1|order.trial), data=NoPARTcorr_Bl)
R2_Bl_nonpartners <- r.squaredGLMM(BlNONPARTNERSMod.ord)
#ID as random factor led to singular fit. Since we are not including conditions (stress/non-stress), we exclude ID

#Check model fit
disp.test.3 <- function() {testDispersion(BlNONPARTNERSMod.ord)}
simulation.3 <- simulateResiduals(fittedModel = BlNONPARTNERSMod.ord, plot = F)
test.plots.3 <- function() {plot(simulation.3)}

Figure2E <- ggplot(NoPARTcorr_Bl, aes(log(FCORT_conc), log(MCORT_conc))) + 
  my_geom_point + my_theme + 
  geom_smooth(method="lm", se=F, col="black", lwd=0.5) +
  xlab("Females' CORT levels (arbitrary units)") +
  ylab("Males' CORT levels (arbitrary units)")


#LMM to test relationship in CORT levels between males and females non-partners in treatment condition (stress/non-stress)
NoPARTcorr_Tr <- CORTcorr %>% 
  filter(MFamiliarity == "non_partner", MTrial == "Treatment") %>% 
  mutate(order.trial = c(1,2,1,2,1,1,1,2,1,2,2,1,1,2,1))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

#Model without including order of test
TrNONPARTNERSMod.ord <- lmer(log(MCORT_conc)~FCORT_conc + MCondition + 
                           + (1|order.trial) + (1|MCoupleID), data=NoPARTcorr_Tr)
R2_Tr_nonpartners <- r.squaredGLMM(TrNONPARTNERSMod.ord)

#Check model fit
disp.test.4 <- function() {testDispersion(TrNONPARTNERSMod.ord)}
simulation.4 <- simulateResiduals(fittedModel = TrNONPARTNERSMod.ord, plot = F)
test.plots.4 <- function() {plot(simulation.4)}

Figure2F <- ggplot(NoPARTcorr_Tr, aes(log(FCORT_conc), log(MCORT_conc))) + 
  geom_point(aes(fill=MCondition), shape=21, size=3.5) + 
  scale_fill_manual(values = c("black", "white")) +
  my_theme + 
  geom_smooth(aes(group = MCondition), method="lm", col="black", lwd=0.5, se=F) +
  xlab("Females' CORT levels (arbitrary units)") +
  ylab("Males' CORT levels (arbitrary units)")

Figure 2: Partner-selective corticosterone matching in pair bonded R. imitator. Male corticosterone (CORT) level matched (A-C) female partners, but not (E-F) familiar female non-partners, across semi-natural and experimental conditions. Regression lines are shown, with shaded regions denoting 95% confidence intervals for statistically significant results. (D) Example of gentle leg-restraint stress treatment.

## Supplementary table 1: Summarized results of linear mixed-effect models showing the relationship in post-treatment corticosterone levels between male-female partners and non-partners

tab_model(TrPARTNERSMod.ord, TrNONPARTNERSMod.ord, digits.p = 3,
                     pred.labels = c("(Intercept)","Female CORT levels", 
                                     "Condition (Stress)"), dv.labels = c("Male CORT levels (partners)",
                                                                          "Male CORT levels (non-partners)"))

	Male CORT levels (partners)			Male CORT levels (non-partners)
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	4.75	3.30 – 6.21	<0.001	5.16	3.01 – 7.30	<0.001
Female CORT levels	0.00	-0.00 – 0.00	0.067	-0.00	-0.01 – 0.00	0.783
Condition (Stress)	-0.81	-1.95 – 0.34	0.129	-2.03	-3.32 – -0.74	0.006
Random Effects
σ2	0.36			0.85
τ00	1.21 MCoupleID			0.93 MCoupleID
	0.00 order.trial			0.54 order.trial
ICC				0.63
N	2 order.trial			2 order.trial
	7 MCoupleID			9 MCoupleID
Observations	11			15
Marginal R2 / Conditional R2	0.629 / NA			0.303 / 0.744

Table 2: Supplementary Table 2. Summarized results of linear mixed-effect models showing the relationship in post-treatment corticosterone levels between male-female partners and non-partners

## Supplementary table 3: Summarized results of linear mixed-effect models showing the relationship in pre-treatment corticosterone levels between male-female partners and non-partners

tab_model(BlPARTNERSMod.ord, BlNONPARTNERSMod.ord, digits.p = 3,
                     pred.labels = c("(Intercept)","Female CORT levels"), dv.labels = c("Male CORT levels (partners)", "Male CORT levels (non-partners)"))

	Male CORT levels (partners)			Male CORT levels (non-partners)
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	2.70	1.64 – 3.75	<0.001	3.62	0.57 – 6.67	0.025
Female CORT levels	0.50	0.30 – 0.71	<0.001	0.21	-0.45 – 0.86	0.490
Random Effects
σ2	0.37			1.02
τ00	0.02 order.trial			0.00 order.trial
ICC	0.06
N	2 order.trial			2 order.trial
Observations	17			13
Marginal R2 / Conditional R2	0.615 / 0.637			0.041 / NA

Table 3: Supplementary Table 3. Summarized results of linear mixed-effect models showing the relationship in pre-treatment corticosterone levels between male-female partners and non-partners

Next, we tested whether male-female state matching predicted the longevity or lifetime reproductive output of partnerships. We ran a separate multiple regression for each condition (housing, experimental pre-treatment, post- stress, and post-non-stress treatments), using male corticosterone levels as the response, and female corticosterone levels, partnership endurance, reproductive output and trial order as predictors. The extent of male-to-female partner corticosterone matching was not explained by the endurance or lifetime reproductive output of partnerships (P-values for all treatment conditions ≥ 0.05; Figure 3; see Supplementary Table 4 for details).

# Housing tanks

House_mod <- Empathy %>% 
  filter(Condition2 == "housing_tank") %>% 
  pivot_wider(names_from = "Sex", values_from = "CORT_conc", 
              id_cols = c("CoupleID", "part_endurance", "reprod")) %>% 
  mutate(male_CORT = log(Male)) %>% 
  mutate(female_CORT = log(Female)) %>% 
  do(tidy(lm(.$male_CORT ~ .$female_CORT + .$part_endurance + .$reprod))) %>% 
  mutate(term = c("Intercept", "female_CORT", "PartnerEndur", "N_offspring")) %>% 
  column_to_rownames("term") %>% 
  select(estimate, statistic, p.value) %>% 
  rename("ß" = estimate, "t" = statistic, "p" = p.value)

# partners-baseline

Baseline_mod <- CORTcorr %>% 
  filter(MTrial == "Baseline" & MFamiliarity == "partner") %>% 
  mutate(male_CORT = log(MCORT_conc)) %>% 
  mutate(female_CORT = log(FCORT_conc)) %>% 
  do(tidy(lm(.$male_CORT ~ .$female_CORT + .$part_endurance + .$reprod))) %>% 
  mutate(term = c("Intercept", "female_CORT", "PartnerEndur", "N_offspring")) %>% 
  column_to_rownames("term") %>% 
  select(estimate, statistic, p.value) %>% 
  rename("ß" = estimate, "t" = statistic, "p" = p.value)

# partners-treatment-stress

Stress_mod <- CORTcorr %>% 
  filter(MTrial == "Treatment" & MFamiliarity == "partner" & 
           MCondition == "Stress") %>% 
  mutate(male_CORT = log(MCORT_conc)) %>% 
  mutate(female_CORT = log(FCORT_conc)) %>% 
  do(tidy(lm(.$male_CORT ~ .$female_CORT + .$part_endurance + .$reprod))) %>% 
  mutate(term = c("Intercept", "female_CORT", "PartnerEndur", "N_offspring")) %>% 
  column_to_rownames("term") %>% 
  select(estimate, statistic, p.value) %>% 
  rename("ß" = estimate, "t" = statistic, "p" = p.value)

# partners-treatment-non_stress

NStress_mod <- CORTcorr %>% 
  filter(MTrial == "Treatment" & MFamiliarity == "partner" & 
           MCondition == "non_stress") %>% 
  mutate(male_CORT = log(MCORT_conc)) %>% 
  mutate(female_CORT = log(FCORT_conc)) %>% 
  do(tidy(lm(.$male_CORT ~ .$female_CORT + .$part_endurance + .$reprod))) %>% 
  mutate(term = c("Intercept", "female_CORT", "PartnerEndur", "N_offspring")) %>% 
  column_to_rownames("term") %>% 
  select(estimate, statistic, p.value) %>% 
  rename("ß" = estimate, "t" = statistic, "p" = p.value)

# All models together 

House_mod %>% 
  cbind(Baseline_mod, Stress_mod, NStress_mod) %>% 
  kable(digits = 2, 
        table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T) %>% 
  add_header_above(c(" ", "Housing" = 3, "Baseline" = 3, "Stress" = 3,
                     "non-stress" = 3))
Figure3 <- Empathy %>% 
  filter(Condition2 == "housing_tank") %>% 
  pivot_wider(names_from = "Sex", values_from = "CORT_conc", 
              id_cols = c("CoupleID", "part_endurance", "reprod")) %>% 
  mutate(male_CORT = log(Male)) %>% 
  mutate(female_CORT = log(Female)) %>%
  mutate(condition = c(rep("Housing",10))) %>% 
  bind_rows((CORTcorr %>% 
  filter(MTrial == "Baseline" & MFamiliarity == "partner") %>% 
  mutate(male_CORT = log(MCORT_conc)) %>% 
  mutate(female_CORT = log(FCORT_conc)) %>% 
    select(CoupleID = MCoupleID, part_endurance, reprod, Male = MCORT_conc,
           Female = FCORT_conc, male_CORT, female_CORT) %>% 
    mutate(condition = c(rep("baseline",17)))),
  (CORTcorr %>% 
  filter(MTrial == "Treatment" & MFamiliarity == "partner" & 
           MCondition == "Stress") %>% 
  mutate(male_CORT = log(MCORT_conc)) %>% 
  mutate(female_CORT = log(FCORT_conc)) %>% 
    select(CoupleID = MCoupleID, part_endurance, reprod, Male = MCORT_conc,
           Female = FCORT_conc, male_CORT, female_CORT) %>% 
    mutate(condition = c(rep("stress",6)))),
  (CORTcorr %>% 
  filter(MTrial == "Treatment" & MFamiliarity == "partner" & 
           MCondition == "non_stress") %>% 
  mutate(male_CORT = log(MCORT_conc)) %>% 
  mutate(female_CORT = log(FCORT_conc)) %>% 
    select(CoupleID = MCoupleID, part_endurance, reprod, Male = MCORT_conc,
           Female = FCORT_conc, male_CORT, female_CORT) %>% 
    mutate(condition = c(rep("non_stress",5))))) %>% 
  ggplot(aes(x=female_CORT, y=male_CORT)) + 
  geom_smooth(method="lm", col="gray55", lty=2, lwd=0.7) +
  geom_point(aes(colour=part_endurance, size=reprod)) + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(size = 15)) +
  xlab("Females' CORT levels\n(arbitrary units)") +
  ylab("Males' CORT levels\n(arbitrary units)") +
  scale_colour_viridis() +
  facet_wrap(.~condition, scales = "free")

	Housing			Baseline			Stress			non-stress
	ß	t	p	ß	t	p	ß	t	p	ß	t	p
Intercept	2.83	1.27	0.25	3.22	5.54	0.00	0.02	0.01	0.99	-0.76	-2.52	0.24
female_CORT	0.44	0.99	0.36	0.53	5.37	0.00	1.11	3.18	0.09	1.37	21.97	0.03
PartnerEndur	0.00	1.37	0.22	0.00	-1.34	0.20	0.00	0.51	0.66	0.04	18.15	0.04
N_offspring	0.00	-0.96	0.37	0.01	0.85	0.41	-0.08	-1.97	0.19	-0.45	-18.32	0.03

Table 4: Supplementary Table 4. Summarized results of multiple regressions of emotional contagion between male-female partners, partnership endurance, and total number of offspring across housing and experimental conditions

Figure 3: Male-female corticosterone state matching is irrespective of the endurance or lifetime reproductive output of partnerships. Shown are multiple regression plots of male-female corticosterone relationship across conditions after controlling for the endurance and total number of offspring of partnerships. Shaded area around multiple regression lines denote 95% confidence intervals.

Stress contagion assay does not provoke a stress response, nor do males behaviorally state match females

We tested for differences in corticosterone concentration between housing and experimental pre-treatment (baseline) condition. For this, we performed a linear mixed-effects model (LMM). We included log-transformed corticosterone concentration as the response variable, and the interaction between condition (housing and experimental pre-treatment baseline), dyad type (partner and non-partner), sex, and trial order as predictors, and the frog identification as a random factor. We further computed estimated marginal means (least-squares means) contrasts using the ‘emmeans’ function within the emmeans package (Lenth 2023). P-values were adjusted using Tukey’s method. Only pair bonded dyads were examined in the housing condition since animals were housed only with partners.
Differences in corticosterone concentration between experimental treatment conditions were examined separately for each dyad type and sex also with LMMs. The models included log-transformed corticosterone concentration as the response variable, the interaction between treatment condition (non-stress and stress) and trial (pre- and post-treatment) as fixed effects, and frog identity and trial order as random factors. Post-hoc analyses were conducted using estimated marginal means with Tukey’s adjustment method.

Subjecting females to the stress treatment caused no increase in corticosterone levels in female demonstrators or male observers across time points, irrespective of whether they were normalized to baseline (planned least-squares means contrast: P ≥ 0.05 in all groups; Supplementary Table 5, Table 6 and Supplementary Figure 1 A, B, D for details).

Figure 4: Supplementary Figure 1. Subjecting females to a leg restraint treatment did not cause a corticosterone (CORT) or behavioral stress response in either female demonstrators or male observers in partnered or non-partner dyads. (A) Water-borne CORT concentration does not differ between conditions (housing vs. experimental pre-treatment (baseline)), dyad (partner vs. non-partner), or sex. Dots and horizontal bars represent the mean, and grey boxes represent 95% confidence intervals. (B, D) Neither water-borne CORT concentration nor (C, E) behavior differs between experimental conditions (pre-treatment baseline vs. post non-stress/stress treatment), in either dyad or sex. Nine behavioral types are represented by 3 non-redundant principal component categories: PC1 = male-female interaction, PC2 = activity, PC3 = stillness. Behavioral values shown are relative to the experimental pre-treatment (baseline) condition. Panels B-E: dots represent the mean and bars represent 95% confidence intervals for the least-squares means. NS = non-significant (P ≥ 0.05)

#####################
# partners - males #
#####################

male_partners <- Empathy %>% 
  filter(Familiarity == "partner" & Sex == "Male") %>% 
  filter(!Condition2 == "housing_tank") %>% 
  mutate(Date2c = dmy(Date2, tz = "UTC")) %>%  # Explicitly parsing dates with a timezone can sometimes help
  arrange(ID) %>% 
  mutate(order.trial =c(1,1,2,1,1,2,2,2,2,1,1,1,2,2,1,1,2,2,2,2,1,
                        2,2,1,1,2,2,1,1,1,1,2,2))#This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

# Including the order of the trial "order.trial" as random factor
malePARTNERSmod.ord <- lmer(log(CORT_conc) ~ Condition * Trial + (1 | ID) + (1 | order.trial), data=male_partners)
MALE_pairW.ord <- emmeans(malePARTNERSmod.ord, pairwise~Condition * Trial, adjust="tukey")
emmeans1 <- plot(MALE_pairW.ord, comparisons = TRUE, col=c('#0979A1', 'black')) + 
  my_theme + coord_flip()

# Check model fit
disp.test.5 <- function() {testDispersion(malePARTNERSmod.ord)}
simulation.5 <- simulateResiduals(fittedModel = malePARTNERSmod.ord, plot = F)
test.plots.5 <- function() {plot(simulation.5)}

# emmeans plot
M_emme = plot(MALE_pairW.ord, comparisons = TRUE, col=c('#0979A1', 'black')) + 
  my_theme + coord_flip() + xlim(2,8)

######################
# partners - females #
######################

fem_partners <- Empathy %>% 
  filter(Familiarity == "partner" & Sex == "Female") %>% 
  filter(!Condition2 == "housing_tank") %>% 
  mutate(Date2c = dmy(Date2, tz = "UTC")) %>%  # Explicitly parsing dates with a timezone can sometimes help
  arrange(ID) %>% 
  mutate(order.trial = c(2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 1))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

# Including order of the trial as random factor

femalePARTNERSmod.ord <- lmer(log(CORT_conc) ~ Condition*Trial + (1 | ID) + (1 | order.trial), data=fem_partners)
FEMALE_pairW.ord <- emmeans(femalePARTNERSmod.ord, pairwise~Condition*Trial, 
                        adjust="tukey")

# Check model fit
disp.test.6 <- function() {testDispersion(femalePARTNERSmod.ord)}
simulation.6 <- simulateResiduals(fittedModel = femalePARTNERSmod.ord, plot = F)
test.plots.6 <- function() {plot(simulation.6)}

# emmeans plot
F_emme <- plot(FEMALE_pairW.ord, comparisons = TRUE, col=c('#DD4E25', 'black')) +
  my_theme + coord_flip() + xlim(2,8)


########################
# non_partners - males #
########################

male_nonpartners <- Empathy %>% 
  filter(Familiarity == "non_partner" & Sex == "Male") %>% 
  filter(!Condition2 == "housing_tank") %>% 
  mutate(Date2c = dmy(Date2, tz = "UTC")) %>%  # Explicitly parsing dates with a timezone can sometimes help
  arrange(ID) %>% 
  mutate(order.trial = c(1, 2, 2, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 2))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

# Including the order of the trial as random factor  
maleNONPARTNERSmod.ord <- lmer(log(CORT_conc) ~ Condition*Trial + 
                              (1 | ID) + (1 | order.trial), data=male_nonpartners)

MALE_nonpairW.ord <- emmeans(maleNONPARTNERSmod.ord, pairwise~Condition*Trial, 
                         adjust="tukey")

# Check model fit
disp.test.7 <- function() {testDispersion(maleNONPARTNERSmod.ord)}
simulation.7 <- simulateResiduals(fittedModel = maleNONPARTNERSmod.ord, plot = F)
test.plots.7 <- function() {plot(simulation.7)}

# emmeans plot
M_emme2 <- plot(MALE_nonpairW.ord, comparisons = TRUE, col=c('#0979A1', 'black')) + 
  my_theme + coord_flip() + xlim(2,6.5)


###########################
# non_partners - females #
##########################

fem_nonpartners <- Empathy %>% 
  filter(Familiarity == "non_partner" & Sex == "Female") %>% 
  filter(!Condition2 == "housing_tank") %>% 
  mutate(Date2c = dmy(Date2, tz = "UTC")) %>%  # Explicitly parsing dates with a timezone can sometimes help
  arrange(ID) %>% 
  mutate(order.trial = c(1,2,1,2,1,2,2,1,2,1,1,2,1,2,2,1,1,2,1,1,2,2,1,1,2,1,2,2,1,1,2,2))# This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

# Including the order of the trial as random factor  
femaleNONPARTNERSmod.ord <- lmer(log(CORT_conc) ~ Condition*Trial + (1 | order.trial) + (1 | ID), data=fem_nonpartners)
FEMALE_nonpairW.ord <- emmeans(femaleNONPARTNERSmod.ord, pairwise~Condition*Trial, 
                           adjust="tukey")

# Check model fit
disp.test.8 <- function() {testDispersion(femaleNONPARTNERSmod.ord)}
simulation.8 <- simulateResiduals(fittedModel = femaleNONPARTNERSmod.ord, plot = F)
test.plots.8 <- function() {plot(simulation.8)}

# emmeans plot
F_emme2 <- plot(FEMALE_nonpairW.ord, comparisons = TRUE, 
                col=c('#DD4E25', 'black')) +
  my_theme + coord_flip() + xlim(2,6.5)



###########################################
#### emmeans table for all conditions #####
###########################################

data.frame(MALE_pairW.ord$contrasts) %>% 
  bind_rows(data.frame(FEMALE_pairW.ord$contrasts)) %>% 
  bind_rows(data.frame(MALE_nonpairW.ord$contrasts)) %>% 
  bind_rows(data.frame(FEMALE_nonpairW.ord$contrasts)) %>% 
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T) %>%
  pack_rows("Males: partner dyad", 1, 6) %>%
  pack_rows("Females: partner dyad", 7, 12) %>% 
  pack_rows("Males: non-partner dyad", 13, 18) %>%
  pack_rows("Females: non-partner dyad", 19, 24)

contrast	estimate	SE	df	t.ratio	p.value
Males: partner dyad
non_stress Baseline - Stress Baseline	-0.34	0.30	20.24	-1.13	0.68
non_stress Baseline - non_stress Treatment	0.05	0.32	20.26	0.14	1.00
non_stress Baseline - Stress Treatment	0.49	0.31	20.21	1.61	0.39
Stress Baseline - non_stress Treatment	0.38	0.33	20.63	1.14	0.67
Stress Baseline - Stress Treatment	0.83	0.30	20.12	2.74	0.06
non_stress Treatment - Stress Treatment	0.45	0.34	20.62	1.31	0.56
Females: partner dyad
non_stress Baseline - Stress Baseline	0.00	0.73	17.88	0.00	1.00
non_stress Baseline - non_stress Treatment	0.80	0.74	17.85	1.08	0.70
non_stress Baseline - Stress Treatment	-0.91	0.79	18.64	-1.16	0.66
Stress Baseline - non_stress Treatment	0.80	0.79	18.13	1.01	0.74
Stress Baseline - Stress Treatment	-0.91	0.81	17.26	-1.12	0.68
non_stress Treatment - Stress Treatment	-1.71	0.83	17.84	-2.05	0.21
Males: non-partner dyad
non_stress Baseline - Stress Baseline	-0.01	0.56	18.55	-0.01	1.00
non_stress Baseline - non_stress Treatment	-0.48	0.61	18.63	-0.79	0.86
non_stress Baseline - Stress Treatment	1.34	0.55	18.25	2.45	0.10
Stress Baseline - non_stress Treatment	-0.48	0.61	19.38	-0.78	0.86
Stress Baseline - Stress Treatment	1.35	0.52	17.55	2.61	0.08
non_stress Treatment - Stress Treatment	1.83	0.60	19.11	3.06	0.03
Females: non-partner dyad
non_stress Baseline - Stress Baseline	0.10	0.59	19.54	0.16	1.00
non_stress Baseline - non_stress Treatment	0.25	0.53	19.43	0.48	0.96
non_stress Baseline - Stress Treatment	0.32	0.53	19.43	0.61	0.93
Stress Baseline - non_stress Treatment	0.16	0.59	19.88	0.27	0.99
Stress Baseline - Stress Treatment	0.23	0.59	19.88	0.39	0.98
non_stress Treatment - Stress Treatment	0.07	0.51	19.16	0.13	1.00

Table 5: Supplementary Table 5. Summarized results of planned least-squares means contrast of corticosterone concentration across experimental treatment conditions. P-values were adjusted using Tukey’s method.

# Baseline comparisons

Treatments <- Empathy %>% 
  filter(Trial == "Baseline") %>% 
  mutate(Condition2 = c(rep("Experiment", 63), rep("Housing", 20))) %>% 
  dplyr::select(Condition2, Familiarity, Sex, CORT_conc, ID, Date2, Condition) %>% 
  mutate(Condition3 = paste(Condition2, Familiarity, Sex, sep = "_")) %>% 
  arrange(ID, Familiarity) %>% 
  mutate(order.trial = c(1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 
                         2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 
                         1, 1, 1, 2, 2, 1, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 
                         1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                         1, 1, 1, 1, 1, 1, 1, 1, 1, 1)) # This is the order of the experiments according to the date. 1 means that the corresponding condition (stress / non-stress) was tested first and 2 that was tested second.

SF1A <- Treatments %>% 
  ggplot(aes(Condition3, log(CORT_conc))) +
  geom_jitter(position=position_jitter(0.05),shape=21, colour="grey", 
              aes(fill=Condition3), alpha=0.3, size=3) +
  stat_summary(fun.data="mean_cl_boot", geom="crossbar", 
               aes(fill=Condition3),shape=21, size=.3, alpha=0.3, 
               col="gray45") + 
  stat_summary(fun.y="mean", geom="point", col="gray45", fill="grey",
               shape=21, size=1) +
  scale_fill_viridis_d() +
  my_theme +
  theme(axis.text.x = element_text(angle = 60))

#Including the order of the trials
comp_treat <- lmer(log(CORT_conc) ~ Condition3 + order.trial + (1 | ID), data=Treatments)

comp_treatpairW <- emmeans(comp_treat, pairwise~Condition3, 
                           adjust="tukey")

# Check model fit
disp.test.9 <- function() {testDispersion(comp_treat)}
simulation.9 <- simulateResiduals(fittedModel = comp_treat, plot = F)
test.plots.9 <- function() {plot(simulation.9)}

comp_treatpairW$contrasts %>% 
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T)

contrast	estimate	SE	df	t.ratio	p.value
Experiment_non_partner_Female - Experiment_non_partner_Male	0.27	0.41	40.44	0.67	0.98
Experiment_non_partner_Female - Experiment_partner_Female	-0.15	0.28	50.84	-0.52	1.00
Experiment_non_partner_Female - Experiment_partner_Male	-0.33	0.40	37.75	-0.82	0.96
Experiment_non_partner_Female - Housing_partner_Female	-0.11	0.43	52.14	-0.26	1.00
Experiment_non_partner_Female - Housing_partner_Male	-0.22	0.43	52.14	-0.52	1.00
Experiment_non_partner_Male - Experiment_partner_Female	-0.42	0.40	38.31	-1.05	0.90
Experiment_non_partner_Male - Experiment_partner_Male	-0.61	0.26	43.74	-2.37	0.19
Experiment_non_partner_Male - Housing_partner_Female	-0.39	0.43	48.71	-0.90	0.94
Experiment_non_partner_Male - Housing_partner_Male	-0.50	0.43	48.71	-1.16	0.85
Experiment_partner_Female - Experiment_partner_Male	-0.18	0.39	35.52	-0.47	1.00
Experiment_partner_Female - Housing_partner_Female	0.04	0.42	50.45	0.09	1.00
Experiment_partner_Female - Housing_partner_Male	-0.07	0.42	50.45	-0.18	1.00
Experiment_partner_Male - Housing_partner_Female	0.22	0.42	47.46	0.52	1.00
Experiment_partner_Male - Housing_partner_Male	0.11	0.42	47.46	0.26	1.00
Housing_partner_Female - Housing_partner_Male	-0.11	0.43	55.56	-0.26	1.00

Table 6: Supplementary Table 6. Summarized results of planned least-squares means contrast between housing and experimental conditions in males’ and females’ corticosterone concentration, within ‘partners’ and ‘non-partners’ condition. P-values were adjusted using Tukey’s method. Results are averaged over the order or the trials.

To account for and eliminate intra-individual sources of variation in behavior, we subtracted the pre-treatment baseline from the post-stress/non-stress treatment conditions. Then, we extracted behaviors, namely activity level, freezing duration, partner approach duration, partner gaze duration, total distance moved, and total area moved. Next, we minimized redundancy of the extracted behaviors using a varimax normalized principal components analysis (PCA) for each dyad type, using the function ‘principal’ within the psych package (2023). Principal component analyses for each dyad type (partners and non-partners) yielded three components that explained more than 70% of the total variance (See Supplementary Table 7).

## keep first order
order.trial.par <- bhs[-118,] %>% 
  dplyr::select(Familiarity, ord.trial, tot_distance_moved_cm) %>% #take total distance as reference for removing NAs
  na.omit() %>% 
  filter(row_number() %% 2 == 0) %>% 
  filter(Familiarity == "Partner") %>% 
  dplyr::select(ord.trial)

order.trial.npar <- bhs[-118,] %>% 
  dplyr::select(Familiarity, ord.trial, tot_distance_moved_cm) %>% #take total distance as reference for removing NAs
  na.omit() %>% 
  filter(row_number() %% 2 == 0) %>% 
  filter(Familiarity == "Non-Partner") %>%
  dplyr::select(ord.trial)

## to select just the individual behaviours 
bhs.1 <- na.omit(bhs[,c(1:7,10,11,13,14,15,16,17,18,19,20)])

## Subset for the categorical variables that will be used later for modeling
bhs_chrs <- bhs.1 %>% 
  filter(State == "Pre") %>% 
  dplyr::select(Familiarity:Subject)

## Subset for baseline
bhs_pre <- bhs.1 %>% 
  filter(State == "Pre") # %>% 
  #select(Subject:Total_area_used_cm2) %>%
  #arrange(desc(Subject)) %>% 
  #mutate_if(is.character,as.numeric) %>% 
  #select(Activity_level:Total_area_used_cm2)

## Subset for treatment
bhs_post <- bhs.1 %>% 
  filter(State == "Post") #%>% 
  #select(Subject:Total_area_used_cm2) %>%
  #arrange(desc(Subject)) %>% 
  #mutate_if(is.character,as.numeric) %>% 
  #select(Activity_level:Total_area_used_cm2)

## Subtracted data
bhs2 <- bhs_post[-56,7:17]-bhs_pre[,7:17]#Female 123 subtracted because It has NAs in movement behaviors in pre condition

## Merge with categorical variables to create subset for every familiarity
bhs2.1 <- cbind(bhs_chrs,bhs2)

## PARTNERS (I removed again the categorical data for the correlations)
p_bhs <- bhs2.1 %>% 
  filter(Familiarity=="Partner") %>% 
  dplyr::select(Activity_level:Total_area_used_cm2)
  

## NON-PARTNERS
Np_bhs <- bhs2.1 %>% 
  filter(Familiarity=="Non-Partner") %>% 
  dplyr::select(Activity_level:Total_area_used_cm2)

## Then, we checked for redundancy among behaviours or mathematically related variables.
## For this, we calculated the variance inflation factors (VIF), where VIFs higher than 5 revealed severe 
## multicollinearity among individual behaviours:

## Correlation matrix
cors_bhs.1 <- cor(p_bhs)

## Variance inflation factors. Here I created a random variable as dependent 
##just to include all behaviours as predictors
vifMOD1.1 <- lm(rnorm(31) ~ Activity_level + FreezingB +  
               PartnerApproachB + PartnerApproachD +
               PartnerGazeB + PartnerGazeD + 
               AV_speed_per_sec_cm +
               tot_distance_moved_cm +
               Total_area_used_cm2 + FreezingD + RefugeUse, data = p_bhs)

vifs1.1 <- vif(vifMOD1.1)

## We can see that 'FreezingB', 'PartnerApproachB', 'PartnerGazeB', 
## 'Speed' and 'Refuge use' were highly collinear and therefore, 
## we excluded from the analysis
vifMOD1.2 <- lm(rnorm(31) ~ Activity_level + 
                  PartnerApproachD + PartnerGazeD + 
                 tot_distance_moved_cm + Total_area_used_cm2 + 
                 FreezingD, data = p_bhs)

vifs1.2 <- vif(vifMOD1.2)

## Data after removing the variables above
p_bhs_2 <- p_bhs[,-c(2,4,6,8,9)]



##NON-PARTNERS

## Correlation matrix
cors_bhs.2 <- cor(Np_bhs)

## Variance inflation factors. Here I created a random variable as dependent 
##just to include all behaviours as predictors
vifMOD2.1 <- lm(rnorm(36) ~ Activity_level + FreezingB +  
               PartnerApproachB + PartnerApproachD +
               PartnerGazeB + PartnerGazeD + 
               AV_speed_per_sec_cm +
               tot_distance_moved_cm +
               Total_area_used_cm2 + FreezingD + RefugeUse, data = Np_bhs)

vifs2.1 <- vif(vifMOD2.1)

## We can see that 'FreezingB', 'PartnerApproachB', 'PartnerGazeB', 
## 'Speed' and 'Refuge use' were highly collinear and therefore, 
## we excluded from the analysis (Activity level was also high, but I decided to
## keep it and see VIF after removing the other variables)
vifMOD2.2 <- lm(rnorm(36) ~ Activity_level + 
                  PartnerApproachD + PartnerGazeD + 
                 tot_distance_moved_cm + Total_area_used_cm2 + 
                 FreezingD, data = Np_bhs)

vifs2.2 <- vif(vifMOD2.2) ## Now Activity level is o.k.

## Data after removing the variables above
Np_bhs_2 <- Np_bhs[,-c(2,4,6,8,9)]

## Then, we checked if the variables are orthogonal (i.e., correlated enough that the PCA makes sense) using a correlation matrix and the Bartlett's test (H0 = variables are not correlated):

## PARTNERS
c.btest.1 <- cortest.bartlett(p_bhs_2, n=31)

cors_bhs.par <- cor(p_bhs_2)

## NON-PARTNERS
c.btest.2 <- cortest.bartlett(Np_bhs_2, n=36)

cors_bhs.Npar <- cor(Np_bhs_2)

######################################
####### VARIMAX Rotated PCA's ########
######################################

## PARTNERS

## To calculate the eigen values: 2 higher than 1, but third one close to 1. 
## So I keep it in 3
eigen.1 <- eigen(cors_bhs.par)

## PCA with three factors according to eigen values
pca.par <- principal(p_bhs_2, nfactors = 3, covar=FALSE, scores=TRUE,
                  rotate="varimax")

## To extract the loadings 
Loadings.1 <- pca.par$loadings[,1:3] %>% 
  data.frame() %>% 
  rename("PC1 (31%)" = RC1, "PC3 (22%)" = RC3, "PC2 (21%)" = RC2)

scores.1 <- as.data.frame(pca.par$scores) #to check the loads of every original variable in every component.


## NON-PARTNERS

## To calculate the eigen values: 3 higher than 1. 

eigen.2 <- eigen(cors_bhs.Npar)

## PCA with three factors according to eigen values
pca.Npar <- principal(Np_bhs_2,nfactors = 3, covar=FALSE, scores=TRUE,
                  rotate="varimax")

## To extract the loadings 
Loadings.2 <- pca.Npar$loadings[,1:3] %>% 
  data.frame() %>% 
  rename("PC1 (34%)" = RC1, "PC2 (29%)" = RC2, "PC3 (20%)" = RC3)

scores.2 <- as.data.frame(pca.Npar$scores) #to check the loads of every original variable in every component.

#Supplementary table 7
Loadings.1 %>% 
  bind_cols(Loadings.2) %>% 
  rownames_to_column("Behavior") %>% 
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  add_header_above(c(" ", "Partners" = 3, "Non-partners" = 3), bold = T) %>%  
  row_spec(0, italic = T, bold = T) 

	Partners			Non-partners
Behavior	PC1 (31%)	PC3 (22%)	PC2 (21%)	PC1 (34%)	PC2 (29%)	PC3 (20%)
Activity_level	0.88	0.00	0.06	0.92	0.14	-0.04
FreezingD	-0.10	0.03	0.91	-0.31	-0.17	0.85
PartnerApproachD	0.30	0.79	-0.17	0.19	0.89	-0.23
PartnerGazeD	-0.12	0.75	0.34	0.06	0.93	0.10
tot_distance_moved_cm	0.57	0.34	0.51	0.56	0.13	0.65
Total_area_used_cm2	0.82	0.11	-0.22	0.87	0.11	-0.12

Table 7: Supplementary Table 7. Standardized loading matrix of behavioral variables represented in rotated principal components. Values in parentheses represent the percentage of explained variance for each component.

We tested whether subjecting females to the stress treatment affected the individual or pairwise behavior of female demonstrators and/or male observers, using simple linear models. For this, we performed each model with each of the rendered components as response variables, and experimental conditions (post-stress/-non-stress), sex, and trial order as predictors. We further computed estimated marginal means (least-squares means) contrasts among conditions and sex, using Tukey’s adjustment method. To determine whether males displayed behavioral indicators of partner-selective emotional contagion, we examined male-female state-matching for stress related individual behaviors within the rendered principal components. For this, we ran simple correlations for each of the four experimental condition*dyad type combinations.

Subjecting females to the stress treatment had no effect on the individual or coordinated stress behavior of female demonstrators or male observers of partner and non-partner dyads in any time point comparison, irrespective of whether responses were normalized to baseline (P ≥ 0.05 for all treatment comparisons; Supplementary Table 8, Table 9, Table 10, Table 11, and Supplementary Figure 1 C, E and Figure 2 D-G for details). Together, this indicates that the treatment did not provoke a stress response as intended.

# Now we merge the scores of the 3 components with the categorical variables and separated males' and females' behaviours:

## PARTNERS
## To merge the categorical variables with the scores and order of the trial
PCA_par <- bhs2.1 %>% 
  filter(Familiarity=="Partner") %>% 
  mutate(scores.1) %>% 
  mutate(order.trial = order.trial.par$ord.trial)


## Subset for males  
bhs.m.par <- PCA_par %>% 
  filter(Sex=="Male") %>% 
  dplyr::select(Familiarity:Sex, PC1.m=RC1, PC2.m=RC2, PC3.m=RC3, order.trial)

## Subset for females
bhs.f.par <- PCA_par %>% 
  filter(Sex=="Female") %>% 
  dplyr::select(PC1.f=RC1, PC2.f=RC2, PC3.f=RC3, order.trial) 

PCA_par.All <- bhs.m.par[-12,] %>% 
  mutate(bhs.f.par)


## NON-PARTNERS
## To merge the categorical variables with the scores 
PCA_Npar <- bhs2.1 %>% 
  filter(Familiarity=="Non-Partner") %>% 
  mutate(scores.2) %>% 
  mutate(order.trial = order.trial.npar$ord.trial)

## Subset for males  
bhs.m.Npar <- PCA_Npar %>% 
  filter(Sex=="Male") %>% 
  dplyr::select(Familiarity:Sex, PC1.m=RC1, PC2.m=RC2, PC3.m=RC3, order.trial)

## Subset for females
bhs.f.Npar <- PCA_Npar %>% 
  filter(Sex=="Female") %>% 
  dplyr::select(PC1.f=RC1, PC2.f=RC2, PC3.f=RC3, order.trial) 

PCA_Npar.All <- bhs.m.Npar %>% 
  mutate(bhs.f.Npar)

#####################
###### emmeans ######
#####################

## Partner
PCA_comp.par <- PCA_par.All %>% 
  pivot_longer(c("PC1.m", "PC2.m", "PC3.m", "PC1.f", "PC2.f", "PC3.f"),
               names_to = "PC.sex", values_to = "PC.beh") %>% 
  mutate(Sex = rep(c(rep("Male",3), rep("Female",3)),15)) %>% 
  mutate(PC = rep(c("PC1", "PC2", "PC3"),30))

PCA_comp.par$Sex <- factor(PCA_comp.par$Sex, levels=c("Male", "Female"))


#Linear models per PC followed by pairwise comparisons 

#PC1
PC1p <- PCA_comp.par %>% 
  filter(PC == "PC1")
PC1mod.par <- lm(PC.beh ~ Condition*Sex*order.trial, data=PC1p)

PC1.pairW.par <- emmeans(PC1mod.par, pairwise~Condition*Sex, 
                           adjust="tukey")$contrast %>% 
  data.frame()

#PC2
PC2p <- PCA_comp.par %>% 
  filter(PC == "PC2")
PC2mod.par <- lm(PC.beh ~ Condition*Sex*order.trial, data=PC2p)

PC2.pairW.par <- emmeans(PC2mod.par, pairwise~Condition*Sex, 
                           adjust="tukey")$contrast %>% 
  data.frame()

#PC3
PC3p <- PCA_comp.par %>% 
  filter(PC == "PC3")
PC3mod.par <- lm(PC.beh ~ Condition*Sex*order.trial, data=PC3p)

PC3.pairW.par <- emmeans(PC3mod.par, pairwise~Condition*Sex, 
                           adjust="tukey")$contrast %>% 
  data.frame()


PCmod.par <- lm(PC.beh ~ Condition*Sex*PC*order.trial, data=PCA_comp.par)

PC_pairW.par <- emmeans(PCmod.par, pairwise~Condition*Sex*PC, 
                           adjust="tukey")


## plot
SF1C <- plot(PC_pairW.par, comparisons = TRUE, col=c('#0979A1', 'black')) + 
  my_theme + coord_flip() + facet_wrap(.~PC, scales = "free") +
  theme(axis.text.x = element_blank())


## NON-Partner
PCA_comp.Npar <- PCA_Npar.All %>% 
  pivot_longer(c("PC1.m", "PC2.m", "PC3.m", "PC1.f", "PC2.f", "PC3.f"),
               names_to = "PC.sex", values_to = "PC.beh") %>% 
  mutate(Sex = rep(c(rep("Male",3), rep("Female",3)),18)) %>% 
  mutate(PC = rep(c("PC1", "PC2", "PC3"),36))

PCA_comp.Npar$Sex <- factor(PCA_comp.Npar$Sex, levels=c("Male", "Female"))

#Linear models per PC followed by pairwise comparisons 

#PC1
PC1np <- PCA_comp.Npar %>% 
  filter(PC == "PC1")
PC1mod.npar <- lm(PC.beh ~ Condition*Sex*order.trial, data=PC1np)

PC1.pairW.npar <- emmeans(PC1mod.npar, pairwise~Condition*Sex, 
                           adjust="tukey")$contrast %>% 
  data.frame()

#PC2
PC2np <- PCA_comp.Npar %>% 
  filter(PC == "PC2")
PC2mod.npar <- lm(PC.beh ~ Condition*Sex*order.trial, data=PC2np)

PC2.pairW.npar <- emmeans(PC2mod.npar, pairwise~Condition*Sex, 
                           adjust="tukey")$contrast %>% 
  data.frame()

#PC3
PC3np <- PCA_comp.Npar %>% 
  filter(PC == "PC3")
PC3mod.npar <- lm(PC.beh ~ Condition*Sex*order.trial, data=PC3np)

PC3.pairW.npar <- emmeans(PC3mod.npar, pairwise~Condition*Sex, 
                           adjust="tukey")$contrast %>% 
  data.frame()



PCmod.Npar <- lm(PC.beh ~ Condition*Sex*PC*order.trial, data=PCA_comp.Npar)

PC_pairW.Npar <- emmeans(PCmod.Npar, pairwise~Condition*Sex*PC, 
                           adjust="tukey")

## plot
SF1E <- plot(PC_pairW.Npar, comparisons = TRUE, col=c('#0979A1', 'black')) + 
  my_theme + coord_flip() + facet_wrap(.~PC, scales = "free") +
  theme(axis.text.x = element_blank())

# Supplementary table 8
PC1.pairW.par %>% 
  bind_rows(PC2.pairW.par) %>% 
  bind_rows(PC3.pairW.par) %>% 
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T) %>% 
  pack_rows("PC1", 1, 6) %>% 
  pack_rows("PC2", 7, 12) %>% 
  pack_rows("PC3", 13, 18)

contrast	estimate	SE	df	t.ratio	p.value
PC1
(Non-Stress Male) - Stress Male	0.62	0.59	22	1.04	0.73
(Non-Stress Male) - (Non-Stress Female)	-0.27	0.47	22	-0.59	0.94
(Non-Stress Male) - Stress Female	0.34	0.59	22	0.57	0.94
Stress Male - (Non-Stress Female)	-0.89	0.59	22	-1.51	0.45
Stress Male - Stress Female	-0.28	0.69	22	-0.40	0.98
(Non-Stress Female) - Stress Female	0.61	0.59	22	1.04	0.73
PC2
(Non-Stress Male) - Stress Male	-0.43	0.68	22	-0.63	0.92
(Non-Stress Male) - (Non-Stress Female)	0.09	0.54	22	0.16	1.00
(Non-Stress Male) - Stress Female	0.05	0.68	22	0.08	1.00
Stress Male - (Non-Stress Female)	0.52	0.68	22	0.76	0.87
Stress Male - Stress Female	0.49	0.80	22	0.61	0.93
(Non-Stress Female) - Stress Female	-0.03	0.68	22	-0.05	1.00
PC3
(Non-Stress Male) - Stress Male	0.39	0.71	22	0.55	0.95
(Non-Stress Male) - (Non-Stress Female)	0.39	0.56	22	0.69	0.90
(Non-Stress Male) - Stress Female	0.40	0.71	22	0.56	0.94
Stress Male - (Non-Stress Female)	0.00	0.71	22	0.00	1.00
Stress Male - Stress Female	0.01	0.83	22	0.01	1.00
(Non-Stress Female) - Stress Female	0.01	0.71	22	0.01	1.00

Table 8: Supplementary Table 8. Summarized results of planned least-squares means contrast between non-stress and stress conditions in males’ and females’ behaviour for partner dyads. P-values were adjusted using Tukey’s method. Results are averaged over the order of the trials.

# Supplementary table 9
PC1.pairW.npar %>% 
  bind_rows(PC2.pairW.npar) %>% 
  bind_rows(PC3.pairW.npar) %>% 
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>% 
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T) %>% 
  pack_rows("PC1", 1, 6) %>% 
  pack_rows("PC2", 7, 12) %>% 
  pack_rows("PC3", 13, 18)

contrast	estimate	SE	df	t.ratio	p.value
PC1
(Non-Stress Male) - Stress Male	-0.06	0.42	28	-0.14	1.00
(Non-Stress Male) - (Non-Stress Female)	-0.90	0.42	28	-2.16	0.16
(Non-Stress Male) - Stress Female	-0.60	0.42	28	-1.43	0.49
Stress Male - (Non-Stress Female)	-0.84	0.42	28	-2.02	0.21
Stress Male - Stress Female	-0.54	0.42	28	-1.29	0.58
(Non-Stress Female) - Stress Female	0.30	0.42	28	0.73	0.89
PC2
(Non-Stress Male) - Stress Male	0.86	0.46	28	1.86	0.27
(Non-Stress Male) - (Non-Stress Female)	0.50	0.46	28	1.08	0.70
(Non-Stress Male) - Stress Female	1.06	0.46	28	2.29	0.12
Stress Male - (Non-Stress Female)	-0.36	0.46	28	-0.77	0.87
Stress Male - Stress Female	0.20	0.46	28	0.44	0.97
(Non-Stress Female) - Stress Female	0.56	0.46	28	1.21	0.62
PC3
(Non-Stress Male) - Stress Male	0.09	0.50	28	0.18	1.00
(Non-Stress Male) - (Non-Stress Female)	-0.22	0.50	28	-0.43	0.97
(Non-Stress Male) - Stress Female	-0.48	0.50	28	-0.97	0.77
Stress Male - (Non-Stress Female)	-0.31	0.50	28	-0.62	0.93
Stress Male - Stress Female	-0.58	0.50	28	-1.15	0.66
(Non-Stress Female) - Stress Female	-0.27	0.50	28	-0.54	0.95

Table 9: Supplementary Table 9. Summarized results of planned least-squares means contrast between non-stress and stress conditions in males’ and females’ behaviour for non-partner dyads. P-values were adjusted using Tukey’s method. Results are averaged over the order of the trials.

##################################################
## Correlations between PCs in every condition ##
##################################################

#We first quickly checked for a mixed model to see whether the behavioural change in males was explained by that of the female partner in two conditions (stress vs. non-stress), and within two familiarities (partners vs. non-partners).
#Although we have 10 pairs, the model summary shows that in fact, female behaviour explains that of the males, but neither familiarity, nor condition or the interaction between them.
#However, the model shows singularity (see ?isSingular), which in this case could indicate that the model is overly complex for the sample size.
#Therefore, we decided to run separate simple correlations for each condition\*familiarity.
#This allow us to compare directly if male partners resemble more the behaviour of their pair bond, after a stressful stimulus, than that of strangers.
#Thus, we proceeded with the simple correlations:


############
#Mixed model
#mm.par <- lmer(PC1.m ~ PC1.f + Condition + (1|Subject), data = PCA_par)
#summary(mm)
############

# Here I mixed databases from partners and non-partners for plotting 
PCA_ALL <- PCA_par.All %>% 
  bind_rows(PCA_Npar.All)

# After separating PCAs based on familiarity, I cannot simply check normality of the scores of each the component, but I have to check based on the discrimination of Familiarity*Condition. I checked Normality for every Familiarity*Condition to see which cor.test to apply: non-normal=Spearman, normal=Pearson.

NORM_test <- PCA_ALL %>% 
  group_by(Familiarity,Condition) %>% 
  normality() # 4 tests indicate non-normality (See bellow).

# Partners-stress
PC1.par.stress <- PCA_par.All %>% 
  filter(Condition=="Stress") %>% 
  cor_test(PC1.f, PC1.m, method = "spearman")

# Partners-Non_stress
PC1.par.nonstress <- PCA_par.All %>% 
  filter(Condition=="Non-Stress") %>% 
  cor_test(PC1.f, PC1.m)

# Non_partners-Stress
PC1.nonpar.stress <- PCA_Npar.All %>% 
  filter(Condition=="Stress") %>% 
  cor_test(PC1.f, PC1.m)

# Non_partners-Non_stress
PC1.nonpar.nonstress <- PCA_Npar.All %>% 
  filter(Condition=="Non-Stress") %>% 
  cor_test(PC1.f, PC1.m, method = "spearman") 
  
# Plot
SF2A <- ggplot(PCA_ALL, aes(PC1.f, PC1.m)) +
  geom_point(shape=21, bg="grey", size=2) +
  my_geom_smooth +
  facet_grid(Condition~Familiarity) + my_theme +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank())


# Partners-stress
PC2.par.stress <- PCA_par.All %>% 
  filter(Condition=="Stress") %>% 
  cor_test(PC2.f, PC2.m)

# Partners-Non_stress
PC2.par.nonstress <- PCA_par.All %>% 
  filter(Condition=="Non-Stress") %>% 
  cor_test(PC2.f, PC2.m) # The p-value is very marginal (0.05) and the rho is moderately high. Thus is weak statistical evidence of behavioural match between f-m.

# Non_partners-Stress
PC2.nonpar.stress <- PCA_Npar.All %>% 
  filter(Condition=="Stress") %>% 
  cor_test(PC2.f, PC2.m, method = "spearman") # This result is unexpected

# Non_partners-Non_stress
PC2.nonpar.nonstress <- PCA_Npar.All %>% 
  filter(Condition=="Non-Stress") %>% 
  cor_test(PC2.f, PC2.m, method = "spearman") # This result is unexpected

# Plot
SF2B <- ggplot(PCA_ALL, aes(PC2.f, PC2.m)) +
  geom_point(shape=21, bg="grey", size=2) +
  my_geom_smooth +
  facet_grid(Condition~Familiarity, scales = "free") + my_theme +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank())

# Partners-stress
PC3.par.stress <- PCA_par.All %>% 
  filter(Condition=="Stress") %>% 
  cor_test(PC3.f, PC3.m)

# Partners-Non_stress
PC3.par.nonstress <- PCA_par.All %>% 
  filter(Condition=="Non-Stress") %>% 
  cor_test(PC3.f, PC3.m)

# Non_partners-Stress
PC3.nonpar.stress <- PCA_Npar.All %>% 
  filter(Condition=="Stress") %>% 
  cor_test(PC3.f, PC3.m) # This result is marginal and unexpected.

# Non_partners-Non_stress
PC3.nonpar.nonstress <- PCA_Npar.All %>% 
  filter(Condition=="Non-Stress") %>% 
  cor_test(PC3.f, PC3.m)


# Plot
SF2C <- ggplot(PCA_ALL, aes(PC3.f, PC3.m)) +
  geom_point(shape=21, bg="grey", size=2) +
  my_geom_smooth +
  facet_grid(Condition~Familiarity, scales = "free") + my_theme +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank())


## We merged the correlation tables among all PC's/condition/familiarity in one table:
PCAsummary <- rbind(
data.frame(rbind(PC1.par.stress,PC2.par.stress[c(1:5,8)],
                 PC3.par.stress[c(1:5,8)]), 
           Condition = rep("Stress",3)),
data.frame(rbind(PC1.par.nonstress[c(1:5,8)],PC2.par.nonstress[c(1:5,8)],
                 PC3.par.nonstress[c(1:5,8)]), 
           Condition = rep("Non-stress",3)),
data.frame(rbind(PC1.nonpar.stress[c(1:5,8)],PC2.nonpar.stress,
                 PC3.nonpar.stress[c(1:5,8)]), 
           Condition = rep("Stress",3)),
data.frame(rbind(PC1.nonpar.nonstress,PC2.nonpar.nonstress,
                 PC3.nonpar.nonstress[c(1:5,8)]), 
           Condition = rep("Non-stress",3))
)


#Supplementary table 10

PCAsummary %>% 
  mutate(Contrasts = paste(var1, var2, sep = " vs. ")) %>% 
  dplyr::select(-var1, -var2) %>% 
  relocate(Contrasts, .before = cor) %>% 
  kable(digits = 3, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>%
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T, align = "c") %>% 
  pack_rows("Partners", 1, 6) %>% 
  pack_rows("Non-partners", 7, 12) 

Contrasts	cor	statistic	p	method	Condition
Partners
PC1.f vs. PC1.m	-0.290	72.000	0.556	Spearman	Stress
PC2.f vs. PC2.m	0.015	0.033	0.975	Pearson	Stress
PC3.f vs. PC3.m	-0.180	-0.415	0.695	Pearson	Stress
PC1.f vs. PC1.m	-0.082	-0.201	0.847	Pearson	Non-stress
PC2.f vs. PC2.m	0.370	0.971	0.369	Pearson	Non-stress
PC3.f vs. PC3.m	0.500	1.433	0.202	Pearson	Non-stress
Non-partners
PC1.f vs. PC1.m	0.140	0.372	0.721	Pearson	Stress
PC2.f vs. PC2.m	0.660	40.669	0.052	Spearman	Stress
PC3.f vs. PC3.m	0.670	2.417	0.046	Pearson	Stress
PC1.f vs. PC1.m	0.068	111.863	0.862	Spearman	Non-stress
PC2.f vs. PC2.m	0.830	20.325	0.006	Spearman	Non-stress
PC3.f vs. PC3.m	0.630	2.146	0.069	Pearson	Non-stress

Table 10: Supplementary Table 10. Correlation estimates of principal components based on behaviours between male and female partners and non-partners, in two conditions (stress/non-stress).

To test pairwise behaviors of dyads, we began by subtracting pre-experimental baseline values from post stress/non-stress values, to eliminate intra-individual sources of variation. Due to resulting zero-inflated data, we then conducted a two-part model following (Boulton and Williford 2018). Briefly, we created two new variables from each original zero-inflated pairwise behavior: 1) a binary variable that indicated whether an observation was zero or non-zero, which was used as a response variable in logistic regression models; and 2) a continuous variable that contained only the non-zero values (zeroes were replaced with NAs), which was used as the response variable in permutation tests. In models, we used the interaction between experimental condition (post-stress/-non-stress), dyad type (partner/non-partner), and trial order as predictors. Because the male-female proximity data contained a large amount of truly missing values, these data were omitted from the former analysis and analyzed separately for pre-treatment baseline and post-treatment condition. In each analysis, we used again a permutation test, with male-female proximity as response variable, and the interaction between experimental condition (stress/non-stress), dyad type, and trial order as predictors.

## New name for behavioral data
bhs.j <- bhs

## To select joint behaviours
bhs.j2 <- na.omit(bhs.j[,c(1:6,8,9,12,24)])

## Subset for the chategorical variables that will be used later for modeling
j.bhs_chrs <- bhs.j2 %>% 
  filter(State == "Pre") %>% 
  dplyr::select(Familiarity:Subject, ord.trial)

## Subset for baseline
j.bhs_pre <- bhs.j2 %>% 
  filter(State == "Pre") #%>% 
  #select(Subject:JointRefuge) %>%
  #arrange(desc(Subject)) %>% 
  #mutate_if(is.character,as.numeric) %>% 
  #select(CordFreezingB:JointRefuge)

## Subset for treatment
j.bhs_post <- bhs.j2 %>% 
  filter(State == "Post") #%>% 
  #select(Subject:JointRefuge) %>%
  #arrange(desc(Subject)) %>% 
  #mutate_if(is.character,as.numeric) %>% 
  #select(CordFreezingB:JointRefuge)

## Subtracted data
j.bhs <- j.bhs_post[,7:9]-j.bhs_pre[,7:9]

## Joint behaviours with categorical variables
J.bhs2 <- data.frame(j.bhs_chrs, j.bhs)

## Since joint baheviours are the same for males and females, we used one sex
J.bhs3 <- J.bhs2 %>% 
  filter(Sex=="Male")

# New binary variables

Bin_bhe <- data.frame(J.bhs3$Familiarity, J.bhs3$Condition,
                      J.bhs3$ord.trial,
                      ifelse(J.bhs3$CordFreezingB != 0,1,0),
                      ifelse(J.bhs3$CordFreezingD != 0,1,0),
                      ifelse(J.bhs3$JointRefuge != 0,1,0))

colnames(Bin_bhe) <- c("Familiarity", "Condition", "order.trial", "CordFreezingB", 
                       "CordFreezingD", "JointRefuge")

# New continuous variables variables

J.bhs3$CordFreezingB[J.bhs3$CordFreezingB == 0] <- NA # CoordFreezingB
J.bhs3$CordFreezingD[J.bhs3$CordFreezingD == 0] <- NA # CoordFreezingD
J.bhs3$JointRefuge[J.bhs3$JointRefuge == 0] <- NA # JointRefuge

# New binary outcomes
Bin_bhe$CordFreezingB[is.na(Bin_bhe$CordFreezingB)] <- 0 # CoordFreezingB
Bin_bhe$CordFreezingD[is.na(Bin_bhe$CordFreezingD)] <- 0 # CoordFreezingD
Bin_bhe$JointRefuge[is.na(Bin_bhe$JointRefuge)] <- 0 # JointRefuge

CordFB.Mod <- glm(CordFreezingB ~ paste(Familiarity, Condition) + order.trial, 
                  data=Bin_bhe, family="binomial")

CordFD.Mod <- glm(CordFreezingD ~ paste(Familiarity, Condition) + order.trial, 
                  data=Bin_bhe, family="binomial")

JointR.Mod <- glm(JointRefuge ~ paste(Familiarity, Condition) + order.trial, 
                  data=Bin_bhe, family="binomial")

log.part <- as.data.frame(rbind(summary(CordFB.Mod)$coefficients,
      summary(CordFD.Mod)$coefficients,
      summary(JointR.Mod)$coefficients))


## Supplementary table 11

log.part %>% 
  rownames_to_column("Predictor") %>% 
  mutate(Predictor = rep(c("Intercept", "Non-partner (Stress)", "Partner (Non-stress)",
                           "Partner (Stress)", "Order trial"), 3)) %>% 
  mutate("P" = .$`Pr(>|z|)`) %>% 
  dplyr::select(Predictor, Estimate, "Std. Error", "z value", P) %>% 
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>%
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T, align = "c") %>% 
  pack_rows("Coordinated freezing (bouts)", 1, 5) %>% 
  pack_rows("Coordinated freezing (duration)", 6, 10) %>% 
  pack_rows("Join refuge", 11, 15)

Predictor	Estimate	Std. Error	z value	P
Coordinated freezing (bouts)
Intercept	-0.48	1.34	-0.36	0.72
Non-partner (Stress)	0.51	1.08	0.47	0.64
Partner (Non-stress)	0.99	1.05	0.94	0.35
Partner (Stress)	0.99	1.05	0.94	0.35
Order trial	-0.51	0.72	-0.71	0.48
Coordinated freezing (duration)
Intercept	-0.92	1.34	-0.68	0.49
Non-partner (Stress)	1.01	1.05	0.96	0.34
Partner (Non-stress)	1.01	1.05	0.96	0.34
Partner (Stress)	1.01	1.05	0.96	0.34
Order trial	-0.22	0.70	-0.31	0.76
Join refuge
Intercept	-1.25	1.95	-0.64	0.52
Non-partner (Stress)	-19.31	5910.12	0.00	1.00
Partner (Non-stress)	0.00	1.14	0.00	1.00
Partner (Stress)	-19.31	5910.12	0.00	1.00
Order trial	0.00	1.14	0.00	1.00

Table 11: Supplementary Table 11. Result summary statistics of logistic regressions for coordinated pairwise stress response behaviors adjusted for the order of the trials.

In brief, males did not behaviorally state match females in any experimental condition, except for gaze and approach towards non-partner females in the non-stress condition (R2 = 0.83, S = 20.32; P = 0.006), and freezing towards non-partner females in the stress condition (R2 = 0.67, S = 2.41; P = 0.04; Supplementary Table 12 and Supplementary Figure 2 A-C for details).

Figure 5: Supplementary Figure 2. Overall, males do not behaviorally state match females, nor do they increase dyadic fear responses with females after females are subjected to a stressor treatment. (A-C) Overall, males do not behaviorally state match female partners or non-partners in any condition. Shown are male-female correlations of 9 behaviors reduced to 3 non-redundant principal component categories: (A) male-female interaction, (B) activity, and (C) stillness levels, after subtracting pre-experimental baseline levels. (D-G) Likewise, the stress treatment does not affect coordinated dyadic fear response behaviors relative to experimental pre-treatment baseline or post- non-stress treatment conditions, in either partner or non-partner dyads. Values shown are relative to the experimental pre-treatment (baseline) condition, except for m-f proximity, which is an absolute value. Bars represent the mean and whiskers denote 95% confidence intervals. NS = non-significant (P ≥ 0.05).

## Function to summarize data
data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      CI = qnorm(0.975)*(sd(x[[col]], na.rm=TRUE))/sqrt(length(x[[col]])))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}


# Data summary
CooB_sum <- data_summary(J.bhs3, varname="CordFreezingB", 
                        groupnames=c("Familiarity", "Condition"))


## independence tests for continuous variables

# To merge familiarity and condition
bhs.j$Ncat <- paste(bhs.j$Familiarity, bhs.j$Condition, sep = "_")

J.bhs3$Ncat <- paste(J.bhs3$Familiarity, J.bhs3$Condition, sep = "_")

## Independence test for Freezing frequency
anco1 <- independence_test(CordFreezingB ~ as.factor(Ncat) * ord.trial, data = J.bhs3, 
                            ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
                            teststat = "quad", distribution = approximate(nresample = 10000))

anco1.out <- data.frame("p-value" = pvalue(anco1)[1], "chi-squared" = statistic(anco1)[1], 
           "Variable" = "Coord. freezing (bouts)")

# The plot
SF2D <- ggplot(CooB_sum, aes(x=Familiarity, y=CordFreezingB, fill=Condition)) +
  geom_col(stat="identity", width=0.7, position = position_dodge(0.7), 
           col="black", lwd=.3) + 
  geom_errorbar(aes(ymin=CordFreezingB-CI, 
                    ymax=CordFreezingB+CI), 
                position = position_dodge(0.7), width=.09, lwd=0.3) +
  my_theme +
  scale_fill_manual(values = c("black", "white")) +
  theme(axis.title.x = element_blank())


# Data summary
CooD_sum <- data_summary(J.bhs3, varname="CordFreezingD", 
                        groupnames=c("Familiarity", "Condition"))

## Independence test for Freezing duration
anco2 <- independence_test(CordFreezingD ~ as.factor(Ncat) * ord.trial, data = J.bhs3, 
                            ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
                            teststat = "quad", distribution = approximate(nresample = 10000))

anco2.out <- data.frame("p-value" = pvalue(anco2)[1], "chi-squared" = statistic(anco2)[1], 
           "Variable" = "Coord. freezing (duration)")

# The plot
CooD <- ggplot(CooD_sum, aes(x=Familiarity, y=CordFreezingD, fill=Condition)) +
  geom_col(stat="identity", width=0.7, position = position_dodge(0.7), 
           col="black", lwd=.3) + 
  geom_errorbar(aes(ymin=CordFreezingD-CI, 
                    ymax=CordFreezingD+CI), 
                position = position_dodge(0.7), width=.09, lwd=0.3) +
  my_theme +
  scale_fill_manual(values = c("black", "white")) +
  theme(axis.title.x = element_blank())

# Data summary
JRef_sum <- data_summary(J.bhs3, varname="JointRefuge", 
                        groupnames=c("Familiarity", "Condition"))

## Independence test for joint refuge duration
anco3 <- independence_test(JointRefuge ~ as.factor(Ncat) * ord.trial, data = J.bhs3, 
                            ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
                            teststat = "quad", distribution = approximate(nresample = 10000))

anco3.out <- data.frame("p-value" = pvalue(anco3)[1], "chi-squared" = statistic(anco3)[1], 
           "Variable" = "Join refuge")

# The plot
SF2F <- ggplot(JRef_sum, aes(x=Familiarity, y=JointRefuge, fill=Condition)) +
  geom_bar(stat="identity", width=0.7, position = position_dodge(0.7), 
           col="black", lwd=.3) + 
  geom_errorbar(aes(ymin=JointRefuge-CI, 
                    ymax=JointRefuge+CI), 
                position = position_dodge(0.7), width=.09, lwd=0.3) + 
  my_theme +
  scale_fill_manual(values = c("black", "white")) +
  theme(legend.position = "none") +
  theme(axis.title.x = element_blank())

#Because the variable male-female proximity is unbalanced because of the missing values, subtracting post - pre treatment would render very few observations to analyse.
#Therefore, we separated the analysis for baseline and treatment:
  
## Baseline (pre)
pre.m_f <- bhs.j %>% 
  filter(State=="Pre", Sex=="Male") %>% 
  mutate(NCat = J.bhs3$Ncat) %>% 
  dplyr::select(State, Familiarity, Condition, AV_M.F_proximity_cm, Ncat, ord.trial)


## Data summary
MF.pre_sum <- data_summary(pre.m_f, varname="AV_M.F_proximity_cm", 
                        groupnames="Familiarity") %>% 
  mutate(Condition = rep("Baseline", 2))

## Independence test m-f proximity pre-stimulus
anco4 <- independence_test(AV_M.F_proximity_cm ~ as.factor(Ncat) * ord.trial, data = pre.m_f, 
                            ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
                            teststat = "quad", distribution = approximate(nresample = 10000))

anco4.out <- data.frame("p-value" = pvalue(anco4)[1], "chi-squared" = statistic(anco4)[1], 
           "Variable" = "m-f proximity (partners)")


## Treatment (post)

post.m_f <- bhs.j %>% 
  filter(State=="Post", Sex=="Male") %>% 
  mutate(NCat = J.bhs3$Ncat) %>% 
  dplyr::select(State, Familiarity, Condition, AV_M.F_proximity_cm, Ncat, ord.trial)

# Data summary
MF.post_sum <- data_summary(post.m_f, varname="AV_M.F_proximity_cm", 
                        groupnames=c("Familiarity", "Condition"))

## Independence test m-f proximity post-stimulus
anco5 <- independence_test(AV_M.F_proximity_cm ~ as.factor(Ncat) * ord.trial, data = post.m_f, 
                            ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
                            teststat = "quad", distribution = approximate(nresample = 10000))

anco5.out <- data.frame("p-value" = pvalue(anco5)[1], "chi-squared" = statistic(anco5)[1], 
           "Variable" = "m-f proximity (non-partners)")


# Merge pre and post summaries
MF.all_sum <- rbind (MF.pre_sum,MF.post_sum)

MF.all_sum$Familiarity <- factor(MF.all_sum$Familiarity, levels=c("Partner",
                                                                  "Non-Partner"))

MF.all_sum$Condition <- factor(MF.all_sum$Condition, levels=c("Baseline",
                                                              "Stress",
                                                              "Non-Stress"))
  
# The plot
SF2G <- ggplot(MF.all_sum, aes(x=Familiarity, y=AV_M.F_proximity_cm, fill=Condition)) +
  geom_bar(stat="identity", width=0.7, position = position_dodge(0.7), 
           col="black", lwd=.3) + 
  geom_errorbar(aes(ymin=AV_M.F_proximity_cm-CI, 
                    ymax=AV_M.F_proximity_cm+CI), 
                position = position_dodge(0.7), width=.09, lwd=0.3) + 
  my_theme +
  scale_fill_manual(values = c("grey","white", "black"))

#Supplementary table 12
anco1.out %>%
  bind_rows(anco2.out) %>% 
  bind_rows(anco3.out) %>%
  bind_rows(anco4.out) %>%
  bind_rows(anco5.out) %>%
  kable(digits = 2, 
      table.attr = 'data-quarto-disable-processing="true"', "html") %>%
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  row_spec(0, italic = T, bold = T, align = "c") 

p.value	chi.squared	Variable
0.58	3.24	Coord. freezing (bouts)
0.90	1.23	Coord. freezing (duration)
1.00	0.00	Join refuge
0.16	6.47	m-f proximity (partners)
0.29	5.12	m-f proximity (non-partners)

Table 12: Supplementary Table 12. Result summary statistics of permutation test for coordinated pairwise stress response behaviors adjusted for the order of the trials. P-values were adjusted for multiple comparisons using the Bonferroni method.

References

Boulton, Aaron J., and Anne Williford. 2018. “Analyzing Skewed Continuous Outcomes With Many Zeros: A Tutorial for Social Work and Youth Prevention Science Researchers.” Journal of the Society for Social Work and Research 9 (4): 721–40. https://doi.org/10.1086/701235.

Burkett, J. P., E. Andari, Z. V. Johnson, D. C. Curry, F. B. M. De Waal, and L. J. Young. 2016. “Oxytocin-Dependent Consolation Behavior in Rodents.” Science 351 (6271): 375–78. https://doi.org/10.1126/science.aac4785.

Hartig, Florian. 2022. “DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models.” https://CRAN.R-project.org/package=DHARMa.

Lenth, Russell V. 2023. “Emmeans: Estimated Marginal Means, Aka Least-Squares Means.” https://CRAN.R-project.org/package=emmeans.

William Revelle. 2023. “Psych: Procedures for Psychological, Psychometric, and Personality Research.” https://CRAN.R-project.org/package=psych.