Linking individual and population patterns of rocky-shore mussels

Study site

The study site, known as Le Petit Minou, is located in the bay of Brest, on the western coast of Brittany, France. Le Petit Minou site constitutes a sheltered macrotidal area, which presents a tidal height range of about 6 m (ca. 1–7 m height). It is made of rocky shores on either side of a beach (Fig. 1). The two shores show contrasted mean orientation, the “East shore” orientated toward the West–Southwest and the “West shore” oriented toward the Southeast (Fig. 1). Mussel patches are sparsely distributed on the West rocky shore whereas they form larger and almost continuous patches on the East rocky shore. Mussel patches consisted of only one layer of individuals in the studied area. Although populations of hybrid individuals are found in Brittany region (Bierne et al., 2003), Le Petit Minou site has been characterized by a dominance of Mytilus galloprovincialis (Simon et al., 2020).

Mussel sampling strategy

Mussel sampling was carried out on the two rocky shores to either side of Le Petit Minou beach, hereafter referred to as the East and West shores. Sampling was performed monthly from March 2019 to February 2020 by sampling three quadrats (25 × 25 cm) in each shore. The quadrats were positioned to represent different height positions, i.e., a total of 3 quadrats to represent the entire intertidal height gradient, on each shore (West and East). We sampled areas with at least 20–25% of mussel cover in low, mean and high height positions along the altitude range of mussel distribution of each shore (about 1.5 m in the study site). Quadrat location coordinates were recorded using a real-time kinematic (RTK) GPS, with Lambert93 projection. Before sampling, a picture of each quadrat was taken, which was later used to calculate the area of the quadrat covered by mussels by manually delineating the aggregation margins using ImageJ software.

A drone survey was performed with a DJI^® Phantom 4 pro in June 2019 to assess the topographic heterogeneity. Drone images were taken at a 17 m altitude with an overlap of 80% between images to optimize the image alignment and produce a relevant digital elevation model (DEM) of the area. We registered multiple reference position points from the mapped area with a RTK GPS to correct for possible distortions in the DEM reconstruction (Jaud et al., 2018). The analysis of drone images and GPS position for the photogrammetric process was performed with Agisoft PhotoScan Pro v1.2.3. Mean aspect, a circular land-surface parameter used to obtain a continuous gradient stressing the north-south or east–west gradient (Olaya, 2009) was calculated in GRASS GIS software. For each shore, aspect was represented as “northness” or “eastness” by deriving the topographic feature orientation (a sine or cosine transformation, respectively) from the DEM (Olaya, 2009). This approach allowed to extract the quadrat samples height position, northness and eastness at the precise coordinates where quadrats were sampled. Height position was referred to the Lowest Astronomical Tide (LAT; the chart datum used for tides predictions and water height measurements).

Percentage immersion time

The mean daily percentage immersion time was determined along the intertidal height gradient in the study site based on the tidal regime measurements over the preceding 5 years (from January 1, 2015 to December 31, 2019) in the study area. Hourly tide gauge measurements of water height at Brest (5 km far from Le Petit Minou) were obtained for the given period from the ‘Service hydrographique et océanographique de la Marine’ (SHOM; https://data.shom.fr/). These measurements were used to fit a cubic polynomial model to compute the mean percentage of immersion time (%) in relation to intertidal height (m). We then determined the immersion percentage of every sample based on its intertidal height position.

Sample treatment and primary data

Biomass, density, crowding and potential fecundity

It is worth mentioning that, in this study, we did not evaluate the covered area as a parameter but we evaluated the structure of the covered area inside the quadrat with diverse population parameters such as biomass, density, crowding, and others derived from biometric traits of individuals (e.g., median shell, see next section). It means these parameters do not represent the real densities or biomass of the entire site but standardized values representing the aggregate composition.

The total biomass of each individual was estimated based on the allometric relationship between SL and TWW, TWW = a SL^b. Parameters a and b were estimated by a nonlinear least square function approach applied to the ensemble of biometric observations (n individuals = 720) (R software, nls function). Then, total biomass per square meter was obtained from the sum of individual TWW (kg) and the percentage of covered area of each sample, determined from the image of the quadrat:

$\text{Total biomass}\left(\mathrm{kg}{\mathrm{m}}^{-2}\right)=\frac{\sum \mathrm{TWW}}{\mathrm{CC}\mathrm{QA}}$

where TWW is the total mass of each mussel in the quadrat, CC is the ‘proportion of covered area’, i.e., the proportion of the quadrat covered by mussels, and QA is the quadrat area = 0.0625 m².

Similarly, density was obtained by dividing the total number of individuals in a quadrat by the mussel-covered area. To examine variations in mussel crowding along the intertidal height gradient, between shores and its relationship with population density and median individual length, a novel “crowding index” was derived. It was defined as the ratio between the observed covered area, corresponding to the visually observed area from the picture analyses (Fig. 2), and the cumulative sum of the basal area of all individuals, derived from the biometric relationship: individual basal area $=\frac{L}{1.8}\frac{L}{2.4}$ (Alunno-Bruscia, Bourget & Fréchette, 2001), without including recruits. This makes it possible to represent the spatial crowding of individuals inside a patch. A ratio equal to one indicates perfect equivalence between the covered basal area and the observed covered area. Values higher than 1 indicate patches with a high degree of crowding among individuals (Fig. 2A), and values lower than 1 correspond to lower crowding and consequently greater space between individuals (Fig. 2B). Our crowding index quantifies the overlapping degree between individuals of a monolayer aggregation occurring in the study site and is original compared to other approaches (e.g., Guiñez, Castilla & Sterner, 1999).

Finally, the potential population reproductive output was estimated as the quantity of eggs released per square meter considering:

• 1:1 sex ratio,

• the mean mussel size distribution on each rocky-shore,

• the contribution of mature individuals (shell length > 12 mm, Van Haren & Kooijman, 1993; Toro, Thompson & Innes, 2002),

• a fecundity of 28% dry mass (Van der Veer, Cardoso & Van der Meer, 2006),

• and the egg quantity - mass equivalence: 10⁶ eggs = 52.5 mg (Bayne, Gabbott & Widdows, 1975).

Statistical analysis

To identify dates where the frequency of spawning individuals differed between the shores, a G-test of independence was performed between West and East rocky shore areas for each sampling date. A Chi-squared test was performed to test whether the population had a 1:1 sex ratio over the entire study site.

Statistical analyses were performed considering the independent variables: (i) West and East shores (i.e., between-shore differences), (ii) date (i.e., temporal variation), and (iii) intertidal height (in meters, continuous variable). We evaluated the effects of these variables on each population parameter and individual trait (listed in the introduction) by applying general linear models (LM) and a generalized linear model (GLM) in the case of condition index. To fulfil the assumptions of homogeneity of variance and normal distribution of residuals for linear models, a log(base10) transformation was applied to the density of individuals. The homogeneity of variance of groups was verified using the Levene test with leveneTest function (car package; Fox & Weisberg, 2019) and visual observation of the residuals versus fitted values. The normal distribution of residuals was verified by their quantile distribution and a Shapiro test of normality using shapiro.test function (stats package; R Core Team, 2020).

A GLM was performed to test the condition index (CI) using the glm function of the stats package (R Core Team, 2020), with a gamma family and the link function: log, which allows fitting count data with high dispersion. The lsmeans function (lsmeans package; Lenth, 2016) was used to perform all pairwise comparisons using a post-hoc Tukey tests with an adjusted p-value and alpha of 0.05.

The spatial variation of the instantaneous mortality rate was tested both along the intertidal height gradient and between the West and East shores, as was the interaction between these factors. This was performed by applying a general linear mixed model (LMM) with log10 transformed response variable and sampling date as a random factor using the lmer function (lme4 package; Bates et al., 2015). The effect of adult density on the recruit density was also tested in addition to the effect of shore, date and intertidal height. One LM was done with all these independent variables and a second with the significant factors, adult density and date to identify the months with significant differences.

Finally, the significance of the relationships between crowding index and the population parameters density, median individual length and shore (East/West), were analyzed using a LMM. Density of individuals and median individual length were considered as fixed factors and date as a random factor. The assumption of the normal distribution of residuals were verified with a Shapiro–Wilk normality test. For this analysis, density and median length of population size distribution were calculated without including recruits because the latter do not contribute greatly to crowding.

Intertidal height range and immersion percentage of mussel samples

The samples covered the entire intertidal height gradient on both West and East shores, ranging from 3.55 to 5.09 m, corresponding to a yearly-averaged immersion time varying from 34.47 to 61.16% (Fig. 3). The West shore showed mussel patches along almost the whole of this intertidal height range, whereas the East shore reached the lower limit of mussel distribution at a higher intertidal height (4 m). The orientation of samples on the two shores presented significant differences in eastness and northness features. The samples from the West shore were mainly orientated toward the south and almost perpendicular to the east, whereas the East shore samples were oriented to the west and almost perpendicular to the north-south axis (Fig. S1).

Condition index, spawning, sex ratio and recruitment

The sex ratio was not significantly different from 1:1 according to a Chi-squared test. The condition index (CI) was only significantly higher on the West shore than on the East shore in August, and there were some significant differences between dates depending on shore (Tables 1, S1, Fig. 4A). CI was significantly lower during March–April 2019 and February 2020 on the West shore, and during March 2019 and February 2020 on the East shore, compared with the period June–September 2019 (Fig. 4A, Table S1). Furthermore, CI was not significantly related to the intertidal height (Tables 1, S2). The temporal variation of CI was consistent with the variation in reproductive stage of individuals, with significantly lower CI during the spawning periods (Fig. 4, Table S2). The proportion of individuals at the spawning stage varied through the year with two main periods of increase during March–May and September–November, with some differences between the West and East shores during the second period (Fig. 4B). Differences between the shores were revealed by G-tests in which there was nonindependence of the frequencies of spawning individuals in July (G = 4.69, p-value = 0.03), November (G = 5.3, p-value = 0.02) and February (G = 4.07, p-value = 0.04).

Recruit density indicated an almost continuous recruitment, with significant differences between dates and the same temporal pattern on both shores (Table 1, Fig. 4C). Recruit density also had a nonsignificant relationship with intertidal height, but a significant positive relationship with the density of adult individuals at both shores considered. The relationship with adult density was non-significantly different between shores, i.e., non-significant interaction (Tables 1, S2, Fig. 5). The highest recruitment occurred from April to September, with a significantly higher density in August relative to most other months. Secondary increases were also seen on the West shore in December 2019 and February 2020, although they were not statistically significant (Fig. 4C). The first recruitment periods followed the increase of the proportion of individuals at the spawning stage (March–May 2019), but the higher recruit density was observed 3 months after the end of this period (August). The second increase in the proportion of individuals in spawning stage, September–November 2019, was followed by a small (non-significant) increase in recruit density in December, in the west shore (Fig. 4).

Individual Growth

The soVBGF parameters indicate differences between the two shores, with K = 0.12, L_inf = 65.25, C = 0.42, t_s = 0.56 and φ = 2.71 in the West, and K = 0.19, L_inf = 62.54, C = 0.51, t_s = 0.5 and φ = 2.87 in the East (Fig. 6). Both curves are shown in Fig. 6, with the estimated parameters and t₀ equal to 0.52 and 0.67 for the West and East, respectively. Moreover, these estimated growth curves matched with the range of the age-at-length data derived from the sclerochronology analysis (Fig. 6).

Mortality

The instantaneous mortality rate (Z) varied from 0.43 to 1.43 yr⁻¹, with a mean of 0.642 (±0.170 SD) in the West and 0.919 (±0.246 SD) in the East. The result from the LMM indicated that the log-transformed mortality rate (Z) was significantly different between the two shores and was positively related to the intertidal height in both, i.e., mortality increases as intertidal height increases (Fig. 7, Tables 2, S2).

Population biomass, size structure, density and crowding

No significant difference in total biomass was found between the East and West shores or between sampling dates (Tables 1, S2). There was a significant interaction indicating a shore-dependent intertidal height effect on biomass, but this was due to two very influential points in the highest topographic position in the entire site. As these points were outside the common range of distribution between West and East, we did not consider them in the final analysis, so there was no significant intertidal height effect on the biomass (Tables 1, S2, Fig. 8A).

A significantly higher density of individuals was found on the West shore than on the East, without any differences between dates (Tables 1, S2, Fig. 8B). In addition, there was a significant intertidal height effect on the density of individuals, with the same positive trend in both rocky-shores (Tables 1, S2–S3, Fig. 8B).

The median length of individuals was significantly higher on the East shore than on the West and was significantly positively related to the intertidal height position on both shores (Tables 1, S2–S3, Fig. 8C). There was a significant higher maximum length in the East than in the West with no differences between dates. Maximum length was significantly negatively related to intertidal height (Tables 1, S2–S3, Fig. 8D).

There were no significant differences in the crowding index either between West and East or between dates (LM results, Table 1, Fig. 8E). Crowding index was positively related to the density of adult individuals and the median length of adult individuals of the population, with a lower slope in the first case (LMM results, Tables 1, S2–S3, Fig. 8F).

Individual and population fecundity

Considering the mean size distribution of the population on each shore, the potential reproductive output was very similar between the two shores: 112.92 (±42, SD) and 106.5 (±32, SD) g m⁻² in West and East, respectively (Fig. 9). This was equivalent to a total of 2.15 × 10⁹ and 2.03 × 10⁹ eggs m⁻² released per covered square meter in the West and East, respectively (Fig. 9). The mean total density of reproductively mature individuals (individuals of size > 12 mm, Toro, Thompson & Innes, 2002) were 8,819 and 7,939 ind m⁻² in the West and East, respectively.

Phenology of spawning and recruitment

Spawning and recruitment phenology determine the conditions to which larvae and recruits are exposed and will, consequently, affect population dynamics. Recruitment occurs after the spawning period, larvae development and settlement. Larvae may recruit in the same area as their progenitors, i.e., self-recruitment, presumably 2–3 weeks after the progenitors spawned. However, larvae can also disperse and settle in other areas and, in this case, the lack of self-recruitment can be observed as local asynchrony between spawning and recruitment dynamics. Other causes of asynchrony are recruitment failure or a primary recruitment on other sites than the mussel bed, as in filamentous algae (Bayne, 1964). Temperature and feeding conditions determine individual spawning and recruitment phenology and success along the coast. In general, there is one main spawning period at high latitudes (e.g., Kautsky, 1982; Antsulevich, Maximovich & Vuorinen, 1999; Thorarinsdóttir & Gunnarsson, 2003) and two periods in temperate regions (e.g., Seed, 1976; Lowe, Moore & Bayne, 1982). In our study site, spawning, although continuous, peaked in spring and autumn in agreement with prior observations in the Brittany region (Lubet, 1959; Lacroix et al., 2017). Similarly, recruitment was also continuous and generally low, with peaks in spring and summer. Thus, we observed synchrony between recruitment and spawning dynamics during spring and asynchrony during summer and autumn. During summer, the asynchrony, which was caused by a peak of recruitment without a preceding peak of spawning, could have resulted from the arrival of larvae from other sites with differing reproductive phenology. In autumn, the asynchrony, which involved a peak in spawning without a following peak in recruitment, may be related to the dispersion of larvae to other sites, recruitment failure or recruitment on other substrates not analyzed in this study. These asynchronies between recruitment and spawning highlight the influence of larval dispersion and connectivity in the population dynamics of the study site. Thomas et al. (2020) pointed out the fundamental role of connectivity between sites for population dynamics along the coast of this region. Thus, the studied population at this site may have a dynamic, undefined role in overall population connectivity, i.e., acting as a larval supplier or receiver, depending on the season and environmental conditions.

Variation in spawning phenology within the study site was also observed, suggesting different temporal dynamics of microclimate conditions between the two shores. Microclimate conditions may change due to the interactions between factors such as wind speed, wave height, solar radiation, time of low tide and angle of substrate (Denny, Miller & Harley, 2006; Seabra et al., 2011; Helmuth et al., 2011; Choi et al., 2019). These factors can be particularly contrasted between West and East shores due to their different orientation. Consequently, individuals showed temporal variations in performance between shores. For instance, both a lack of spawning during July and lower condition index during August suggest a decrease in the performance of individuals on the East shore during summer compared with the West. Summer conditions were characterized by a high maximum aerial temperature (maxima of 27 °C and 25 °C in July and August, respectively, https://www.historique-meteo.net/france/bretagne/brest/2019/07/), high SST (almost 3 °C higher than the mean of the last nine years, data not shown) and low tide periods at midday. Choi et al. (2019) estimated that higher variation in substrate temperature (ranging from 25.7 °C to 41.4 °C) could be due to co-occurring high aerial temperature and high solar radiation, which is likely what took place in our study site. Gametogenesis can be slowed by high temperatures of water (Fearman & Moltschaniwskyj, 2010) and seawater temperatures higher than 20 °C can accelerate gamete release (Múgica et al., 2015), which may explain the observed discontinuity of spawning between month, and presence of oocyte atresia (data not shown). Also, a low condition index, i.e., loss of individual weight, which is mainly related either to spawning or physiological stress periods, was observed during August (when no spawning occurred) suggesting physiological stress.

Based on these observations, we hypothesize that the observed shore differences in the reproductive stage and condition of individuals were related to a temporal change in the microclimate, specifically the high temperatures recorded. The temperature also drives the reproductive phenology at within-site scale. Our work suggests the likely role of the different shore orientations within the studied site. That is, the different shore orientations potentially create refuges during stressful conditions. The different microclimates between shores were fundamental during summer because they allowed a second spawning period on the West shore, where in general conditions individuals grow slower. This is an example of how, during fluctuating conditions, individuals living in suboptimal conditions can show a higher performance than individuals living in more optimal conditions, i.e., the “suboptimal can be optimal” (Martin & Huey, 2008), which ensures their contribution to the gametes that will make the future cohort.

Variability of traits contrasts with stability of biomass

Unexpectedly, biomass was stable in both space and time in our study. This pattern of biomass stability results from the interaction between the patterns of individual traits and population parameters and is mediated by a relatively constant crowding of individuals. The aggregation structure itself results from a trade-off between the density and size of individuals. More individuals can occupy a given area if they are smaller. This pattern has already been observed in single-layer mussel aggregations (Petraitis, 1995). Population size structure and density are affected by both individual- and population-level processes, such as growth, recruitment and mortality. Petraitis (1995) pointed out that these processes can drive the capacity of a mussel bed to keep the covered area constant. Keeping the covered area constant is also associated with keeping the crowding constant to maintain the aggregation structure of a 100% covered area. The aggregation structure can remain constant if individual mussels move and tend to segregate, which has been observed in sandy environments (Liu et al., 2014) and may play a fundamental role in the aggregation pattern observed on rocky shores. Nevertheless, movement behavior under stressful wave conditions can increase the probability of dislodgement/mortality, as suggested in a M. edulis population by Schneider et al. (2005). We hypothesize, therefore, that growth, mortality and recruitment rates can also play a role in the stability of the crowding structure of an aggregation. This is supported by the observed patterns of individual growth rate, mortality and recruitment throughout the site, between shores and along the intertidal height gradient. The higher individual growth rate observed on the East shore could compensate for the higher loss of covered area or decrease of crowding due to mortality, forming mussel aggregations of larger individuals with a lower density, as observed. Conversely, on the West shore, vacant space can be filled due to higher recruitment, despite lower growth rates forming highly dense aggregations of small individuals. A similar pattern emerged along the intertidal height gradient. The higher mortality observed at increasing intertidal height could be compensated by the higher recruitment, but not likely by the growth rate since individuals grow slower at high intertidal height due to lower proportion of time in immersion. The density-dependent recruitment process observed is fundamental for maintaining the crowding level, particularly under lower growth rate conditions. This is because, when growth is slow, several individuals will fill a determined gap caused by mortality of an individual much faster than a few individuals growing slowly. Density-dependent recruitment has already been observed in other studies (e.g., McGrorty, Goss-Custard & Clarke, 1993; Blanchette, Broitman & Gaines 2006; Dolmer & Stenalt, 2010).

For sessile species on rocky-shores, space is a limiting factor, and the dynamics of space availability are fundamental for understanding population and community dynamics. The crowding structure of the population in our study indicated a highly compacted positioning of individuals. The sum of the basal area of individuals was as much as 50% higher than the observed covered area (crowding of 1.5), which indicates that the general hypothesis that the sum of area occupied by the sum of individuals is equal to the covered area (crowding equal to 1) is not corroborated in wild populations. This assumption was considered for modeling the space occupation in monolayer mussel beds (e.g., Roughgarden, Iwasa & Baxter, 1985; Petraitis, 1995), and its violation could mean an overestimation of the occupied area and consequently the possible limitation of available space for expansion of a population. Furthermore, the aggregation structure is important for patch stability and population persistence (Dugatkin, Mesterton-Gibbonsand & Houston, 1992; Guichard et al., 2003; Liu et al., 2014; Guichard, 2017). For instance, in our study site, the aggregation structure may increase patch stability, as indicated by the longevity of some individuals, found to be around 11 years old, which may underpin the persistence of this population. Thus, representing the mechanisms underlying the spatial occupation and crowding structure of mussel aggregations in rocky environments could allow an exploration of the population response under changing conditions (e.g., Zardi et al., 2021).

Interestingly, despite the different population size distribution between the two shores, our estimates indicate similar potential reproductive output (eggs per square meter), suggesting a density-size trade-off. This is to say that high-density patches of small mussels may have a similar egg production as sparser patches of large mussels, under the assumption of similar individual reproductive output (but see Sukhotin & Flyachinskaya, 2009) and 1:1 sex ratio observed in the study site. Thus, we hypothesize that the structure of the patches would tend to maintain or carry a maximum of biomass due to trade-offs between density and size while maintaining and maximizing the reproductive output.

Environmental heterogeneity underlying individual traits and population parameters patterns

The patterns of spatial variability of individual traits and population parameters in our study site suggest that heterogeneous environmental conditions influence variation in individual- and population-level processes. In our site, shore orientation likely determined different temperatures (e.g., Helmuth, 1998; Choi et al., 2019), whereas intertidal height gradient determine the duration of emersion. During emersion time, orientation heterogeneity and intertidal height could determine the maximum temperature and the duration of exposition to such temperature, respectively. Combining the temperature and duration of exposure is fundamental in regulating the mortality by heat stress (Harley, 2008; Mislan & Wethey, 2015; Seuront et al., 2019). This could be linked to the observed patterns of increasing mortality when increasing intertidal height due to a higher time of exposure. In addition, energetic demand depends on temperature, which could explain the observed differences on growth rate.

Nevertheless, determining the main driver of the observed patterns could be difficult due to the combination of multiple environmental factors. For instance, topographic orientation combined with slope can increase the wave exposition and intensity (Guichard, Bourget & Robert, 2001; Denny et al., 2004), causing higher mortality by dislodgement (McQuaid & Lindsay, 2000) and determine the swash and potential differences in immersion/emersion time. Furthermore, individuals exposed to similar thermal stress conditions (same temperature and duration) show higher survival if they have more food available, i.e., more energy to support the stress (Fitzgerald-Dehoog, Browning & Allen, 2012). This reinforces our interpretation of the gradients observed along the intertidal height gradient. Individuals located at low height, with longer immersion times, have greater access to food and consequently lower mortality rates.