Negative and positive interspecific interactions involving jellyfish polyps in marine sessile communities

Sessile marine invertebrates on hard substrates are one of the two canonical examples of communities structured by competition, but some aspects of their dynamics remain poorly understood. Jellyfish polyps are an important but under-studied component of these communities. We determined how jellyfish polyps interact with their potential competitors in sessile marine hard-substrate communities, using a combination of experiments and modelling. We carried out an experimental study of the interaction between polyps of the moon jellyfish Aurelia aurita and potential competitors on settlement panels, in which we determined the effects of reduction in relative abundance of either A. aurita or potential competitors at two depths. We predicted that removal of potential competitors would result in a relative increase in A. aurita that would not depend on depth, and that removal of A. aurita would result in a relative increase in potential competitors that would be stronger at shallower depths, where oxygen is less likely to be limiting. Removal of potential competitors resulted in a relative increase in A. aurita at both depths, as predicted. Unexpectedly, removal of A. aurita resulted in a relative decrease in potential competitors at both depths. We investigated a range of models of competition for space, of which the most successful involved enhanced overgrowth of A. aurita by potential competitors, but none of these models was completely able to reproduce the observed pattern. Our results suggest that interspecific interactions in this canonical example of a competitive system are more complex than is generally believed.

(S1) subcomposition onto a subspace that preserves variation between A. aurita, bare panel and the geometric mean 21 of potential competitors, but ignores variation among potential competitors (Pawlowsky-Glahn et al., 2015, 22 p. 51). Ternary diagrams in such subspaces are the preferred way to represent marginal relationships among 23 compositions with more than three parts, because this approach is consistent with the Aitchison geometry 24 (van den Boogaart and Tolosana-Delgado, 2013, section 4.2.1). In addition, the choice of v 2 as a logcontrast 25 between geometric means results in a lower-dimensional representation in which differences have a population-26 dynamic interpretation. For two subcompositions c (1) , c (2) , a difference in the coefficient of v 1 is proportional 27 to log(c 2 ), which has the form of a difference in proportional growth rates between A.
(S2) 39 We coded treatment effects in the compositional manova using orthogonal contrasts. We represented treat-40 ment and depth effects by the matrices F, G respectively: normal prior centred on the posterior mean vector from Chong and Spencer (2018). However, because the experiment was done in a different year and on a different substrate, and the taxa in Chong and Spencer (2018) 68 were amalgamated, we made our prior more diffuse and with a weaker correlation structure ( Figure S1, orange): 69 α ∼ N (µ α , diag(2s)R diag(2s)), 70 where µ α = (2.8, −2.6, −0.1, 0.1) T is the sample posterior mean vector, s is the vector of sample posterior 71 standard deviations, and R is a correlation matrix whose off-diagonal elements are all set to half the sample 72 mean posterior correlation among all pairs of components, giving the prior covariance matrix 73 Σ α = diag(2s)R diag(2s)) =          12.7 3.2 × 10 −4 5.5 × 10 −5 8.5 × 10 −5 3.2 × 10 −4 7.1 4.1 × 10 −5 6.4 × 10 −5 5.5 × 10 −5 4.1 × 10 −5 0.21 1.1 × 10 −5 8.5 × 10 −5 6.4 × 10 −5 1.1 × 10 −5 0.5 This is close to a diagonal covariance matrix. Although it is a reasonable description of most aspects of the 75 posterior distribution from Chong and Spencer (2018), it does not capture the strong correlation between α 1 76 and α 2 . This is probably an advantage: we are not confident that data from an observational study in a different 77 year are sufficient to support very strong prior information.

78
For the covariance matrices Σ and Z of the panel effects ε and block effects δ respectively, we used generic of the corresponding effect for dock wall images from Chong and Spencer (2018) because our panels were 83 smaller than the dock wall images, and in past experiments (Maxatova, 2016;Edney, 2017;Sharpe, 2020), there 84 appeared to be more variation among panels than adjacent areas of dock wall. We used the same prior for block 85 effects because previous experiments (Maxatova, 2016;Edney, 2017;Presser, 2019;Sharpe, 2020) suggested that 86 among-block differences could plausibly be as large as among-panel differences.

87
For the overall mean µ in ilr coordinates, we chose independent univariate N (0, 4) priors on each coordinate 88 (where the second parameter is the variance We compared versions of the compositional manova using leave-one-cluster-out cross-validation. The loss func-135 tion was the Bayesian leave-one-cluster-out estimate of out-of-sample prediction error (expected log predictive 136 density) elpd loco , where a cluster corresponds to a block of panels: Here, L is the number of blocks and f (y l |y −l ) is the posterior density of the lth block y l , given the data set 139 y −l in which block l is excluded (Vehtari et al., 2017). The leave-one-cluster-out posterior density is where θ is the parameter vector {µ, α, β 1 , β 2 , γ 1 , γ 2 , Σ, Z}, f (y l |θ) is the density of the lth block given param- f (δ|θ)f (ε|θ)f (y jkl |θ, δ, ε) dε dδ,

146
where f (δ|θ) and f (ε|θ) are the multivariate normal densities of block and panel effects, with mean vector 0 147 and covariance matrices Z and Σ respectively, and f (y jkl |θ, δ, ε) is the multinomial probability for the counts 148 from the panel at depth j, removal treatment k in block l (Equation S2). We estimated the integral in Equation

149
S3 using samples of size 1000 from the posterior density of f (θ|y −l ), obtained using In the compositional manova, let φ jk = α j ⊕ β k ⊕ γ jk ∈ S 4 denote a treatment effect (the effect of the projection of a treatment effect onto a subset of these basis vectors with indices S is given by where i∈S denotes repeated perturbation, ·, · a denotes the Aitchison inner product (Egozcue et al., 2003) and 162 represents the powering operator (Aitchison, 1986, p. 120). With S = {1, 2}, we obtain the projection of treat-163 ment effects onto the 2-simplex with parts representing A. aurita, bare panel and gm (potential competitors), 164 where gm () denotes the geometric mean. With S = {3, 4}, we obtain treatment effects on the subcomposition 165 of potential competitors, with parts Botrylloides spp., Bugula spp. and Molgula tubifera. 166 We represented each of these projections of treatment effects in a ternary plot. We calculated 95 % highest 167 posterior density credible regions for the projections of treatment effects using the algorithm in Hyndman 168 (1996), implemented in the hdr.2d() function in R package hdrcde 3.4 (Hyndman, 2018). We plotted the 169 corresponding observations in the formρ ijkl μ , whereρ ijkl denotes the sample estimate of the relative 170 abundance vector in panel i from depth j, treatment k, block l (with zeros replaced by 1/2, which is half the 171 detection limit),μ denotes the estimated posterior mean overall metric centre from the manova, and is the 172 compositional difference operator (Egozcue et al., 2003). 173 We visualized the panel and block covariance matrices (Σ and Z respectively) by constructing ellipses of 174 unit Mahalanobis distance around the origin in the first two and last two ilr coordinates. We constructed 175 100 such ellipses using draws from the posterior distributions of the corresponding submatrices of Σ and Z, 176 back-transformed and plotted them in ternary plots.

177
S7 Measuring treatment effects 178 We assessed the effects of potential competitors on A. aurita using differences in logit (A. aurita) between 179 potential competitor removal (O) and control (C) treatments. Similarly, we assessed the effects of A. aurita 180 on potential competitors using differences in logit (potential competitors) between A. aurita removal (A) and 181 control (C) treatments.

182
Differences in treatment metric centres from the manova are not directly relevant to our hypotheses, because 183 we expect the metric centre to be affected by a removal treatment, even if the non-target taxon does not respond 184 to this removal. In contrast, logit differences take the value zero if the non-target taxon does not respond. For The removal treatment has no direct effect on logit(c 1 ) = log(c 1 /(1 − c 1 )), because we are simply turning some In the main text, we do not include polyp growth on potential competitors, because this was very rarely 212 observed in the experimental data (although it is common in some years). Here, we describe and analyze the 213 basic properties of a model with a third state variable representing polyps on potential competitors (the model 214 we specified before collecting data), because it may be useful for future work. We show that in this model, it is

217
To the model in the main text we add a third state variable, the density y 2 of A. aurita polyps on potential 218 competitors, per unit surface area of substrate (numbers L −2 ). We separate polyp density into those on substrate 219 and those on potential competitors because polyps are able to settle on some organisms such as the ascidian 220 Ascidiella aspersa and the bivalve Mytilus edulis. In this model, increases in the proportion of substrate filled by potential competitors could potentially have either a negative or a positive net effect on polyps, depending 222 on the amount of new surface area created and the ability of polyps to settle on this new area.

223
Our updated model is The processes included in this model are sketched in Figure S3.

248
We determine the signs of the elements of the community matrix for this model below (section S10). and 4, we obtain the dimensionless system x, y * 1 . Then the community matrix C is (1, 2) and (2, 1)). 281 S10 Dimensionless form and community matrix for model including 282 polyp growth on potential competitors 283 We take the same approach to nondimensionalizing the model with polyp growth on potential competitors 284 (section S8) as for the model in the main text. Let y * 2 = δy 2 be the area filled by polyps on potential competitors 285 per unit substrate area: this may be greater than 1, because ψ ≥ 1. Substituting this and the other dimensionless 286 variables from the main text into Equations S6, S7 and S8, we obtain the dimensionless system

297
We use the community matrix to measure interaction strengths, as for the model in the main text, addition-298 ally defining h(x, y * 1 , y * 2 ) = 1 y * 2 dy * 2 dt * to be the proportional rate of change of the dimensionless variable y * 2 . Then 299 the community matrix is As in the main text, 0 ≤ x ≤ 1, 0 ≤ y * 1 ≤ 1, 0 ≤ y * 2 , 0 ≤ x + y * 1 ≤ 1, 0 ≤ ψx − y * 2 ≤ 1, ψ ≥ 1, and Π 3 is the only 302 negative parameter, so the signs of the elements in the matrix are In addition to the effects discussed in the main text, polyps on substrate and polyps on potential competitors 305 have positive effects on each other (elements (2, 3) and (3, 2)), and polyps on potential competitors have no 306 direct effect on potential competitors (element (1, 3)). The sign of the effect of potential competitors on polyps 307 growing on potential competitors (element (3, 1)) is unknown until parameter values are known. However, the 308 sign does not depend on ψ, the surface area available for occupation by polyps per unit substrate area occupied 309 by potential competitors, although the magnitude does. Furthermore, the sign is more likely to be negative if 310 y * 2 , the area filled by polyps on potential competitors per unit substrate area, is large, the ratio Π 7 of settlement 311 rates of polyps on potential competitors to substrate is small, the ratio Π 10 of polyp budding rate onto potential  For reference, the model with facilitation of settlement is The dimensionless form of Equation S13 is

340
For reference, the model for facilitation of growth is The dimensionless form of Equation S14 is

S11.3 Overgrowth of polyps by potential competitors
For reference, the model for overgrowth of polyps by potential competitors is The dimensionless form of Equations S15 and S16 is

S11.4 Protection from predators 358
For reference, the model for protection from predators is The dimensionless form of Equation S17 is where m 2 is already dimensionless. The community matrix is 365 S12 Fitting dynamic models to experimental data 366 We fitted versions of Equations 3 and 4, with each of the modifications in section 2.3.2 in turn, to the experi-367 mental data from all weeks and panels. As explained below, our underlying model for dynamics is deterministic, 368 so we will only obtain a caricature of the true dynamics. Nevertheless, this will give a qualitative understanding 369 of the interactions in the community. We wrote Equations 3 and 4 in terms of proportions of space occupied, x and y * 1 , but retained time in dimensioned form (measured in weeks since panels were put in the water). Thus 371 the basic model is: In the version with overgrowth of polyps by potential competitors, the overgrowth effect is then +a 1,y * 1 xy * 1 in 376 dx dt , and −a 1,y * 1 xy * 1 in dy * 1 dt , where a 1,y * 1 = a 1,y1 /δ.

377
Let y Bugula spp and M. tubifera, as in the models for final composition, but amalgamated into a single part) on the 379 panel from depth j, treatment k, block l at time t, pre-treatment in those cases where a treatment was applied.

380
We model these counts using: where x jkl (t) and y * 1,jkl (t) are the solutions to Equations S18 and S19 (or modified versions as in section 2.3.2) 385 for the panel with depth j, treatment k in block l, obtained using a fourth and fifth order Runge-Kutta method, 386 with initial conditions x jkl (0) = 0, y * 1,jkl (0) = 0 (the empty panel). We estimated all the parameters (other 387 than removal effects, as described below) separately at depths 1 m and 3 m. 388 We did not include block and panel effects. These should really be the outcome of block-and panel-specific 389 temporal variation in parameters such as settlement and growth rates in a stochastic differential equation version 390 of Equations S18 and S19. We did not attempt to fit such a model because it presents substantial technical 391 difficulties. First, each parameter has a fixed sign (for example, settlement rates a 0 and δb 0 must always be 392 positive). Care must be taken to respect such conditions when specifying a stochastic differential equation.

393
The usual stochastic differential equation models for community dynamics avoid the problem by allowing only 394 proportional population growth rates at low density (which do not have a sign constraint) to be stochastic 395 (Kloeden and Platen, 1999, p. 254). Second, we cannot obtain an explicit solution for a stochastic differential 396 equation version of Equations S18 and S19, and would therefore have to rely on numerical methods such as fall outside the simplex. These problems are not insoluble, but are outside the scope of this paper.
Where a treatment was applied, we model the post-treatment counts y   Figure S4a).

431
Settlement and growth rates of potential competitors may be somewhat higher than those of A. aurita.

432
In a previous experiment, potential competitors filled a mean of 0.06 of the available space after 1 week, and 433 mean 0.24, range 0.04 to 0.97 after 5 weeks (Maxatova, 2016). A positive half-normal N (0, 0.05) prior for a 0 , 434 positive half-normal N (0, 1) for a 1 and negative half-normal N (0, 0.05) for a 2 gives mean 0.06 after 1 week, and proportional cover over almost the entire range from 0 to 1 (with mean 0.42) after 5 weeks ( Figure S4c).

436
In the removal treatments, we aimed to remove half of the target taxon. However, this was done by eye, so 437 the actual proportion removed may differ. We treated r A and r O as parameters to be estimated, with beta(2, 2) 438 priors, which are moderately concentrated around 1/2. We did not allow r A and r O to differ between depths, 439 because we have no reason to believe that the proportions removed will depend on depth.

440
In the settlement facilitation model, it seems plausible that facilitation might double the settlement rate of 441 potential competitors when polyps are very abundant (i.e. y * 1 close to 1), so m 0 could be of similar size to a 0 .

442
We therefore chose a positive half-normal N (0, 0.05) prior for m 0 . In the growth facilitation model, it seems 443 plausible that facilitation might double the proportional growth rate of potential competitors when polyps are 444 very abundant, so m 1 could be of a similar size to a 1 . We therefore chose a positive half-normal N (0, 1) prior 445 for m 1 . In the overgrowth model, it seems plausible that the rate at which potential competitors overgrow space 446 occupied by polyps might be similar to the rate at which they grow into empty space, so a 1,y * 1 could be of a 447 similar size to a 1 . We therefore chose a positive half-normal N (0, 1) prior for a 1,y * 1 . In the protection model, it  Table S1: Manova parameter estimates for the selected model (with no interaction): overall mean µ, depth effect α, removal treatment effects β 1 , β 2 , rows of lower triangle of panel covariance matrix Σ, rows of lower triangle of block covariance matrix Z. Columns are ilr coordinates. Each cell contains the posterior mean, with marginal 95 % credible highest density regions in parentheses. For some elements of Z, the highest density region consists of multiple disjoint intervals.  Table S2: Parameter estimates for ordinary differential equation models. Columns are models. Each cell contains the posterior mean, with marginal 95 % credible highest density regions in parentheses. For some parameters, the highest density region consists of multiple disjoint intervals. Where a boundary of the highest density region has a different sign from the posterior mean, this is an artefact of smoothing in density estimation, and all sample values had the same sign. All parameters other than r A and r O had separate estimated values at 1 m and 3 m. Spaces indicate parameters that were not included in a given model.  Figure S2: Choice of scale factor c for covariance matrix cI 4 of scraping treatment effects, scraping treatment by depth interaction and block effects. Each line represents a change corresponding to multiplying or dividing a single component by 2 (green), 3 (orange), 4 (purple) or 5 (pink). The y-axis is the probability of drawing a random vector at least as far from the origin as this, for scale factor c on the x-axis.
free space x potential competitors    (right, yellow). Some polyps can be seen apparently partially underneath the translucent margin of the colony. This is a closeup of the lower right corner of panel 3 (1 m depth, potential competitor removal treatment, 2019-09-24, pre-treatment, for which the whole panel is shown in Figure 2b in the main text.  Figure S16: Posterior predictive simulation from the overgrowth model for proportional cover of (a) A. aurita, (b) bare panel and (c) potential competitors. Each line is a single posterior predictive simulation for a single panel, with parameters drawn from their posterior distributions, and observations generated by multinomial sampling. Dashed lines represent panels at 1 m, and solid lines panels at 3 m. Colours represent treatments: control (C) green, A. aurita removal (A) orange, potential competitor removal (O) purple. 100 simulated data sets are shown. Lines are drawn partially transparent, so that regions with more intense colour represent higher posterior predictive density. Figure S17: Performance of the overgrowth model on simulated data sets. Each panel represents a single parameter at a single depth, as labelled on the x-axis. Green vertical lines: true parameter values (posterior mean estimates from the real data, overgrowth model). Orange dashed lines: prior densities. Thin grey lines: posterior densities from each of 10 simulated data sets simulated under the overgrowth model, with the same structure as the real data. Numbers on panels: proportion of simulated data sets for which the posterior 95% HPD interval for the parameter contained the true value. The parameters k A and k O in panels o and p are the proportions not removed, 1 − r A and 1 − k A respectively. Note that both x-and y-axis scales differ between panels.