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.