Chapter 8 Do populations mantain a different GM after acclimation?
load("data/data_03122025.Rdata")
load("data/objects_03122025.Rdata")
load("data/alpha_16092025.Rdata")
load("data/beta_16092025.Rdata")8.1 Shapiro test
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="Acclimation") %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.2524848alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="Acclimation") %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.4351805alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point =="Acclimation") %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.32429058.2 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(method="t.test",size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_blank(),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
8.3 Beta diversity
Number of samples used
samples_to_keep_accli <- sample_metadata %>%
filter(time_point == "Acclimation") %>%
select(Tube_code) %>%
pull()
subset_meta_accli <- sample_metadata %>%
filter(time_point == "Acclimation")
length(samples_to_keep_accli)[1] 26Richness
richness_accli <- as.matrix(beta_q0n$S)
richness_accli <- as.dist(richness_accli[rownames(richness_accli) %in% samples_to_keep_accli,
colnames(richness_accli) %in% samples_to_keep_accli])
betadisper(richness_accli, subset_meta_accli$species) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.094207 0.094207 11.879 999 0.004 **
Residuals 24 0.190338 0.007931
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.003
Podarcis_muralis 0.0021024 adonis2(richness_accli ~ species,
data = subset_meta_accli %>% arrange(match(Tube_code,labels(richness_accli))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.651732 | 0.1945425 | 5.796731 | 0.001 |
| Residual | 24 | 6.838608 | 0.8054575 | NA | NA |
| Total | 25 | 8.490340 | 1.0000000 | NA | NA |
Neutral
neutral_accli <- as.matrix(beta_q1n$S)
neutral_accli <- as.dist(neutral_accli[rownames(neutral_accli) %in% samples_to_keep_accli,
colnames(neutral_accli) %in% samples_to_keep_accli])
betadisper(neutral_accli, subset_meta_accli$species) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.05661 0.056610 4.2136 999 0.048 *
Residuals 24 0.32244 0.013435
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.04
Podarcis_muralis 0.051158 adonis2(neutral_accli ~ species,
data = subset_meta_accli %>% arrange(match(Tube_code,labels(neutral_accli))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.965807 | 0.2500091 | 8.00039 | 0.001 |
| Residual | 24 | 5.897134 | 0.7499909 | NA | NA |
| Total | 25 | 7.862941 | 1.0000000 | NA | NA |
Phylogenetic
phylo_accli <- as.matrix(beta_q1p$S)
phylo_accli <- as.dist(phylo_accli[rownames(phylo_accli) %in% samples_to_keep_accli,
colnames(phylo_accli) %in% samples_to_keep_accli])
betadisper(phylo_accli, subset_meta_accli$species) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03621 0.036209 2.2972 999 0.148
Residuals 24 0.37829 0.015762
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.138
Podarcis_muralis 0.14267 adonis2(phylo_accli ~ species,
data = subset_meta_accli %>% arrange(match(Tube_code,labels(phylo_accli))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.2126877 | 0.1387856 | 3.867626 | 0.016 |
| Residual | 24 | 1.3198034 | 0.8612144 | NA | NA |
| Total | 25 | 1.5324911 | 1.0000000 | NA | NA |
NMDS
richness_accli %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_accli, by = join_by(sample == Tube_code))%>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species)) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080',"#d57d2c")) +
scale_shape_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("circle","square","triangle")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
neutral_accli %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_accli, by = join_by(sample == Tube_code))%>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species, shape=type)) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080',"#d57d2c")) +
scale_shape_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("circle","square","triangle")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
phylo_accli %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta_accli, by = join_by(sample == Tube_code))%>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species)) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080',"#d57d2c")) +
scale_shape_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold_intervention"),
labels=c("Cold_control", "Warm_control", "Cold_intervention"),
values=c("circle","square","triangle")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (),x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none") 
8.4 Functional capacity
8.4.1 Metabolic capacity index at functional level
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="Acclimation") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.348 0.0158
2 Podarcis_muralis 0.329 0.0320MCI_element_acclimation<- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="Acclimation")
shapiro.test(MCI_element_acclimation$value) #normality test
Shapiro-Wilk normality test
data: MCI_element_acclimation$value
W = 0.92669, p-value = 0.06465
F test to compare two variances
data: value by species
F = 0.24311, num df = 8, denom df = 16, p-value = 0.04866
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.0777983 0.9909228
sample estimates:
ratio of variances
0.2431059
Two Sample t-test
data: value by species
t = 1.6212, df = 24, p-value = 0.1181
alternative hypothesis: true difference in means between group Podarcis_liolepis and group Podarcis_muralis is not equal to 0
95 percent confidence interval:
-0.005046497 0.042002844
sample estimates:
mean in group Podarcis_liolepis mean in group Podarcis_muralis
0.3478327 0.3293545 8.4.2 Differences in functional pathways
element_gift_acclimation<- element_gift %>%
filter(time_point=="Acclimation") %>%
select(-time_point,-type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where(~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 4)], by = "Tube_code")significant_elements_acclimation <- element_gift_acclimation %>%
pivot_longer(-c(Tube_code,Population), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ Population)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
select(-Elements)
element_gift_sig_acclimation <- element_gift_acclimation %>%
select(Tube_code, all_of(intersect(
significant_elements_acclimation$trait,
colnames(element_gift_acclimation)
))) %>%
left_join(., sample_metadata[c(1, 4)], by = join_by(Tube_code == Tube_code))
difference_table_acclimation <- element_gift_sig_acclimation %>%
select(-Tube_code) %>%
group_by(Population) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Cold-Warm)%>%
mutate(group_color = ifelse(Difference <0, "Warm","Cold")) difference_table_acclimation %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c('#008080',"#d57d2c")) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
legend.title = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Function") +
ylab("Mean difference")