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 operating at different scales (Azovsky 2000;Hewitt et al. 2007). Elucidating the particular scales of variability is essential to distinguish between different factors that influence species abundance and distribution, and determine 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 time-series 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 have 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 paid to meiofauna, which is the most abundant and diverse animal group in marine sediments (Giere 2009). High abundance and diversity, widespread distribution, and rapid generation times 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, 67 Page 2 of 12 including temporal dimensions, which has increased the difficulty 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. These are small organisms approximately 1 mm long or less, and short lived, with development times of approximately 2 or 3 weeks (Giere 2009). Similar to any other group of organisms, harpacticoids demonstrate pronounced spatiotemporal variability 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 the latter, sediment properties are usually considered the main factor that determines harpacticoid distribution and assemblage composition (Coull 1970;Williams 1972;Willems et al. 1984;Coull and Dudley 1985;Bodin and Leguellec 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 flattened and shortened stout appendages facilitate digging into the 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 flattened lanceolate bodies with short, nonprotruding appendages. These features give them the high flexibility necessary to live in the interstices between sediment particles. Phytal-dwelling species often have a stout, sometimes depressed or patelliform body with richly setose, often sturdy legs well adapted to clinging to plants. Planktonic or bentho-planktonic species have elongated streamline body and long swimming legs. Interstitial species commonly inhabit coarse-to-medium sands, epibenthic species commonly live in muddy habitats, whereas burrowing species usually prefer poorly sorted heterogeneous sediments (e.g., Bell et al. 1987;Giere 2009). However, the affinities of these groups to certain substrate types are often noticed but rarely specifically quantified.
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; see Chertoprud and Azovsky 2006 for a review). Although long-term data on harpacticoids are rare, a few examples include a 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 10-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 25-year (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 (1) to estimate the contribution of temporal (intra-and interannual) and spatial (within-and between-habitat) variations to the overall variability of the harpacticoid community and (2) to search for common long-term trends and peculiarities in harpacticoid diversity and composition in loci with different sediment properties.

Site description
The study was carried out on the gently sloping, low-energy intertidal beach situated in the semi-exposed shallow Chernaya Bay (Kandalaksha Bay, the White Sea, 66°52′ N, 32°54′ E) (Fig. 1a). Water temperatures ranged between 9 and 17 °C, and salinity ranged from 16 to 22. The tides were regular semidiurnal, with a maximum tidal range of 2.6 m.
As a part of a long-term research on benthic communities, harpacticoids were studied at three sites on the beach (A, B and C, Fig. 1b). Site A is a 50 m wide tidal flat situated in the northern part of the beach, and is protected from strong waves by numerous islands. Sediments are poorly sorted, silty fine-grained sands, with an oxic-anoxic boundary layer at a depth of 1-2 cm. Site B is situated in the central 35 m wide part of the beach, with slightly silty fine-grained sands. Site C is 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 are approximately 60-70 m 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-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 C 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 icefree 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 per site 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-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 fixed in 4% formaldehyde solution in 32 μm filtered seawater. The harpacticoid copepods were extracted by washing and sieving through a 63 μm sieve, picked out under a dissecting microscope and identified to the species or genus level in hanging-drop mounts or by dissection and enumeration. Nauplii and copepodites were not counted. To avoid possible misidentifications, part of the material collected in earlier surveys was later re-examined, and the taxonomic assignments were unified 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 weightbased 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 five 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 intra-annual or interannual dynamics, respectively.
A nonparametric Mann-Kendall trend test was used to check the significance 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 (based on the abundance data averaged over all surveys in a given year) to examine the temporal changes in species structure. A similarity profile test (SIMPROF) was applied to determine significant differences between the clusters and indicate nonrandom multivariate structures on the dendrogram . A significance level of 0.05 was chosen to stop unwarranted subdivisions. We also applied the non-metrical Multidimensional Scaling (nMDS) to ordinate the surveys and visualize the interannual trajectories in the species structure.
A distance-based nonparametric multivariate analysis of variance (PERMANOVA) was conducted to test the sources of significant differences in total abundance, diversity, species structure and life-mode composition. In the first 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 fixed factors and 'plot' as a random factor nested in the 'site'. A significant 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 five abovementioned groups (epibenthic, burrowing, interstitial, planktonic or phytal), and their percentage of total abundance was calculated. Taxa not identified 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 fine 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 2000s. In contrast, at site A, the silt content dropped from 2006 to 2018 by half again while the MPS remained steady. From 2019 to 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 ice-free 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 intra-annual dynamics of related species (e.g., five 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 significant (verified 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 significance of a positive trend was confirmed by Mann-Kendall trend test, S = 26, p = 0.003). Multivariate analysis of variance (PERMANOVA) confirmed these results ( Table 2). The seasonal differences were the only significant 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 significant, 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 significantly 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 to this figure, while the role of within-habitat variations was rather small and insignificant. 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 insignificant (p = 0.105), indicating that the harpacticoid assemblages in the three habitats differed consistently in their species structure but not in 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 confirmed the significant differences in species structure between these periods (p < 0.05) and indicated that nonrandom shifts in community structure occurred between 2013  and 2017. The results of nMDS ordination corresponded well with those of cluster analysis, indicating the considerable shifts in community structure occurred at each site between 1998-2013 and 2017-2020 periods (Online Resource 6)". 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, five 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 first 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,

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 defined in the previous section. In 1996-2013, the assemblages tangibly differed in their "ecomorphological profiles" (Fig. 7). Epibenthic harpacticoids prevailed in silty sediments at site A, and epibenthic and burrowing forms prevailed in fine 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 around 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 endopsammicdominated assemblages occurred in very fine silty sands (φ = 2.5-3.5). Burrowers are most abundant in coarse and medium sands (φ = − 0.5 to 1.5). Our data from this intertidal zone in the White Sea fitted 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
Analysis of an extensive dataset covering a 25-year period enabled us to reveal the patterns and main sources of variability for the intertidal harpacticoid community. The total abundance was relatively stable. We did not identify any systematic variations in this variable, and did not find dynamically stable spatial differences nor long-term trends. Species richness, however, demonstrated a longterm increasing trend. The between-habitat and seasonal differences, while statistically significant, 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 to 2013, while 7 species were first 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.
In total, 37 species were registered during the study. To date, diversity of harpacticoid copepods has been investigated on sandy beaches in several regions around the world, allowing for comparison of our result with others. For example, McIntyre and Murison (1973) reported 29 harpacticoid species from the small beach at the Aberdeen coast, Scotland; Coull and Dudley (1985) found 47 and 56 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 Leguellec 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); although much lower diversity values have been occasionally reported (e.g., George and Rose (2004) found only one single species during a 15-months study on a sandy beach on Chiloé Island, Chile). The faunistic results of this study are in the range of these figures; therefore, the species richness of Harpacticoida in Chernaya Bay can be evaluated as typical for beach habitats, considering the limited area of study.
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. H. minuta, H. jadensis, P. kliei, N. procerus, D. palustris and T. discipes were the most abundant species, accounting for 60-98% of the total abundance; however, the order of dominance was specific 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 pronounced clearer in the ecomorphological profiles ( Fig. 7): epibenthic species were most abundant in fine 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 Leguellec 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 finer 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 confirms this conclusion  . 8). According to these data, the transition from epibenthic-dominated to endopsammic-dominated assemblages commonly occurred in very fine 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 efficiency, 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 classification 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 fine-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 siltation of the sediments, with a decline in dominant interstitial forms, while the epibenthic species did not decline significantly. 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 10-year period on a Scottish sandy beach with relatively stable granulometric composition of sediments.
Interestingly, a similar pattern was observed in a longterm (1991-2011) study on the ciliate community carried out at the same locality close to site B Mazei 2010, 2017). The ciliate community was temporally (at an inter-year scale) stable over the first 15 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 significant, contribute little to the total variability of community structure (Table 2). Individual species also showed no definite 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 m) spatial heterogeneity. Recently, we found that after removing the effects of environmental differences, pure spatial distance decay was negligible for harpacticoids (Azovsky et al. 2022). These organisms are highly mobile and are able to leave sediments into near-bottom 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 at 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 fine 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 profiles" for a number of assemblages from geographically remote loci corroborates the generality of this pattern. Therefore, we suggest that these profiles are useful and rather universal features for describing and typifying benthic harpacticoid communities.