Appendix 1
Validation of the collection design
The collection design of this study had several limitations, which will now be addressed to show that they did not influence the validity of our conclusions.
(1) The material was collected once a year. This periodicity, however, is appropriate for testing our hypothesis. In the White Sea mussels usually settle in late August–early September (Dobretsov & Miron, 2001; V.Kh.’s observations). Accordingly, September, the first month of mass recruitment, can be considered as the beginning of the “biological” year for a mussel bed. Our sampling was done during the first two weeks of August, just before crucial changes in the size structure of the mussel beds were to occur. Therefore, the mussel bed demographic structure revealed in this study should reflect the results of all the changes in the population structure that have happened at the site after the previous settlement, that is, during the “biological” year of the mussel bed, which started in the September of the previous calendar year.
Special care was taken to identify primary settlers in the samples. For this, we carefully examined the external view of all mussels when measuring them under a dissecting microscope. Primary settlers are easy to identify as they have no rings of winter growth stops, only a ring marking the embryonic shell. In 18 years of studies, very few primary settlers were found. The overwhelming majority of mussels with a size lessthan 5 mm had at least 1 (usually, 3-4) rings. This means that most recruits were at least 1 year old.However, during occasional visits to the area in late August and early September primary settlers were found in large numbers, especially on filamentous algae in the subtidal. This confirms our suggestion that some fluctuations in the time of settlement did not shift the observed demographic structure.
(2) Sampling was performed with a small core (55 cm2). The choice of the core size was associated with the fact that the collected mussels were rather small (shell length 1-55 mm), and a sufficient number of them was present in a single core (Mean = 155.9 ind. core-1; SD = 111.89; Min = 20; Max = 827; total N = 314). Samplings with cores of a similar size provided relatively good estimations of abundance in other studies (Thiel & Ullrich, 2002; Tsuchiya, 2002).
However, still concerned about the small core area used, we took additional samples in 2006 and 2012. In 2006, two additional samples (by the same core) were taken near each standard core and mussel abundance from two additional cores was summarized, resulting in a doubled core area. In 2012, we additionally sampled 6 cores of the much larger area (256 cm2), randomly placed in the same way as the standard ones.
Multiplying the abundance from the standard cores and the “extra-large” ones by the corresponding coefficients (x182 for standard cores, x91 and x40 for the cores used in 2006 and 2012, respectively), the mussel abundance per 1 m2 in 6 standard cores and in 6 “extra-large” ones for each mussel bed in a given year was assessed. These data were used in two-factor ANOVA with “mussel beds” (3 levels) and “core area” (2 levels) as orthogonal fixed factors. The analyses were performed separately for 2006 and 2012. The abundance was significantly influenced by the factor “mussel bed” (F=15.44 in 2006 and F=19.9 in 2012; p<0.001 for both years) but was not influenced either by the factor “core area” (F= 2.37; p=0.134 and F=1.07; p=0.309) or by the interactions between the factors (F=2.37; p=0.111 and F=0.52; p=0.600). Thus, we may conclude that the assessment of mussel abundance did not critically depend on the core area.
(3) We sampled only mussel patches, not the entire mussel bed, which also includes mussel-free areas. As mussel population parameters are highly variable on the scale of several meters (Commito et al., 2006), the employed collection method allowed the description of only a fraction of the complex mussel bed system (the centres of large, dense patches). Therefore, our design would have been inappropriate for the assessment of the total mussel stock. This, however, was not the purpose of our study.Our aim was to assess the structure of mussel patches (in terms of mussel population density and demographic composition). At the same time, this method made it possible to standardise the samples and to minimise spatial variation in the population parameters. In this way, we could avoid the confounding effect of spatial variation, concentrating on the examination of the temporal one.
Another confounding factormight be associated with the possible spatial heterogeneity in the distribution of different size classes. However, it has been shown that the abundance of juveniles, the most variable demographic group, does not crucially depend on the position within a patch, that is, at the edgesorin the middle (Svane & Ompi, 1993).Therefore, we believe that by sampling in the middle of the patches we did not underestimate the abundance of juveniles in comparison with other mussel-covered areas inside the bed. Besides, the mussel beds in our study were very small (see above) and the distribution of juveniles across the mussel patcheswas more or less uniform in each particular year.
This collection design is unsuitable fordescribing the dynamics of some mussel bed characteristics e.g. mussel total stock.However, it is useful for describing the mussel population structure in terms of density or demographic composition. Analysing the dynamics of these values, we attempted to reveal the relative changes, i.e. the shape of thecurve describing the changes in population density and demographic structure. We hope that the dynamical properties of the time-series obtained by sampling only the central parts of mussel patches (whichseem to be representative of the entiresmall mussel beds) reflect long-term tendencies presentin other parts local populations.
Calculation of Euclidian distance matrices for presentation of long-term changes in multivariate and univariate population parameters
To represent quantitative long-term changes in mussel populations, we calculated two types of matrices, Euclidean Distances Multivariate (EDM) and Euclidean Distances Univariate (EDU) matrices. These Euclidian distances matrices were calculated separately for each mussel beds on the base of multivariate (log-transformed elements of ASCAM or elements of AMT matrix) or univariate (log-transformed yearly averaged total mussel abundances or yearly averaged regional temperature) datasets. Thereafter these matrices will be designated as EDU (triangular matrix of Euclidean distances between objects calculated on the basis of univariate data, the elements of this matrix could be calculated using Eq. 1) and EDM (the same calculated on the basis of multivariate, Eq. 2).
(Eq. 1)
where Xi and Xj – either log-transformed mean mussel abundance in years i and j or yearly averaged temperature in years i and j.
(Eq. 2)
where Yi,k and Yj,k – either log-transformed mean abundance of k-th size class in years i and j or averaged monthly air temperature in k-th month in years i and j.
We used Euclidian distances but not some other dissimilarity coefficient because the analyses based on this distance metric in the case of a singular variable produce the same or very closely related statistics as in the case of common parametric analyses such as ANOVA (Anderson, 2005) or autocorrelation (Borcard & Legendre, 2012). So, dealing with Euclidean distances, we could use a comparable logic in univariate as well as in multivariate analyses.
Cyclic matrix calculation
Briefly, the cyclic matrix is the matrix of Euclidean distances between points evenly distributed around a circle. If the time process is a stationary one (i.e. the general mean of the process is unchanged) we could imagine it as a continuous movement between these points. In this case aperfect cyclic time series could be presented as a two-dimensional array containing the coordinates of these points in the temporal order. The first coordinate could be calculated as X=cos(2πVi) and the second one as Y=sin(2πVi), where i is the number of a time point and Vi is the cyclic factor. The Vi for the cycle with the period T varies from 0 to (T-1)/T with the step 1/T. For example, for the cyclic process with the period of T=4 and number of time points I=18 the set of cyclic factors would be as follows: V1=0, V2=1/4, V3=2/4, V4=3/4, V5=0, V6=1/4 … V16=3/4, V17=0, V18=1/4. Basing on this approach, 16 cyclic matrices of different periods (T1=3…T16=18) were calculated.
PRCF analysis
Partial rate correlation function is one of the diagnostic tools allowing to reveal density dependence in the system (Berryman & Turchin, 2001). Before the analysis, all univariate time series (total mussel abundance) were detrended by subtracting the linear regression line (Legendre & Legendre, 2012). Then the log-transformed mean abundances in year (t-1) were subtracted from log-transformed mean abundances in year t, the resultant values representing the per-capita rate of change in population density, Rt (Berryman & Turchin, 2001). Further construction of a partial rate correlation function is equivalent to performing stepwise regression for the linear model
Rt=Lt-Lt-1=a0 + a1Lt-1 + a2Lt-2 +….+ adLt-d + ε ,
where Lt=lnNt, Nt – mean total mussel abundance in the year t (see Berryman & Turchin, 2001). The first step of the analysis was the calculation of correlation between Rt and Lt-1 (the first partial correlation, PRCF[1]). On the second and subsequent steps the PRCF[2]…PRCF[d] values are equal to the corresponding values of the partial autocorrelation function (Berryman & Turchin, 2001). The PRCF[1] was assessed as Pearson correlation between Rt and Lt-1 and PRCF[2] – PRCF[9] were assessed with the function “pacf()” from the package “stats”.
Literature
Anderson, M. J. , 2005. Permutational multivariate analysis of variance. Dep. Stat. Univ. Auckl. Auckl.
Berryman, A. & P. Turchin, 2001. Identifying the density-dependent structure underlying ecological time series. Oikos 92: 265–270.
Borcard, D. & P. Legendre, 2012. Is the Mantel correlogram powerful enough to be useful in ecological analysis? A simulation study. Ecology 93: 1473–1481.
Commito, J. A. , W. E. Dow & B. M. Grupe, 2006. Hierarchical spatial structure in soft-bottom mussel beds. Journal of Experimental Marine Biology and Ecology 330: 27–37.
Dobretsov, S. V. & G. Miron, 2001. Larval and post-larval vertical distribution of the mussel Mytilus edulis in the White Sea. Marine Ecology Progress Series 218: 179–187.
Legendre, P. & L. Legendre, 2012. Numerical ecology. Third English edition. Elsevier.
SvaneI. & M. Ompi, 1993. Patch dynamics in beds of the blue mussel Mytilus edulis L.: effect of site, patch size, and position within a patch. Ophelia 37: 187-202.
Thiel, M. & N. Ullrich, 2002. Hard rock versus soft bottom: the fauna associated with intertidal mussel beds on hard bottoms along the coast of Chile, and considerations on the functional role of mussel beds. Helgoland Marine Research 56: 21–30.
Tsuchiya, M. , 2002. Faunal structures associated with patches of mussels on East Asian coasts. Helgoland Marine Research 56: 31–36.