Statistical modelling
We used generalized linear mixed models (GLMMs) to analyze the effect of
the temperature treatment on performance linked traits of the host and
compare differences in performance between ranges. All traits were
measured at the end of the experiment (t12) and
therefore had, in contrast to the microbial count data, no repeated
measures. The models were fitted as a function of temperature,
range, their interaction and the nested structure of population
identity inside range . The random intercept individual
identity was included in each model to represent the parent individual
from which the thallus in the 15 °C group and the thallus in the 22 °C
group originated. However, when individual identity did not explain
variation, it was excluded and the model was run without random effects.
We used a Gaussian distribution (with identity in the link
function) to fit GLMMs on the RGR and the photosynthetic yield data. To
satisfy the model assumptions photosynthetic yield was rank transformed.
The disease symptoms (thallus brittleness and tissue decay) were
modelled using a binomial distribution with a logit in the link
function.
Diversity was analyzed as rarefied OTU richness (all samples were
subsampled to 1171 counts, the number of reads of the smallest sample
size) and evenness (measured as the probability of interspecific
encounter, PIE), calculated with the R package mobr (McGlinn et al.,
2019). GLMMs were then fitted on two subsets of the data. The first
subset comprised the samples from the field (tf), before
the disturbance treatment (t0) and the samples from the
first timepoint (t1) after the disturbance treatment at
15 °C. This model included time as factorial variable,range , the corresponding interaction, and individual
identity nested in population identity as random intercepts. The
second subset comprised all post-disturbance timepoints
(t1, t2, t4 and
t12) and was modelled as a function of time ,temperature , range , the interactions, and individual
identity nested in population identity as random intercepts.
We explored variation in community composition with non-metric
dimensional scaling (nMDS) on the full dataset (including all time
points) and the post-disturbance dataset (t1,
t2, t4 and t12), based
on Bray-Curtis and Euclidean distances. Using the R package mvabund (Y.
I. Wang, Naumann, Wright, & Warton, 2012), multivariate generalized
linear models (mGLMs) were fitted on the count matrix with the log
transformed sequencing depth (LSD) as an offset prior to the scaling
procedure to adjust for the effect of different read depths per sample.
We then fitted mGLMs on the post-disturbance data, including different
combinations of the variables time, temperature, and range to visualize
their partial effects. All mGLMs assumed a negative binomial
distribution with a natural logarithm in the link function. The nMDS was
then conducted on the residuals of the models that were back transformed
to the original scale. Group centroids were computed and their 95%
confidence regions were computed using the R package vegan (Oksanen et
al., 2017). Statistics for the community response were obtained from the
mGLM fitted on the post-disturbance dataset that included all variables
(i.e., time , temperature , range , the corresponding
interactions and population identity nested in range ),
using the anova.manyglm function from the mvabund package, bootstrapping
the univariate models with 500 iterations and individual identity as a
blocking factor.
To analyze beta diversity within-populations we used Bray-Curtis and
Euclidean distances calculated from the same community matrix as was
used for the nMDS (i.e., adjusted for the effect of the sequencing
depth). After computing the distances, we analyzed two parts of our
dataset. We only considered distances between samples from the same
range and the same temperature regimes. The first subset was prepared to
analyze the effect of the disturbance on the beta diversity and included
the pre-disturbance samples (tfield, t0)
and the first post-disturbance time point (t1) in the 15
°C treatment. The GLMM used for the first model was a function oftime as factor, range, and their interaction. The
multilevel factors population and individual combinationwere included as random intercepts to account for dependency of
distances calculated from the same individual pair. The post-disturbance
subset of the data was used to analyze the effect of range and
temperature on the microbiota dispersion. The GLMMs included the
variables time , temperature , range , and all
possible interactions.