Spatiotemporal dynamics of subarctic intertidal macrobenthos: going their own ways

Abstract

Intertidal macrobenthos at the small Chernaya Bight (the White Sea) was surveyed six times during 1993–2018 in order to study spatiotemporal variability. Distributions of sediments and macrophytes were highly variable in both space and time, as were most macrofaunal community attributes. Biomass slightly increased with time, while no long-term trends were found in total abundance, diversity, or functional structure. All community attributes were patchily distributed across the beach, and their patterns were not spatially autocorrelated and poorly associated with sediment properties, but changed considerably from year to year. Temporal changes in the community composition were considerable but less substantial compared with the spatial variations. The overall dynamics of species structure did not show any regular trend-like pattern but formed quasicyclic trajectories in ordination space, with nondirectional, spatially noncorrelated fluctuations around some relatively stable state. Comparison with two other neighbouring intertidal sites, studied annually in 1987–2017, showed that macrofauna at every site had similar average biomasses and common dominant species; however, the communities maintained their specificity in structure and exhibited distinct types of dynamics. In particular, the communities demonstrated different long-term trends in total biomass and diversity and followed their own paths in dynamics, appearing as differently oriented interannual trajectories. Nine most abundant species revealed no significant among-site correlations in abundance, and only two bivalve species showed good intersite agreement in dynamics of biomass. We suggest that local benthic communities are largely influenced by site-specific environmental conditions, resulting in independent and even opposite patterns of dynamics in neighbouring localities.

Introduction

Understanding the processes and mechanisms driving the changes in communities is necessary for assessing the present state and predicting the future of particular ecosystems. This is especially important at present, when many ecosystems are globally undergoing rapid alterations due to ongoing climate changes and increasing anthropogenic pressure (Knights et al., 2017; Deb and Bailey, 2023).

Observed changes in community structure and function are often interpreted as a result of broad-scale extrinsic forces, such as direct anthropogenic disturbances or climate changes (Kröncke et al., 2011, 2019; Schückel and Kröncke, 2013; Jørgensen et al., 2019). Spatiotemporal variability is, however, a complex phenomenon, which could be driven by multiple intrinsic processes such as intra- and interspecies interactions or local environmental variations (Fortin and Dale, 2005; Barnes and Ellwood, 2012; Gerwing et al., 2016). These processes are often site-specific, proceed nonsynchronously, and thus cannot be extrapolated to larger spatial scales (Fromentin et al., 1997; Arribas et al., 2019; Azovsky and Kokarev, 2019). In particular, most authors agree that small-scale variability is a general property of marine benthic assemblages, and site-specific processes prevail in determining their dynamics (Hewitt and Thrush, 2009; Gerwing et al., 2016; Arribas et al., 2019; Azovsky et al., 2022b, 2023). The relative importance of these multiple drivers acting in different scales is a fundamental question of ecological research, but disentangling their effects is a challenging task requiring long-term studies in multiple loci (Hewitt and Thrush, 2009; Soltwedel et al., 2016; Azovsky, 2019).

In this context, the White Sea represents a suitable ‘model region’ to investigate the natural variability of marine communities. Due to isolation from the direct input of Atlantic or Arctic waters, the regime of the White Sea is largely determined by regional-scale processes and is physically resistant to climate change. Until the early 2020s, no significant trends were observed in the ice cover and seasonal pattern of water temperature (Drozdov and Usov, 2014; Naumov, 2019). Compared with high Arctic areas, climatic changes are relatively weak here and have less pronounced effects on benthic communities (Solyanko et al., 2011; Chikina et al., 2020). In contrast to many other coastal areas of the Northern Hemisphere, the anthropogenic load is minor in most parts of the sea, as well as evidence of alien species invasions.

In a series of previous publications (Naumov, 2013; Varfolomeeva and Naumov, 2013; Savchenko and Naumov, 2020; Azovsky et al., 2023), we described the structure and long-term dynamics of macrobenthic communities at two geographically proximate bights that differed in geomorphology, sediment properties, and hydrological conditions. The 31-year-long regular surveys showed that both communities were similar in species composition and functional structure but substantially differed in their dynamics, without any correlation with climatic variables. The contrasting types of dynamics were suggested to be largely determined by differences in local topography and disturbance regime.

In this paper, we present a long-term study of macrobenthic community structure from another White Sea intertidal site, relying on a series of surveys from 1993 to 2018. The objectives of the study were as follows:

1. (1) Evaluate the spatial–temporal variability in species diversity, composition, and functional structure of benthic invertebrate communities, considering (a) a spatial gradient along the intertidal zone and (b) a time period of 25 years.

2. (2) Investigate the influence of sedimental features in community spatial–temporal patterns.

3. (3) Assess the among-site consistency in long-term macrobenthic dynamics in adjacent intertidal localities.

Materials and methods

Site description

The study was carried out on the gently sloping, low-energy beach situated in the sheltered shallow Chernaya Bight (CB; Kandalaksha Gulf, the White Sea, 66°32′ N, 32°58′ E) (Figure 1A, B). From late November until May, the bight is covered by ice. In the ice-free period, the surface water temperature ranges between 7 and 23°C, and the salinity varies from 16 to 24. The tides are regular semidiurnal, with a mean tidal range of 1.8 m (maximum of 2.6 m). The intertidal zone is graduated from a 90 m wide tidal flat in the northern part, with the salt marsh at the upper horizon, to a 50–80 m wide sandy beach in the southern part. Sediments vary from very silty fine-grained sands to coarse sands. According to the morphodynamic classification of McLachlan and Dorvlo (2005), it is generally typified as a tide-modified intermediate beach with a moderately steep slope (beach index 2.2). The site was described in detail by Burkovsky (2006).

Figure 1. Maps of study sites (A, B; rectangle – CB; circle – MB and SB) and locations of sampling stations in CB (C; only the outer boundaries of sampling grids are shown).

Sample collection and processing

Five surveys were carried out from July to beginning of August in 1993, 1994, 1996, 2006, and 2018. In each survey, equally spaced sampling stations were arranged in rectangular grids parallel to the shoreline. The grids varied between years in the number of stations and configuration but always covered high-, middle-, and low-intertidal zones (and, in some cases, the supralittoral zone and upper sublittoral zone below the low-tide waterline) (Figure 1C). Despite these differences, the grids widely overlapped, and the distance between nearest-neighbour stations from different surveys usually did not exceed 3–5 m, thus allowing a reliable comparison of spatial patterns. The detailed maps of stations for each survey are presented in Supplementary Figure S1.

A single sample was taken at each station with a 0.0625 m² hand box-sampler to a depth of 15–20 cm, down to the underclay layer, and a 0.01 m² subsample was taken from a total sample. The sampling protocol differed in 2006, when the large sample and small subsample were 1 and 0.04 m², respectively. Sediment from the large sample was sieved through a 3 mm mesh, sediment from the small subsample was sieved through a 0.75 mm mesh, and the retained material was fixed in 5% buffered formalin. In the laboratory, the organisms were sorted, identified to the lowest possible taxonomic level, counted, and weighed (wet weight). Animals larger than 3.0 mm in size (bivalves, priapulids, crustaceans, most gastropods, and large polychaetes) were recorded and counted from a total sample, and the smaller forms (mud snail Peringia ulvae, small polychaetes, oligochaetes, chironomids, etc.) – in a subsample only. The abundance and biomass of each species per m² were calculated. Maritime macrophytes from the whole sample were identified and weighed.

Sediments for particle-size analysis were collected at each station except for five stations from the highest and five stations from the lowest transects in 1994. Samples were oven dried at 70°C, passed through a series of sieves and weighed, and the median particle size (MPS) and silt–clay content (as the per cent of the fraction <100 μm) were calculated.

Additionally, eight randomly located samples of 0.09 m² and 12 samples of 0.016 m² were collected from the studied area in 1997 and processed as described above. The exact coordinates of these samples were not recorded, the sampling design differed, and the data on sediment composition were not available. For all these reasons, the 1997 dataset was not used for spatial distribution analysis but was used jointly as an additional datapoint in the analysis of interannual dynamics.

For the among-site comparative analysis, we also used the data on the intertidal macrobenthos collected in two small bights in the Chupa Inlet (Kandalaksha Gulf, White Sea). These bights, the Medvezhya (steeply sloping sandy beach, 66°21′ N, 33°36′ E) and Seldyanaya (sheltered silty flat, 66°20′ N, 33°37′ E), are situated 30 km to the southeast of the CB. At each site, sampling was conducted in 1987–2017 four times a year during hydrological winter, spring, summer, and autumn. Three replicate samples were collected at four permanent stations arranged in a line across the beach from the high- to low-intertidal horizon. Overall, 48 samples were collected annually at each site, so the sampling effort was comparable to that in the CB (21–72 samples, at an average of 44 samples per year). Detailed descriptions of these sites and sampling and processing methods are provided in Naumov (2013) and Azovsky et al. (2023).

Data analysis

The resulting species list was carefully checked for validity, and taxonomic assignment was unified according to the World Registry of Marine Species database (WoRMS, https://www.marinespecies.org).

To estimate the contribution of each species in a community, we used its relative respiration rate R:

$$R_i = A_i\;B_i^{0.75} \;N_i^{0.25};$$

where B_i and N_i are biomass (g wet wt m⁻²) and abundance (ind m⁻²) of the i-th species, and A_i is the taxon-specific coefficient of respiration rate. This measure combines the abundance and biomass of each species and could be treated as a rate of population energy consumption. The rationale behind using this metric is that it is more balanced than abundance (overestimating the role of small but numerous organisms, e.g. hydrobiids or oligochaetes) or biomass (overrating the role of rare but large organisms, e.g. bivalves) (Kokarev et al., 2017; Azovsky et al., 2023). Prior to the analysis, the respiration values were standardized by sample totals to use percentages. Each species was assigned to a specific environmental position (epifaunal or infaunal modalities) and to five feeding-type modalities: surface deposit feeder, subsurface deposit feeder, suspension feeder, carnivore/omnivore (including predators and scavengers), or herbivore. The feeding modalities of the species were scored using ‘fuzzy coding’ to assess the trophic structure as percentage of different feeding modalities in total community respiration (Kokarev et al., 2017; Azovsky and Kokarev, 2019) based on data from the Marine Life Information Network (MarLIN database, available at www.marlin.ac.uk) and on information found in the literature.

Alpha diversity was estimated using the average number of taxa at a sampling point, S _sample, and Shannon's index, H (to the base e). Gamma diversity was estimated as the total number of species recorded in a year (S _total). Species accumulation curves were constructed by sample-based rarefaction to visualize different diversity components standardized by sampling effort (Gotelli and Colwell, 2011). Two beta-diversity measures were used, the Whittaker β_W = S _total/S _sample − 1 and the exponent of the species accumulation curve approximated by a power function, β_SLOPE (Scheiner, 2003).

The similarities between samples were estimated by the Bray–Curtis index. Nonparametric distance-based multivariate regression analysis (DistLM) was used to explore how much of the variation in community attributes was explained by variation in the six predictor variables: calendar year of survey, local coordinates along (x) and across (y) the beach, MPS, silt content, and total phytomass (log-transformed wet weight). Euclidean distances were used for univariate metrics, and Bray–Curtis similarity was used for multivariate metrics. The stepwise selection procedure, with Akaike's information selection criterion based on 999 permutations, was used to find the best-fit explanatory regression model (Anderson et al., 2008).

Moran's I index was calculated to inspect the spatial autocorrelation in MPS values, total biomass, and alpha-diversity measures. Distance decay in the similarity relationship (i.e. the Bray–Curtis similarity values regressed against the spatiotemporal distances between samples, Azovsky et al., 2022a) was used to analyse the spatiotemporal autocorrelation in species structure.

The significance of long-term trends was estimated by a nonparametric Mann–Kendall trend test. To analyse changes in community composition among years and sites, we used metric multidimensional scaling (mMDS ordination) based on Bray–Curtis similarity. This method allows the distances on the ordination plot to be fitted to the estimated dissimilarities and thus reflects the structure of the input resemblance matrix.

To estimate the level of among-site agreement in the interannual variations of individual species, we calculated intraclass correlation coefficients (ICCs) for the abundance and biomass of the most common species. Samples collected in the summers of 1993, 1994, 1996, 1997, and 2006 in three localities were used. At each locality, average abundance and biomass values were standardized by subtracting the mean and dividing by the standard deviation, and ICC values were calculated using a two-way mixed-effects ICC model for multiple measurements (Koo and Li, 2016). An ICC >0.6 was interpreted as good intersite agreement in species dynamics.

Maps were created by regular grid interpolation of 2D datapoints (inverse distance weighting algorithm). The statistical analyses were performed using PRIMER 7 with PERMANOVA⁺ add-on (PRIMER-E Ltd., Plymouth) and PAST 4.12 software (Øyvind Hammer, Natural History Museum, University of Oslo).

Results

Environmental settings at the monitored site

Regarding sedimentology (median grain size and mud content), the studied area was rather heterogeneous, ranging from very fine silty sands (MPS = 0.10–0.15 mm, with 80–93% of silt-clay fraction) to coarse sands (MPS = 0.50–0.59 mm, 5–2.6% of silt-clay). In general, the sediments became finer and siltier from the upper to lower horizon and from the southern beach side northwards to the tidal flat.

Comparative analysis revealed the redistribution of sediments over the site, resulting in considerable variations in spatial patterns from year to year. The maps for MPS distribution are shown in Figure 2 (first column), and the maps for silt content are shown in Supplementary Figure S2A. In 1993 and 1994, the cross-beach zonation was more clearly defined, but in 1994, the sediments were coarser (on average, MPS was 0.292 and 0.366 mm, respectively). In 1996, the sediments became on average siltier, with a strong longitudinal gradient, because the sampling grid did not cover the uppermost part of the beach this year but was extended further northwards to the marsh part. In 2006, the average MPS value did not change, but the pattern reshaped to cross-beach zonation again. Finally, in 2018, both cross-beach and alongshore gradients were strongly pronounced, resulting in a diagonally oriented pattern. The Moran's correlograms indicated a spatially structured distribution of sediments in each survey. Both MPS and silt content were significantly positively autocorrelated (I > 0, P < 0.05) at short distances up to 40–45 m and negatively autocorrelated at most longer distances (Supplementary Figure S3a). The correlograms for 1993, 1994, and 2006 can be interpreted as an indication of patchiness with a few areas of finer and coarser sediments, 60–70 m in size, while in 1996 and 2018, the downward trends in the correlograms with increasing distance are indicative of a gradient in the sediment composition (Fortin and Dale, 2005).

Figure 2. Distribution maps of mean particle size (mm, left column), total macrofaunal biomass (g m⁻², middle column), and diversity (Shannon's H index) (rightmost column). Area and form of maps correspond to sampling grids. The axes are coordinates in metres; vertical and horizontal zero-reference points correspond to the uppermost intertidal boundary (spring tide level) and northwesternmost sampling station, respectively. Colour calibration bars are shown at the bottom of the figure.

Vegetation cover was also unevenly distributed over the beach and variable in time (Supplementary Figure S2B), with cover being generally much higher in the seagrass bed in the lower intertidal and subtidal zones. In the lower intertidal zone, vegetation coverage and total phytomass decreased gradually from 1460 g m⁻² (wet weight) in 1993 to 620 g m⁻² in 1994 to 300 g m⁻² in 2018; in the subtidal zone in 1994, they reached 2220 g m⁻². The vegetation in the bed was mainly composed of the eelgrass Zostera marina (80–100% of total phytomass), with the addition of the ditch grass Ruppia maritima and the green algae Cladophora fracta. In the middle intertidal zone, more or less dense vegetation of Z. marina, R. maritima, and C. fracta was observed in 1993 only, with an average phytomass of 320 g m⁻². The next year, eelgrass disappeared, and the phytomass decreased tenfold. In 1996, only a few sporadic patches of halophytic plants occurred, and later, this zone remained mostly unvegetated. In the upper shore zone, vegetation was represented by halophytic plants (R. maritima, Tripolium pannonicum, Salicornia europaea, Puccinellia sp., and Glaux maritima), which were especially abundant in the salt marsh at the northernmost part of the beach. The average phytomass in the upper zone was 200–400 g m⁻² in 1993–1997 but dropped to 49 g m⁻² in 2006 and practically disappeared in 2018.

Macrofauna: integral characteristics

Overall, 32 distinct taxa were recorded, 27 of which were identified to the species level. The richest groups were Polychaeta (ten species), Gastropoda (five), Oligochaeta (four species + some undetermined taxa), Crustacea (three), and Bivalvia (three). Two chironomid and two priapulid species were also found, as well as a few undetermined representatives of nemerteans and insect larvae.

Biomass showed a weak but statistically significant upward long-term trend (Mann–Kendall trend test, P = 0.028). No obvious long-term trends were found in total abundance or diversity measures (Table 1). The average alpha-diversity metrics (S _sample and H-index) were stable in 1993–1996 and then showed peaks in 2006 and troughs in 1997 and 2018. Beta diversity changed in phase opposite to alpha diversity, being lowest in 2006 and high in 1997 and 2018 (Table 1).

Table 1. Annual averaged community attributes in the CB

In each survey, the total biomass was patchily distributed across the beach, and the size and position of the patches changed from year to year (Figure 2, middle column). Alpha diversity (Shannon's H) was more evenly distributed, but its spatial patterns also strongly varied with time (Figure 2, rightmost column). Moreover, the among-station variations in abundance, biomass, and diversity, estimated as standard deviations, were much greater than the variations year by year, indicating that spatial variability prevailed over the interannual variability.

Although the abundances of individual taxa were patchily dispersed across the locality, the overall trophic structure of the assemblage was basically constant. In general, the community was strongly dominated by surface deposit feeders and suspension feeders, which at average accounted for 44 and 43% of total respiration, while other feeding modalities contributed much less. The proportion of epifauna to infauna in terms of respiration varied from year to year, but on average, it was 40:60. In 1994 and 1997, the infauna accounted for as much as 80% of total respiration, while at the low and high horizons in 1996 and the low horizon in 2018, the epifauna prevailed. No regular pattern in trait composition, i.e. the proportions of feeding modalities or environmental position, was observed across the beach horizons or among the years (Figure 3A, B).

Figure 3. Trait composition of community by beach horizons and years (as % of total respiration): (A) feeding modalities and (B) environmental position. Horizons labelling: high (H), middle (M), and low (L).

Main sources of variability

Nonparametric multiple regression analysis (DistLM, Table 2) showed that only for alpha diversity the predictor variables did explain a considerable part of the total variability. For other characteristics of macrofauna, only a minor part of the variability was explained by these predictors. Interannual differences were among the most relevant predictors for diversity and species composition, i.e. for characteristics that are concerned with the contribution of individual species. In contrast, for the integral community properties, such as total biomass, trophic structure, or epi- to infauna ratio, no indication of long-term trends was revealed. Vertical position across the beach (Y-coordinate) had a considerable effect on species richness (slightly increased from high to low horizon) and on total biomass, which was lowest at the high horizon. The position along the beach (X-coordinate), by contrast, had no effect on macrofauna. As a result of this spatial anisotropy, the Moran's correlograms, which operated with two-dimensional distances, did not reveal any spatial patterns in total biomass or diversity; the autocorrelation coefficients were low at all distances, except for a weak patchiness of biomass in 1994, with a 70 m distance between patches (Supplementary Figure S3b, c). Sediment properties were relevant for the functional characteristics (trophic structure or epi-/infauna ratio) while explaining only 16 and 10% of their variation. The biomass of macrophytes had no considerable effect on any community property.

Table 2. Results of DistLM analysis

Per cent of explained variation (Var Expl) and significance levels (P-values) are shown for predictors included in the best-fit models. The terms that contributed >5% are given in bold.

Species structure

The list of dominant species remained the same throughout all the observation periods, while their dominance order alternated over time. In most surveys, two species, the mud clam Mya arenaria and mud snail P. ulvae, took the first place in turn, accounting for 50–85% of total respiration. The Baltic clam Macoma balthica, blue mussel Mytilus edulis, and oligochaetes (taken together here) usually occupied the next three positions, while all other species were less abundant (Figure 4).

Figure 4. Contribution of dominant species (% of total respiration).

Regarding the dominant species, the community in 2018 was most similar to the first community state (1993). In other respects, however, the latest community showed some distinctness in composition. The abundance of oligochaetes was high in 1993–2006 (1800–3660 ind m⁻²) but dropped tenfold (to 242 ind m⁻²) in 2018; correspondingly, their role in the community decreased sharply (Figure 4). The lugworm Arenicola marina was also common in 1993–2006, with an average abundance of 7–21 ind m⁻², but only ten specimens in seven of 72 samples were found in 2018 (average abundance 0.14 ind m⁻²). Furthermore, several species that were not numerous but regularly occurred before, such as chironomids, priapulids (Priapulus caudatus, Halicryptus spinulosus), polychaetes (Eteone longa, Fabricia stellaris, Spio filicornis), and tubificid Paranais litoralis, were absent in the 2018 samples.

Some species showed certain preferences for the horizons (e.g. E. longa, oligochaetes, and chironomids were more abundant at high horizon, M. arenaria – at middle horizon, while priapulids, periwinkles, and mussels – at low intertidal). Even for these species, however, the within-horizon and temporal variability exceeded the differences among the horizons. In general, species structure was weakly related to spatial position, as revealed by DistLM analysis. We checked these spatial relationships directly, analysing the distance decay in assemblage similarity. All similarity indices were regressed against the spatial distance between samples and the time interval between surveys. The results of multiple regression analysis are presented in Supplementary Table S1. Both the effects of distance and time interval appeared to be very weak even though they were statistically significant because of the large sample size (n = 26,796). The similarity between samples decreases, on average, by 13% with a more than 100-fold increase in distance and by 3.5% only with an increasing interval between surveys from 0 to 25 years. Together, these effects explained no more than 0.3% of variations in similarity, so they should actually be considered negligible. The distance-similarity plot (Figure 5) showed the average Bray–Curtis similarity between samples as a function of spatial distance between plots and time interval between surveys. Two domains of high similarity can be distinguished: the first domain is among the samples of the same survey (zero-time interval), and the second domain is between the samples collected in 1993 and 2018 (25 year interval, upper edge of the figure); both domains form long but narrow bands covering a wide range of spatial distances, with low-similarity area between them. These results indicated the lack of remarkable autocorrelation in community spatiotemporal dynamics, except for a return to the initial state after the 25 year period.

Figure 5. Average similarity between the samples as a function of spatial distance and time lag between the surveys. Similarity values are denoted by colour (calibration bar is shown on the right).

The mMDS ordination was used to display and compare the community dynamics at different intertidal horizons. The obtained stress level was 0.1, indicating a good representation of the resemblance matrix that allowed reliable interpretation of the overall pattern. In the MDS diagram (Figure 6), the interannual trajectories showed a weak directionality and mismatched between the horizons. There was a relatively clear separation of high-intertidal horizon in 1994–2006 from other samples, while the datapoints from 1993 and 2018 clustered close to each other, independently from the horizon. The distinct position of high-intertidal samples in the middle years was due to very low abundance of M. arenaria in these samples (at average tenfold lower, as compared with middle and low horizons); this also resulted in a decreased contribution of suspension feeders to trophic structure (Figure 3A). In 2018, the population here recovered to that of the 1993 level. Thus, while some spatial differentiation of the community existed in 1994–2006 (at least between high- and mid-intertidal areas), the initial and final years sampled demonstrated the convergence of all trajectories to the same, spatially homogenized state. In general, the points representing different horizons were scattered on the ordination plot and widely overlapped, indicating pronounced interannual fluctuations but a lack of stable habitat-specific composition within the site.

Figure 6. Ordination mMDS plots based on species composition. Each point represents the community state of an intertidal horizon in a given year.

Among-site comparisons

We compared the data on intertidal macrobenthic communities in the CB with those obtained in two other sites in the Kandalaksha Gulf, Medvezhya Bight (MB), and Seldyanaya Bight (SB), studied in 1987–2017. The last two sites were described in detail earlier (Azovsky et al., 2023, and references therein). The main habitat characteristics and community attributes of these sites are presented in Supplementary Table S2.

With respect to beach geomorphology and sediment properties, CB is more similar to MB, whereas in SB, the intertidal flat is wider and more sheltered, with a muddier bottom. The three communities were rather similar in species composition. The Sörensen index of similarity of the CB species list with the MB and SB lists was 0.69 and 0.77, respectively. Only three species (Marenzelleria arctia, Cylichna alba, and Tubifex tubifex) recorded in CB were absent in MB and SB. Biomasses at all three sites were very close on average (170–210 g m⁻²), with similarly high within-site variability (Supplementary Table S2). Species richness, estimated at different levels (as the total number of taxa registered during the six surveys, per year, and per sample), was slightly but nonsignificantly lower in CB than in SB and MB. The differences in functional structure were also inessential. Epifauna and infauna were almost equally abundant in CB and SB, while in MB, epifauna prevailed. Regarding the trophic structure, surface deposit feeders were the dominant group at all sites. Suspension feeders contributed as much in the CB macrofauna but were twice as important in the MB and SB macrofauna.

At each site, species richness and diversity significantly depended on the vertical (across beach) position, increasing from high horizon seawards. For other community features, i.e. total biomass, species, and trophic structure, this factor had a weak effect. Interannual variations in diversity were also highly significant, but the pattern of change was site-specific, namely, an increasing trend in MB and a decreasing trend in SB, while no temporal trend was observed in CB. Interannual variations in biomass were also site-specific, showing a weak tendency to increase in CB, a strong upward trend with 4–5 year-long quasiperiodic oscillations in MB, and strong aperiodic fluctuations with no trend in SB. Regarding the species and trophic structure, vertical zonation contributed most in SB and MB, and the interannual changes contributed most in CB, but these effects explained only a minor part of the total variability.

Year-to-year comparison

Samples collected in the SB and MB in the summer months of 1993, 1994, 1996, 1997, 2006, 2016, and 2017 were chosen to compare the dynamics of the communities. The 2016–2017 data were used since no observations at these sites were conducted later. For each site, all samples were averaged within a year. Joint MDS ordination showed that the clusters of points did not overlap among the sites (Figure 7). Moreover, the trajectories of benthic communities were distinct from one another in shape, size, and direction. In MB, there were signs of directional shifts in community composition from 1993–1997 towards 2006–2017. This shift corresponds to a more general pattern discovered earlier for a full dataset, namely, the transition between two stable states occurred in 1993–1998 (Azovsky et al., 2023). In SB, the ordination plot did not show any clear temporal trends but looked like irregular fluctuations in different directions. The CB community showed a semiclosed trajectory, widely diverged from the starting point but returned at the end (see also Figure 6).

Figure 7. Joint mMDS ordination of community dynamics in CB, MB, and SB. Each point represents the average community state in a given year.

In the CB community, M. arenaria, P. ulvae, M. balthica, and M. edulis occupied the first place in turn, accounting together for 70–98% of community respiration; oligochaetes were also abundant here in three early surveys (1993, 1994, 1996). The same taxa also dominated in SB and, except for M. arenaria, in MB, constituting 70–90% of total respiration at both sites. A single polychaeta species, Scoloplos armiger, entered the list of dominants in MB in 1996. All these species alternated in the order of dominance from year to year, but P. ulvae was in the lead most often at every site (Supplementary Table S3). These rearrangements seem to occur at each site independently.

To estimate the among-site agreement in individual species dynamics, we calculated the ICCs for the abundance and biomass of nine common species in 1993, 1994, 1996, 1997, and 2006 across the sites. Only two bivalve species showed good intersite agreement (ICC > 0.6) in biomass dynamics: M. balthica (at all sites, the lowest biomass was recorded in 1993, and the highest biomass was recorded in 2006) and M. arenaria (high local biomasses in 1993 and 2006, low in 1996). No significant correlations were found for the abundances. The level of spatial synchrony was not related to the taxonomic position or developmental mode of a species (Table 3).

Table 3. Among-site ICCs for annual abundance and biomass values of common species

Development: P, pelagic larvae; SP, pelagic (short spawning); Bt, benthic.

^a Significant ICC (P < 0.05).

Discussion

In this study we have described a spatial–temporal variability in species diversity, composition, and functional structure of macrobenthic community from the White Sea intertidal beach, and compared the macrobenthic dynamics in this site and in the adjacent intertidal localities, in order to reveal the common trends and peculiarities. Similar trends are often interpreted as a result of broad-scale extrinsic forces (e.g. climate changes) (Kröncke et al., 2011, 2019; Schückel and Kröncke, 2013; Jørgensen et al., 2019), while the independent dynamics indicated the prevalence of local (site-specific) processes. Our analyses revealed significant variability of the community, both in space and time. Despite the similarity in the overall features (diversity, biomass, and species composition), the communities in different sites demonstrated different patterns in long-term dynamics.

General community features

In general, the macrobenthic community in the CB is typical for the soft-bottom intertidal of the White Sea and other subarctic seas. The average abundance and biomass values fell in the range of values reported for similar habitats. Overall, 32 distinct taxa (including those not identified to the species level) were recorded in this study, and this figure was also close to those reported from other studies of the White Sea soft-bottom intertidal, ranging from 30 to 45 species (Azovsky et al., 2000; Chertoprood and Azovsky, 2000; Burkovsky, 2006; Khaitov et al., 2007; Filippova et al., 2015; Naumov et al., 2018; Naumov, 2019). The only exception was the abovementioned MB and SB (75 and 73 species, respectively). Due to extensive sampling spanning a 31 year period (over 1300 samples at each site), a number of rare and incidental species were occasionally found (Azovsky et al., 2023). The composition of dominating species was also quite typical for soft-bottom intertidal habitats of the region (Naumov, 2019).

Within-site spatiotemporal patterns

Biotic zonation of the intertidal zone, appearing as a belt-like or gradient spatial arrangement of assemblages, is an established phenomenon. Known to be a distinctive and well-described phenomenon on rocky shores, vertical zonation is also often present in soft-bottom habitats (beaches and tidal flats) but is less conspicuous and more variable (Peterson, 1991). In these habitats, it is caused by species-specific responses to wave and tide regimes and sedimentology, so the grain size and silt-mud content are often reported as the most relevant environmental variables (McLachlan and Jaramillo, 1995; Degraer et al., 2008; Reise et al., 2008).

This was not the case, however, in our study. First, although some vertical gradation in sediments persisted (i.e. they generally became finer seawards), both the mapping and analysis of Moran's I correlograms detected visible changes in distribution from year to year. Furthermore, all community characteristics were highly variable in space and time. Total biomass, for instance, varied among samples almost 300 times, abundance – 175 times, number of species – tenfold. The abundance of individual species, including the dominating ones, varied even wider, resulting in relatively low among-station similarity in species structure (average value of Bray‒Curtis index was 0.44). In addition, the DistLM revealed a rather weak dependence of the community variations on spatiotemporal distances or environmental dissimilarity, except for the species richness and diversity that showed considerable temporal variations (Table 2). DistLM analysis implies a monotonous response across the whole range of predictor values (Anderson et al., 2008). Therefore, the low explanatory power of the model may indicate either a strong effect of some other unmeasured environmental factors that influenced the community, or the prevalence of small-scale, short-term variations in community descriptors. The first explanation seems less probable in our case. No significant trends were observed over the last three decades in ice cover, seasonal courses of water temperature and salinity, or tidal regime (Drozdov and Usov, 2014; Naumov, 2019). Average content of organic carbon in the CB sediments also changed negligibly, indicating the permanent rate of organic input to the beach (Azovsky et al., 2022b). The sediments are constantly redistributed across the beach by wave and tidal currents, with 10–15 year periods of local siltation or sanding (Burkovsky, 2006; Azovsky et al., 2022b; see also Figures 2 and S2A). Contrary to our expectations, however, the macrobenthic community features were poorly associated with sediment properties or vertical beach zonation. Some spatial differentiation between high and middle intertidal during middle years (1994–2006) arose from local decline in the only species, the mud clam M. arenaria. This decline is most likely resulted from long-term intrapopulation dynamics and did not affect the other species. Analysis of Moran's I correlograms and distance decay relationships indicated the lack of remarkable spatiotemporal autocorrelation in community biomass, diversity, and composition. Therefore, we suppose that the small-scale, short-term variations contributed most to the macrobenthic distribution in our case study. As a result, the distribution of macrofauna, in contrast to sediments and vegetation, looked like a random fine-grained mosaic of patches, without stable spatial structuring across the beach.

Similar patterns have also been reported for other intertidal and subtidal communities. In Plymouth Sound, UK, the mean similarity between the samples up to 50 m apart persisted more or less constant (i.e. not distance-dependent) but consistently dropped at 500 m and larger scales (Kendall and Widdicombe, 1999). Compton et al. (2008) found that the diversity and distribution of bivalve species in temperate and tropical tidal flats were not associated with sediment heterogeneity, and bivalve distributions showed a large degree of spatial overlap. They suggested that other physical and/or biological factors determined the bivalve distribution. Barnes and Ellwood (2012) studied intertidal macrofauna at multiple spatial scales. Individual species were patchily dispersed, and their assemblages were randomly structured at all studied spatial scales, but the structural variations were largest at the smallest scales (0.5–5 m). This variability is likely to result from factors such as localized recruitment of juveniles, aggregation of mobile consumers in areas of high food availability, and small-scale variation in habitat structural complexity. Gerwing et al. (2016) also found weak relationships between community variation and biotic and abiotic variables in intertidal mudflats in the Bay of Fundy, Canada. They hypothesized that the distributional patterns of intertidal communities primarily reflected stochastic small-scale spatial events. Our results are in line with all these findings.

The temporal changes in the community composition were also considerable but less substantial compared with the spatial variations. The observed community trajectories lacked any regular (trend-like) pattern but formed more or less closed loop paths in the multidimensional space (Figures 6 and 7). These changes were similarly pronounced at short (a few years) and at larger (decadal) timescales and mainly consisted of rearrangements within one and the same set of dominating species. The overall dynamics can therefore be interpreted as quasicyclic, with nondirectional, spatially noncorrelated fluctuations about some relatively stable state. The time gaps in our data after 1997 were, however, too long to reconstruct the dynamic regime with certainty.

The possible reason for these changes is the surprisingly variable distribution of sediments instead of the expected permanent zonation. Earlier, Burkovsky (2006), based on observations in the same locality from the early 1970s to 2004, suggested cyclic community fluctuations with a periodicity of approximately 22–25 years, expressed in the alternating dominance of several common species. He related these changes with the cyclic changes in silt and organic matter content in the sediments, as well as with evolution in salt marsh vegetation. Wave- and wind-generated sediment redeposition and transport may cause the changeable pattern of settlement and active migration of post-larval juveniles (Butman, 1987; Burkovsky et al., 1997). In addition, ice scouring may also have a local influence on the macrofauna. Ice observations were not conducted at this site, but severe ice-scouring events occurred irregularly in SB and profoundly affected the local community (Azovsky et al., 2023). Similar events may also be suspected at this site; the diversity drops in 1997 and 2018 and successive degradation of vegetation coverage could be the consequences of such an impact.

Among-site comparison: species-level dynamics

Spatial synchrony, i.e. the concordance of fluctuations among spatially distinct local populations of a species, has been detected in many animal and plant taxa (Azovsky et al., 2022b, and references therein). To test this type of synchrony, we used data on several common species from the three sites. Out of the nine species examined, none showed intersite consistency in abundance dynamics, and only two bivalves (M. balthica and M. arenaria) demonstrated correlated changes in biomass. Both of these species have relatively long spawning periods and long stages of pelagic larvae. The spatially synchronous population dynamics of these species have previously been reported in several locations in the White Sea intertidal zone (Maximovich and Guerassimova, 2003; Genelt-Yanovskiy et al., 2018). Spatial synchrony among local populations is often attributed to either dispersal processes (e.g. larval dispersal by along-shore currents) or regionally correlated environmental variations (the so-called Moran effect). In the above cited articles, in particular, water temperature and ice conditions were suggested as the main factors shaping the recruitment and spat survival of bivalves. Both mechanisms, however, act most of all at the earliest (larval and post-larval) stages and should therefore affect mainly abundance rather than biomass. Furthermore, dispersal-induced synchrony is expected to be strongly influenced by the developmental mode of organisms, being more pronounced for long-spawning species with pelagic larvae. These hypotheses, however, were not confirmed in our data. Quite the contrary, the intercorrelated fluctuations in biomass but not in abundance were found for both M. balthica and M. arenaria, indicating synchronization in the dynamics of the older, larger-sized cohorts. Indeed, in both localities from the Chupa Inlet, the peaks of M. arenaria biomass (in 1993 and 2006) followed the peaks in abundance (in 1988–1989 and 2002) with a 4–5 year lag. Moreover, all other tested species did not reveal any among-site synchrony, irrespective of their taxonomic position and mode of development (Table 3). Thus, intersite synchrony was revealed for a few bivalve species only and most likely reflected the correlated survival and growth rates of adults, rather than recruitment or juvenile survival. Hence, further study is required to disentangle the causes of (a)synchrony in the population dynamics of different species.

Among-site comparison: community-level dynamics

A comparative analysis showed that macrofauna at every site not only maintained their specificity in community structure but also exhibited distinct types of dynamics. Each of the three communities demonstrated different long-term trends in biomass (weak upward trend in CB, significant upward trend in MB but strong aperiodic fluctuations with no trend in SB) and in diversity (irregularly fluctuated in CB, stable or slightly increased in MB, and decreased in SB). The epifauna to infauna biomass ratio also differed: epifauna biomass was generally lower (compared with infauna) but slightly increased in CB, it was steadily higher in MB, and it decreased, becoming much lower in SB. Furthermore, each community not only maintained its specificity in structure (site trajectories occupied clearly different domains on the ordination diagram) but also followed its own path in dynamics, appearing as differently oriented trajectories (Figure 7). These trajectories could be described as quasicyclic fluctuations around a single stable state in CB, two stable states, each over 10 years long, with a relatively short transition phase in MB, while gradual degradation was observed in SB, with rapid and abrupt fluctuations occurring in years of ice scour impact (Azovsky et al., 2023). Thus, all three communities, although occupying similar habitats, have much in common (e.g. similar average biomasses and common list of dominant species), but in time, each goes its own way. This indicates that local site-specific processes prevail in determining the long-term variability.

There is sufficient evidence for such discordant dynamics of benthic communities. Subtidal assemblages, monitored over a 20 year period in four sites off northern Europe, showed different time courses (Fromentin et al., 1997). One assemblage showed large short-term oscillations, two showed long-term changes, and one remained more or less constant, and the individual species showed different trends from site to site. The changes that occurred in these four assemblages were not related to mesoscale climatic events but mainly to local environmental factors and biotic interactions. A 15 year-long study at six intertidal sandflats in Manukau Harbour, New Zealand (Hewitt and Thrush, 2009) showed that changes in the abundance of many species were asynchronous and not consistent across sites. Moreover, the sets of important location-specific predictors varied among the sites, and no variable was important across all sites for any species. Mussel bed communities in the southwestern Atlantic rocky shores also showed asynchronous dynamics at any spatial scale, even among sites located a few 100 m apart (Arribas et al., 2019). Recently, Toumi et al. (2023) applied community trajectory analysis to quantify and compare an extensive dataset on benthic coastal communities collected in a 14 year monitoring programme covering a total of 26 sites in four distinct habitats in the Northeast Atlantic. Community dynamics were not consistent across habitats and sites, and the observed trajectories were nondirectional at all monitored sites. The authors qualified this dynamic regime as the short-term, local-scale fluctuations around the site-specific equilibrium points, rather than a long-term directional trend at a regional scale.

All these examples, together with the results of this study, suggest that local benthic communities are largely influenced by site-specific environmental conditions, such as local weather patterns, hydrodynamics and sedimentation regimes, food availability, or local disturbances (e.g. ice scouring or storm actions). The joint effect of these factors may be considerable, resulting in independent and even opposite patterns of their dynamics. These local site-specific processes can shade or outweigh the common response to broader-scale factors at the regional level, especially if the response is not yet strongly pronounced. These findings highlight the caution needed in generalizing from spatially limited time-series data, which are not always able to consistently detect a response to broad-scale impacts, such as climate change or anthropogenic pressure.