Environmental Stability and Long-term Variability of Harpacticoid Copepod Assemblages

Harpacticoid copepods of the Chernaya Bay (White Sea) intertidal zone were collected in 45 surveys carried out from spring to autumn over a 25-year period (1996–2020) at three sites that differed in sediment properties. There were no significant long-term trends or seasonal cycles in total abundance. Regarding the species composition, the differences between sites were the most important source of variability over the whole period while the fine-scale (within-habitat) variability was low. Epibenthic species prevailed in fine silty sand, both burrowing and epibenthic species prevailed in medium sand, and interstitial and burrowing species prevailed in coarse sand. A comparison of the data on harpacticoid assemblages from a number of geographically remote loci corroborated the generality of this pattern. In a temporal dimension, the structure of each community was stable until the early 2000s, when the proportion of epibenthic, burrowing and interstitial species changed following changes in sediment properties (increasing siltation at sandy sites and decreasing siltation at the silty site). At each site, there was an increasing long-term trend in diversity (both in total richness and in expected species number). This increase was particularly apparent at sandy sites because of the appearance of epibenthic species. We suggest that sediment composition is the key factor determining the composition of harpacticoid assemblages in space and time. The “ecomorphological profile”, i.e., the proportion of species with different lifestyle and morphological traits, is a useful and informative indicator for describing and typifying these assemblages.


Introduction
Populations and communities vary in time and space in response to multiple interacting factors and processes operated at different scales (Azovsky 2000;Hewitt et al. 2007). Elucidating the particular scales of variability is essential to distinguishing between different factors that influence species abundance and distribution and determining their roles in assembling communities. This task, however, is far from trivial and requires extensive datasets collected on a broad range of spatiotemporal scales. Despite the general importance of the problem, investigative efforts are unequally arranged across realms and taxa (Dornelas et al. 2014;Snell Taylor et al. 2018). The majority of timeseries datasets and community time series come from terrestrial bird and plant assemblages, with fewer datasets from marine and freshwater systems. Among marine datasets, the fewest gave been collected on benthic assemblages, which are generally the richest in species. Furthermore, the vast majority of these studies have focused on relatively large benthic organisms (macrofauna), while less attention has been focused on meiofauna, which is the most abundant and diverse animal group in marine sediments (Giere 2009). High abundance and diversity, widespread distribution, and rapid generation time make meiofaunal organisms excellent candidates to test general ecological hypotheses (Zeppilli et al. 2015). However, only a limited number of studies have performed direct comparisons across multiple scales, including temporal dimensions, which has increased the di culty of performing generalizations on the relationships between meiofaunal communities and environmental conditions (Rosli et al. 2018;Soltwedel et al. 2020).
Harpacticoid copepods are among the dominant groups of marine meiobenthos. They are small organisms at approximately 1 mm or less long and short-lived, with development times of approximately two or three weeks (Giere 2009). Similar to any other group of organisms, harpacticoids demonstrate pronounced spatiotemporal variability manifested at various scales. Small-scale spatial heterogeneity (cm to a few meters) is mostly attributed to sediment microtopographic features or biotic interactions (Sun et al. 1993;Azovsky et al. 2004;Smirnova and Azovsky 2020), while at larger spatial scales (tens of meters to kilometers), environmental factors prevail. Among these factors, sediment properties are usually considered a main factor that determines the harpacticoid distribution and assemblage composition (Coull 1970;Williams 1972;Willems et al. 1984;Coull and Dudley 1985;Bodin and Guellec 1992;Rybnikov et al. 2003;Chertoprud et al. 2007). Several authors have suggested that correlations occur among species lifestyle, preferred biotopes and their morphological traits, i.e., body size, shape and leg morphology (Hockin 1982;Bell et al. 1987;Giere 2009).
Among the benthic harpacticoids, several groups can be distinguished based on their morphology and lifestyle (Hicks and Coull 1983;Chertoprud et al. , 2007. Epibenthic species, which crawl on the surface of sediments, are typically large in size and have well-developed long limbs supplied with long bristles. Species that burrow through sediments have cylindrical or spindle-shaped bodies, and their attened and shortened stout appendages facilitate digging in sediment, and the bristles are partly reduced or transformed into thick prickles. Interstitial species are small and have thin, uniformly metameric, almost vermiform, or dorsoventrally attened lanceolate bodies with short, nonprotruding appendages. These features give them the high exibility necessary to live in the interstices between sediment particles. Interstitial species commonly inhabit coarse-to-medium sands, epibenthic species commonly live in muddy habitats, whereas burrowing species usually prefer poorly sorted heterogeneous sediments. However, the a nities of these groups to certain substrate types are often noticed but rarely speci cally quanti ed. In the temporal dimension, most harpacticoid studies tend to be of short duration (< 2 years) and mainly focus on seasonal dynamics. The observed seasonal patterns often differ either among species or among sites (Castel and Lasserre 1979;Heip 1979;Coull and Dudley 1985;Ansari and Parulekar 1993;Chertoprud and Azovsky 2006; see the last-named reference for review). Although long-term data on harpacticoids are rare, a few examples include 6-year study at the Grevelingen estuary (North Sea) (Willems et al. 1984), 7-year datasets collected at the Southern Bight of the North Sea (Heip 1980;Heip 1983, 1987), and 11-year datasets collected at two estuarine sites in South Carolina, USA (Coull 1985a(Coull , b, 1986Coull and Dudley 1985). In particular, Coull (1986) reported that harpacticoid assemblages in sandy and muddy sites demonstrated different temporal patterns and were apparently driven by different environmental factors. McIntyre and Murison (1973) also reported relatively constant meiofaunal compositions over a ten-year period on a Scottish sandy beach.
Although a number of studies have been performed on the spatiotemporal dynamics of meiofauna, some important questions remain unanswered. Here, we report the results of a quarter of a century-long (1996-2020) study on harpacticoid assemblages performed at three different sites in the White Sea intertidal zone. The main objectives of the present study are to (1) estimate the contribution of temporal (intra-and interannual) and spatial (within-and between-habitat) variations to the overall variability of the community and (2) search for common long-term trends and peculiarities in diversity and composition in loci with different sediment properties.

Site description
A study was carried out on the gently sloping, low-energy intertidal beach situated in the semiexposed shallow Chernaya Bay (Kandalaksha Bay, the White Sea, 66° 52' N, 32° 54' E) (Fig. 1a). Water temperatures were between 9 and 17 °C, and salinity was from 16-22. The tides were regular semidiurnal, with a maximum tidal range of 2.6 m.
As a part of long-term research on benthic communities, harpacticoids were studied at three sites on the beach (A, B and C, Fig. 1b). Site A was a 50 m wide tidal at situated in the northern part of the beach, and it was protected from strong waves by numerous islands. Sediments were poorly sorted, silty negrained sands, with an oxic-anoxic boundary layer at a depth of 1-2 cm. Site B was situated in the central 35 m-wide part of the beach, with slightly silty ne-grained sands. Site C was situated in the southern most exposed part of the beach, with well-washed medium-to-coarse sands, and the oxidized layer was usually more than 3 cm thick. The distances from sites A-B and B-C were approximately 60-70 meters each. Photographs illustrating these sites are presented in Online Resource 1 (Electronic Supplementary material).

Sampling
Sampling was performed irregularly in summer (ice-free) periods in 1996, 1998-2001, 2003, 2013, and 2017-2020. At each site, the top 5 cm of sediments were collected using plastic tube corers at several (1 to 3) 1-m 2 plots chosen randomly in the middle intertidal zone at a distance of 10-15 m apart. In the early period, the sampling design varied from year to year. In 1996, sixteen 5-cm 2 replicated samples were taken in August from each of two plots at sites A and B. In 1998, three 2-cm 2 replicated samples were taken in July from single plots at sites A, B and C. In 1999, three plots at site B and one plot at site B were repeatedly sampled on July 1, 2, 8, 15 and 31 (three 2-cm 2 replicated samples per plot). In 2000, one plot at site B was sampled on July 1, 3, 7, 12 and 25 and October 5 (three 1.8-cm 2 replicated samples in each survey). In 2001, three plots from every site were sampled during all ice-free seasons (10 surveys from May to October). In 2003, two plots from site B were sampled once in July. Later, in 2013 and 2017-2020, the sampling design was the same, with three plots from every site were sampled monthly. The sampling details are presented in Online Resource 2. Despite the differences in design, the overall sampling effort (total volume of samples per site) was similar enough (30 to 90 cm 2 , except the two earliest surveys) to provide comparability of the results across sites and years. Overall, 463 samples in 45 surveys were collected during a 25-year period.
The samples were xed in 4% formaldehyde solution in 32-μm ltered seawater. The harpacticoid copepods were extracted by washing and sieving through a 63-μm sieve, picked out under a dissecting microscope and identi ed to the species or genus level in hanging-drop mounts or by dissection and enumeration. Nauplii were not counted. To avoid possible misidenti cations, part of the material collected in earlier surveys was later re-examined, and the taxonomic assignments were uni ed across all the datasets.
To determine the grain-size distribution, sediment samples were taken from the upper 3 cm and dried. During the 1993-2000 period, sediments were collected 9, 12 and 7 times at sites A, B and C, respectively. The air-dried weight-based particle size distribution was determined via the dry sieving method; then, the mean particle size (MPS, in phi units) and silt-clay content (as the percent of the fraction < 0.1 mm) were calculated.

Data analysis
The abundance of harpacticoids was expressed as the number of individuals cm -2 . All species were grouped into ve lifestyle categories according to their morphological features (epibenthic, burrowing, interstitial, planktonic or phytal) (Chertoprud et al. , 2007. All samples from a given site were pooled together by month or by year to analyze the intra-annual or interannual dynamics, respectively. In addition to the number of species, the expected number of species ES(n) was calculated to standardize the diversity values for the different sampling efforts (Colwell et al. 2012). Based on the minimal sample counts, the randomized subsample value n of 50 or 100 specimens was chosen to analyze the intraannual or interannual dynamics, respectively.
A nonparametric Mann-Kendall trend test was used to check the signi cance of a trend.
A Bray-Curtis similarity matrix based on the relative (%) species abundances was used as a measure of similarity between sites and years. A hierarchical cluster analysis (complete linkage method) was performed separately for each site (basing on the abundance data averaged over all surveys in a given year) to examine the temporal changes in species structure. A similarity pro le test (SIMPROF) was applied to determine the signi cant differences between the clusters and indicate nonrandom multivariate structures on the dendrogram . A signi cance level of 0.05 was chosen to stop unwarranted subdivisions.
A PERMANOVA (distance-based nonparametric multivariate analysis of variance) was conducted to test the sources of signi cant differences in total abundance, diversity, species and life mode composition.
In the rst two cases, the PERMANOVA test was based on Euclidean distances, in the last two (multivariate) cases, it was based on the Bray-Curtis similarity matrices. The PERMANOVA model included 'site', 'year' and 'month' variables as xed factors and 'plot' as a random factor nested in the 'site'. A signi cant result for a given factor from PERMANOVA could signify that the groups differ either in their location (i.e., in the structure of assemblages) or in their dispersion (i.e., in the variability of assemblages) (Anderson et al. 2008). To discern these hypotheses, the PERMDISP routine was performed to test the homogeneity of multivariate dispersions.
Additionally, we used literature data to test the relationship between sediment properties and harpacticoid life-mode composition. We compiled the data on the local assemblages inhabiting intertidal and upper subtidal soft sediments. Overall, 66 quantitative datasets from 25 sources were used, listed in Online Resource 3. The species (or genera in rare instances) in each dataset were assigned into ve abovementioned groups, and their percentage of total abundance was calculated. Taxa not identi ed down to the genus level or indicated without abundance values were not considered. The MPS (in phi units) was used to characterize the sediment composition.

Sediment properties
The three sites differed distinctly in their sediment properties, from very ne silty sand (site A) to medium sand (site B) and coarse sand (site C). The sediment composition, however, changed with time ( Fig. 2). At site B, in particular, the silt content increased by 8-15% from 1993-2001 to 18-25% from 2006-2018. At site C, the silt content increased up to 10-11% from 2001-2013 compared to 2-5% before and after this period.
Thus, there is a marked tendency to silting since the early two thousandths. In contrast, at site A, the silt content dropped from 2006-2018 by half again while the MPS remained steady. From 2019-2020, the sediment composition at each site returned back to the initial state.

Abundance and diversity
In total, 37 harpacticoid species were found during the study (Table 1). Halectinosoma (7 species) and Mesochra (3 species) were the richest genera. A full species list with annual occurrences is given in the supplement (Online Resource 4). The total abundance varied widely among the samples and surveys, although the overall average abundance at the three sites was almost the same, with values of 10.5 ± 6.9, 11.0 ± 8.7 and 14.1 ± 11.4 ind/cm 2 for sites A, B and C, respectively.
No clear seasonal trends in total abundance or diversity (estimated as ES (50)) were found during the icefree period (May-October). The seasonal trends of the abundance and diversity varied between sites and years in an irregular manner (Online Resource 5). Moreover, the seasonal dynamics of individual species were also very variable, with peak abundance occurring in different months in different years. The intraannual dynamics of related species (e.g., ve Halectinosoma species, two Heterolaophonte, or two Mesochra representatives) were uncorrelated, with unmatched peaks. The only exception was Huntemannia jadensis, whose abundance permanently reached its peak in July.
In the long-term dimension, no interannual trends in abundance were observed, and the annual average abundance at each site varied around the same level of 11-13 ind/cm 2 (Fig. 3a). In contrast, both diversity measures, the total number of species per site and the expected number of species per 100 individuals, increased over time (Fig. 3b, c). This increasing trend was steepest in coarse sand (site C) and slightest in silty sediments (site A), while all trends were highly signi cant (veri ed by nonparametric Mann-Kendall trend test, p < 0.005 in all cases). The average spatial occupancy (frequency of a species occurrence in samples) also increased with time ( Fig. 3d) (data for 1998 and 2000 were not included because of the low number of samples; for other data, the signi cance of a positive trend was con rmed (Mann-Kendall trend test, S = 26, p = 0.003).
Multivariate analysis of variance (PERMANOVA) con rmed these results ( Table 2). The seasonal differences were the only signi cant source of variation for total abundance but explained only 13% of the total variability. In the case of species diversity estimated as ES(50), differences between years, months and sites were statistically signi cant, but interannual differences explained approximately half of the total variability, while the contribution of other sources was substantially smaller.

Species composition and structure of assemblages
The PERMANOVA indicated that three factors contributed signi cantly to variations in species structure (Table 2). Differences between sites and interannual changes equally contributed and explained together approximately half of the total variation. Seasonal variability added one-quarter of this gure more, while the role of within-habitat variations was rather small and insigni cant.
To estimate the role of the dispersion effects, i.e., the differences in the variability of assemblages), the PERMDISP test was performed. The overall effect of between-site differences was insigni cant (p = 0.105), indicating that the harpacticoid assemblages in the three habitats differed consistently by their species structure but not by the level of its multivariate inhomogeneity. At the same time, the dynamics of dispersion from year to year showed a distinct temporal pattern (Fig. 4). In the silty habitat (site A), high dispersion values were observed in 2001 and 2013, indicating a sharp increase in either the seasonal or microspatial heterogeneity of species structure. In the sandy habitats (sites B and C), a similar increase in dispersion occurred later, in 2017-2020.
A cluster analysis was performed separately for each site and yielded similar results (Fig. 5). On the dendrograms, two distinct periods in assemblage evolution (before and after 2017) were clearly distinguished. The SIMPROF test con rmed the signi cant differences in species structure between these periods (p < 0.05) and indicated that nonrandom shifts in community structure occurred between 2013 and 2017.
To consider species turnover, we determined whether the species occurred in only one of the two periods based on the presence-absence data (Online Resource 4). Occasional species that were found less than three times at different stations or years were not considered. Among the others, ve species were only found in 1998-2013 (Ameira parvula, Amphiascoides dispar, Vermicaris minuta, Parastenhelia spinosa, and Ectinosoma melaniceps), and all these species except E. melaniceps were interstitial. Four species appeared rst in 2017 or later, including three burrowing species (Halectinosoma chislenki, H. elongatum, and H. herdmani) and one interstitial species (Proameira hiddensoensis).
The composition of the most abundant species also changed noticeably between these periods (Fig. 6). At site A, Heterolaophonte minuta and Huntemannia jadensis dominated before 2017, and these species were replaced later by Delavalia palustris and Tachidius discipes. At sites B and C, H. jadensis and H. minuta remained among the dominant species, but after 2017, D. palustris, T. discipes and Nannopus procerus were added at site B, while at site C, Paraleptastacus kliei was replaced by N. procerus and two Nitokra species.

Life mode composition
Out of 37 species registered during the study, 12, 11, and 12 species were assigned to interstitial, burrowing, and epibenthic modes, respectively. Additionally, few specimens of one phytal species (Tisbe sp.) and one planktonic species (Microsetella norvegica) were found occasionally. We estimated the average contribution of three main lifestyle groups to the total abundance. The calculations were performed separately for three sites and for two time periods de ned in the previous section. In 1996-2013, the assemblages tangibly differed in their "ecomorphological pro les" (Fig. 7). Epibenthic harpacticoids prevailed in silty sediments at site A, and epibenthic and burrowing forms prevailed in ne sands at site B, and interstitial species were the most abundant group in coarse sands at site C. In 2017-2020, the percentage of interstitial harpacticoids increased at site A but decreased at sites B and C.
Hence, the lifestyle composition of the White Sea harpacticoid assemblages correlated with the sediment composition. Thus, we sought to determine whether this result illustrates some general pattern or a special circumstance that varies case by case, by analyzing sixty-six published datasets on the local assemblages that inhabit soft sediments from various regions all over the world. Data on the separate assemblages were ranked by MPS values and averaged within each one-unit phi interval. The results presented in Fig. 8a showed that the relative abundance of epibenthic species decreased while that of interstitial species increased in more coarse-grained sediments. The transition from epibenthic-dominated to endopsammic-dominated assemblages occurred in very ne silty sands (φ = 2.5-3.5). Burrowers are most abundant in coarse and medium sands (φ = -0.5-1.5). Our data on the White Sea intertidal zone tted this general pattern. Notably, the species richness (i.e., the total number of putative species) did not depend appreciably on the sediment composition (Fig. 8b) except for a slightly lower richness in very coarse sediments with φ < 0.

Discussion
By analyzing an extensive dataset (463 samples in 44 surveys during a 25-year period), the main patterns and sources of variability were estimated for intertidal harpacticoid communities.
In total, 37 species were registered during our study. To date, many studies have been carried out around the world to determine the diversity of harpacticoid copepods inhabiting sandy beaches. For example, McIntyre and Murison (1973) reported 29 harpacticoid species from the small beach at the Aberdeen coast, Scotland; Coull and Dudley (1985) and 73 species in an 11-year period of sampling at two sites in South Carolina. Sixteen species were present in the intertidal zone of Arcachon Bay, France (Bodin and Guellec 1992); and 25 species were present in the Mandovi estuary, west coast of India (Ansari and Parulekar 1993). Sixty species were reported by Mielke (1984) from beaches of the Galapagos Islands. An analysis of the literature showed that the average richness of local harpacticoid assemblages was generally the same in the Arctic (29 sites, 48.2 species) and in the tropics (49 sites, 46.4 species) (Azovsky et al. 2016). The faunistic results of this study are in the range of these gures; therefore, the species richness of Harpacticoida in Chernaya Bay can be evaluated as typical for beach habitats, considering the limited area of study.
The total abundance was relatively stable. We did not identify any systematic variations in this parameter or nd dynamically stable spatial differences or long-term trends. Species richness, however, demonstrated a long-term increasing trend. The between-habitat and seasonal differences, while statistically signi cant, explained only a minor part of the variation in diversity. The total species pool, however, did not increase. Overall, 11 species were found only from 1998-2013 while 7 species were rst found after 2013. Thus, the turnover was balanced and the increase in alpha diversity in recent years was not caused by the mass appearance of new species. At the same time, the mean spatial occupancy, i.e., frequency of species occurrence, doubled with time ( Fig. 3d), indicating that many species spread farther across the whole beach. This spatial expansion resulted in a higher probability of their discovery in samples from different sites and therefore was the more likely reason for the enrichment of local assemblages.
Regarding the species structure, the between-habitat differences and interannual changes were the main sources of the variation, explaining together approximately half of the total variability. Heterolaophonte minuta, Huntemannia jadensis, Paraleptastacus kliei, Nannopus procerus, Delavalia palustris and Tachidius discipes were the most abundant species, accounting for 60 to 98% of the total abundance; however, the order of dominance was speci c for each site (Fig. 6). Despite the short distance between sites (50-60 m), the assemblages maintained their compositional distinctness over a long period and demonstrate independent dynamics. These spatial and temporal structural differences are clearer pronounced in the ecomorphological pro les (Fig. 7): epibenthic species were most abundant in ne silty sand, burrowers were most abundant in medium sand, and interstitial forms were most abundant in coarse sand. Similar patterns in relation to sediment properties have been found earlier in some other case studies (Coull 1970;Williams 1972;Willems et al. 1984;Coull and Dudley 1985;Bodin and Guellec 1992;Rybnikov et al. 2003;Chertoprud et al. 2007). In particular, Moore (1979a, b) found that in the Irish Sea, well-developed interstitial harpacticoid fauna were most abundant in sands with MPSs up to 2.5-2.6 φ, whereas the group was effectively absent in ner sediments, with burrowing copepods prevailing. In the sense of Moore (1979a), the latter group also included the surface swimming (epibenthic) forms. Therefore, Moore (1979a, b) suggested 2.5-2.6 φ (165-177 µm) as the critical grain size demarcating the two types of assemblages. Our analysis based on the extensive worldwide dataset con rms this conclusion (Fig.  8). According to these data, the transition from epibenthic-dominated to endopsammic-dominated assemblages commonly occurred in very ne silty sands (according to the Wentworth grade scale, φ = 2.5-3.5).
Grain size is a key factor that determines the structure and dimensions of the sediment pore system, which is directly correlated with the anatomy of the inhabitants. In particular, both burrowing and interstitial organisms are physically restricted in their distribution by the interstitial space available to move. This factor also closely correlates with or even indirectly determines the physical and chemical milieu of the sediment, e.g., wave action, sorting e ciency, oxygen concentration or organic matter content (Giere 2009). Regardless, the distribution of life forms, distinguished by their morphology, shows an even clearer correspondence with sediment properties than the distribution of particular species. The importance of sediment properties in determining the status of benthic communities is well known, and therefore, this result might seem to be trivial. However, we emphasize it as a relevant remark since this classi cation still has not been properly documented or widely adopted in meiobenthology.
In the long-term dimension, all surveys at each site were grouped into two clusters by date: the periods of 1996-2013 and 2017-2020 (Fig. 5). These results indicated that considerable shifts in community structure occurred between 2013 and 2017. In the later period, the relative abundance of interstitial forms increased at site A but decreased at sites B and C while epibenthic and burrowing forms showed the opposite trend (Fig. 7). Most species that disappeared after 2013 were interstitial, while most species not recorded before 2017 were burrowing, indicating a compositional shift from interstitial to burrowing fauna.
These changes were evidently codirected with changes in sediment composition induced by alterations to the sedimentation regime and redistribution of the ne-grain fractions across the beach. From 2001 to approximately 2013-2018, the silt-clay content decreased in the site A sediments but then reverted back ( Fig. 2b). At sites B and C, on the contrary, the increasing siltation led to doubling of the silt-clay content, which also returned back to the initial level in recent years. These two processes, however, proceeded nonsynchronously. The changes in assemblage structure followed the changes in sediment properties with a time lag of 5-10 years. Therefore, a long response delay is quite unexpected because harpacticoids are considered a fast-responding group. This internal inertia of response likely indicates that local meiobenthic communities exhibit self-sustaining units, which maintain their structure until some threshold will be reached in changing physical conditions. Earlier, Coull (1985a, b) and Coull and Dudley (1985) found that a similar shift in community composition occurred after sharp silti cation of the sediments, with a decline in dominant interstitial forms, while the epibenthic species did not decline signi cantly. At another site with no sedimentological or apparent physical changes, no apparent changes in species diversity or composition were observed. McIntyre and Murison (1973) also reported a relatively constant meiofaunal composition over a ten-year period on a Scottish sandy beach not subjected to sedimentological dynamics.
Interestingly, a similar pattern was observed in a long-term (1991-2011) study on the ciliate community carried out at the same locality close to site B (Burkovsky & Mazei 2010, 2017. The ciliate community was temporally (interyear scale) stable over the rst fteen years but drastically transformed later. Clearly, these tendencies began in 2005 and resulted in a change in the dominant species. The authors related these transformations to increases in the silt content in the sediments, which also induced changes in the microalgal composition, an increasing amount of dead organic matter and a decrease in the oxygen content. Thus, the protistan community demonstrated a faster and sharper response to environmental changes than the harpacticoids.
Seasonal variations, while statistically signi cant, contribute little to the total variability of community structure (Table 2). Individual species also showed no de nite seasonal patterns. Huntemannia jadensis was the only exception annually and showed a peak in July, which is consistent with our earlier observations (Chertoprud and Azovsky 2006). Most of the other common species were present throughout the summer but reached a peak in different months in different years, most often from July to September, when minimal daily temperatures were permanently above zero. Such variable temporal patterns with unstable positions of the peaks are consistent with the data observed in other areas. The majority of harpacticoids are spring and summer breeders, although the reproductive patterns of species are indistinct or variable (different between years) and seem to depend on locality (Castel and Lasserre 1979;Huys et al. 1986;Webb and Parsons 1992;Chertoprud and Azovsky 2006). In particular, time series (spectral) analyses of common species at two South Carolina sites indicated no long-term population cycles with periodicity greater than 1 year. Both mud and sand sites had more year-to-year variability than seasonal variability (Coull 1985a;Coull and Dudley 1985). These authors suggested that the mechanisms controlling diversity differed in different habitats, i.e., varied microhabitats in sand and distinct seasonal suites of species in mud. Indeed, according to our data, the silty site was the only site that showed a more or less pronounced seasonal pattern in diversity, with maxima in July (Online Resource 5).
Within-habitat differences (between replicated samples within a site) counted very little in total variability, indicating low small-scale (1-10 meters) spatial heterogeneity. Recently, we found that after removing the effects of environmental differences, pure spatial distance decay was negligible for harpacticoids (Azovsky et al., in press). These organisms are highly mobile and are able to leave sediments into nearbottom water and move there a distance of tens of meters during a rising tide (Kern 1990;Walters and Bell 1994). As a result, a relatively homogeneous part of the intertidal beach constitutes a single whole biotope for most harpacticoids.

Conclusions
Harpacticoid assemblages inhabiting three spatially close sites with different granulometries exhibited clearly distinct ecomorphotypical structures, and this distinctness was maintained during the 25-year period of observations. The assemblage composition at each site was stable for a long (at least decennial) time while the habitat was unchanging but evolved rapidly after the sediment properties changed. It can be concluded that major variations in the meiofaunal community are controlled physically, mainly by the sediment composition. Epibenthic species prevailed in ne silty sand, both burrowing and epibenthic species prevailed in medium sand, while interstitial and burrowing species prevailed in coarse sand. A comparison of the "ecomorphological pro les" for a number of assemblages from geographically remote loci corroborates the generality of this pattern. Therefore, we suggest that these pro les are useful and rather universal features for describing and typifying benthic harpacticoid communities.

Declarations Statements and Declarations
Funding This work was supported by the Russian Science Foundation  Competing Interests The authors have no competing interests to declare that are relevant to the content of this article.

Ethics approval
The manuscript had not been submitted to any other journal for simultaneous consideration. The submitted work is original and was not been published elsewhere in any form or language (partially or in full).

Data availability statement
The datasets generated and analyzed during this study, if not included its supplementary information les, are available from the corresponding author on reasonable request.

Authors contribution statements
All authors contributed to the study conception and design. Material preparation and data collection were performed by Elena Chertoprud and Lesya Garlitska, the analysis was performed by Andrey Azovsky. The rst draft of the manuscript was written by Andrey Azovsky and all authors commented on previous versions of the manuscript. All authors read and approved the nal manuscript. VarExpl is the percent of total variation explained by a given factor, the p values are signi cance levels estimated by Pseudo-F test. The "-" sign indicates that the term had negative effect value and was excluded from the nal reduced model. Interannual changes in species structure (results of cluster analysis of data averaged over all surveys by a year) for sites A (a), B (b) and C (c). The nonsigni cant differences (distinguished at p > 0.05, SIMPROF test) are indicated by dashed lines Figure 6 Composition of the most abundant species at sites A (a), B (b) and C (c)

Figure 7
Page 25/26 Contribution of main life-style groups in the total abundance ("eco-morphological pro les") in 1996-2013 and 2017-2020 periods at sites A (a), B (b) and C (c) Figure 8 Average composition of life-style modes (percent of total abundance) (a) and species richness (b) in sediments with different mean particle size (MPS, phi units). Data compiled from literature