Population structure and habitat assessment for two commercial clam species exploited in small-scale fisheries

Small-scale fisheries play a crucial role in providing food and jobs in local communities worldwide. Nonetheless, their environmental impact remains poorly understood. To assess the effect of different levels of harvesting pressure on clam population dynamics, we selected three areas for study within each of three intertidal shellfish beds (in NW Spain) on the basis of historical harvesting pressure. The abundance (up to 149 ind m−2) of the introduced clam Ruditapes philippinarum of marketable size was much greater than that of the native R. decussatus (up to 20 ind m−2) in all three beds, which is consistent with the low level of recruitment of the native species. Our results suggest that the harvesting pressure did not significantly affect reproduction, which was asynchronous across beds, and that the reproductive period was longer in R. philippinarum than in R. decussatus. Nonetheless, the intertidal system was strongly affected by harvesting, as bare sediment was typical in frequently harvested areas, while sparse or dense patches of the seagrass Zostera noltei occurred in areas where harvesting pressure was scarce or null. The abundance and diversity of non-commercial infaunal species were greatest in unharvested areas. However, commercial clams were not abundant in these areas, possibly due to natural habitat preferences or artificial seeding and movement of clams along the beds. Management plans based on local knowledge of ecosystems are needed to maintain sustainable stocks of R. decussatus and also to limit the effects of harvesting on the systems.


Introduction
Coastal areas around the world support small-scale fisheries (SSFs), which mainly target sedentary and low-mobility resources (Orensanz et al. 2005) and which employ more than 90% of the world's fish workers (https:// www.fao.org).However, traditionally, greater attention has been given to industrial fisheries than to SSFs (Mills et al. 2011) and knowledge of the spatial dynamics, productivity and value of these local activities remains limited (Defeo and Castilla 2005).Global efforts to develop SSF regulations Vol:.( 1234567890) have increased in the last few decades, with the aim of preventing overexploitation of commercially exploited stocks, although the efforts have not always been successful (Selig et al. 2017).The lack of success can be attributed to the fact that management plans are often of general application and overlook local social-ecological aspects, including fishers' needs, local traditions, social cohesion, the particular characteristics of the harvesting areas, and the stocks and population dynamics of commercial and noncommercial species (De la Torre-Castro et al. 2014;Ruiz-Díaz et al. 2020).Thus, scientific knowledge about diverse social, economic and ecological aspects of populations at local scales is needed to improve management strategies and regulations and to make them acceptable to fishers.
Galicia (NW Iberian Peninsula) is the main fishing region in Spain (Freire and García-Allut 2000) and one of the most fishing-dependent regions in the European Union (Villasante et al. 2016).In the region, SSFs employ around 10,000 fishers, including 4,000 shellfishers, mainly women, working in productive intertidal areas (Villasante et al. 2021).Catches of the native clam Ruditapes decussatus (Linnaeus 1758) and the introduced clam Ruditapes philippinarum (Adams and Reeve 1850) provided total incomes of 8.5 and 17 million € respectively in 2021 (http:// www.pesca degal icia.gal).Nonetheless, populations of these species are changing rapidly due to diverse anthropogenic drivers, including climate change, overexploitation and diverse pollution events (Freire and García-Allut 2000;Loureiro et al. 2006;Domínguez et al. 2020Domínguez et al. , 2021;;Vázquez et al. 2021).A good understanding of the population dynamics of the target species is essential to improve management decisions and conservation plans.
The grooved carpet shell clam (R. decussatus) inhabits an area along the Atlantic coast ranging from the British Isles to Senegal and the Mediterranean Sea (Breber 1985).This clam was the target of intertidal shellfishing throughout centuries, until the beginning of the 1970s, when the Manila clam (R. philippinarum) was introduced to Europe for aquaculture (Flassch and Leborgne 1992).The introduced species rapidly became naturalized and is nowadays also commercially exploited (Dang et al. 2010).In fact, although catches of R. decussatus have declined significantly in Europe in the last few decades (Caill-Milly et al. 2006), catches of R.
philippinarum have increased sharply (https:// www.fao.org).This is particularly true in Galicia where catches of R. decussatus decreased from 830 t in 2000 to 278 t in 2021, while those of R. philippinarum increased from 721 t in 2000 to 2041 t in 2021 (http:// www.pesca degal icia.gal).Larger stocks of the introduced species than of the native clam may be due to greater resistance of the former to adverse environmental conditions such as those produced by climate change in recent years, e.g.heatwaves and low salinity events (Macho et al. 2016, Domínguez et al. 2020, 2021, Woodin et al. 2020, Vázquez et al. 2021, He et al. 2022).In addition, the greater abundance of R. philippinarum may be a consequence of a longer reproductive period (Genez et al. 2015), seeding of larger densities of juveniles in the shellfish beds by the shellfishers (Juanes et al. 2012), and the continual overexploitation of the more economically valuable native clam.
Harvesting pressure not only affects populations of the target species but also those of coexisting non-commercial species, as occurs with infaunal invertebrates (Kaiser et al. 1996, Griffiths et al. 2006) and with the seagrass Zostera noltei Hornemann (Cabaço et al. 2005;Román et al. 2020;Garmendia et al. 2021), which is very common in intertidal beds (e.g.Cochón andSánchez 2005, Cacabelos et al. 2015).Seagrass meadows can provide several benefits to intertidal beds, acting as a thermal buffer by providing shade or retaining water (Crespo et al. 2017;Román et al. 2022b), enhancing the deposition of particles and food availability near the bottom (Peterson et al. 1984, Gonzalez-Ortiz et al. 2014b) and decreasing the predatory pressure because the above and belowground structures hinder prey encounters and capture (Bartholomew et al. 2000;Wong 2013).Despite their importance in providing ecosystem services, seagrass habitats have declined globally in extension during the last few decades, mainly because of anthropogenic stressors such as dredging, filling and fishing activities (Lotze et al. 2006;Orth et al. 2006).The coexistence of seagrass meadows and fisheries (mainly shellfisheries) has led to the need to apply appropriate management plans to preserve the presence of meadows in commercially exploited areas.In NW Spain, there has been a conflict between shellfishing interests and seagrass conservation in the past 160 years (Fernández et al. 2022) when the seagrasses form large, dense patches that make clam harvesting difficult (Bas Ventín 2015).In a social study of Galician shellfishers, all interviewees admitted that they removed or damaged seagrasses during their daily work.The density of Z. noltei plants along the intertidal bed is thus highly dynamic and depends on the harvesting pressure (Bas Ventin 2015), resulting in different types of habitat within shellfish beds (Román et al. 2022a).
In this study, field sampling was conducted at different times in three of the most productive shellfish beds in Galicia, where shellfishing activity has modified intertidal beds during decades.The study had two main objectives: (1) to assess the population structure, including the reproductive patterns, of two commercially important clam species, R. philippinarum and R. decussatus, in relation to harvesting intensity at the shellfish bed scale, and (2) to characterize the habitat, infaunal assemblages and structural characteristics of Z. noltei meadows in relation to the harvesting intensity.Our main hypotheses were that commercial clam species would be larger and more abundant in areas with low harvesting pressure within each shellfish bed.We also expected the biodiversity of infauna assemblages within each shellfish bed to be greater in areas with low harvesting pressure.In addition, we predicted that intra-specific spatial asynchrony in reproduction would differ in relation to the presence of seagrass, as previously observed for other bivalves (Allen and Williams 2003).Finally, we predicted that such patterns would be general and consistent across shellfish beds.To our knowledge, no previous studies have tried to explain the population dynamics of commercial clams at the shellfish bed scale, while including a broad description of the ecosystem, which is crucial for developing integrated, long-term conservation plans.

Study area
The shellfish beds under study, located in the Galician Rías Baixas, NW Iberian Peninsula (Fig. 1), are among the most productive in the region.In the NW Iberian Upwelling System adjacent to the Rías Baixas, strong upwelling events mainly occur in spring and summer during northerly winds (Álvarez et al. 2005, 2008).Sampling was carried out from north to south in the Testal (Ría de Noia), O Sarrido (Ría de Arousa) and A Seca (Ría de Pontevedra) shellfish beds, which are located in intertidal areas influenced by the rivers Tambre, Ulla and Lérez respectively.The three beds are located in mesotidal estuarine environments exposed to a semidiurnal tidal regime with a long emersion period (approx.3.5 h), which leads to high temperatures in the sediment during summer.In the shellfish beds, meadows of the seagrass Z. noltei are common in the upper and mid intertidal zones, where harvesting pressure is usually less intense.Harvesting in the Testal, O Sarrido and A Seca beds is carried out by respectively 425, 211 and 438 shellfishers, mainly women (http:// www.pesca degal icia.gal).These areas are co-managed between fishers' guilds and the Galician Fisheries Authority.The exploitation and management plans are based on territorial use rights for fishing (TURFs) (Molares and Freire 2003).The management plans annually specify the authorized number of fishers, fishing areas permitted inside the fishing grounds, general objectives, state of the fishery and stock assessment analysis, harvesting and trade plans, and actions for stock enhancement.Although poaching was previously common in the shellfish beds, it is now controlled by guards hired by the fisheries associations and by the Galician Coast Guard.

Sample collection
In collaboration with technical assistants from the fishers' guilds, we selected three types of area in each shellfish bed on the basis of historical harvesting pressure: areas with no harvesting pressure in the last 5 years (hereafter NH), areas harvested only on few days per year, i.e. low harvesting pressure (hereafter LH) and, finally areas where harvesters worked monthly or weekly, i.e. intense harvesting pressure (hereafter IH).We attempted to locate these areas at the same tidal level in order to maintain the emersion period constant.However, in shellfish beds exploited in NW Spain, the low intertidal zones have been harvested continuously throughout decades and, therefore, NH and LH areas are only found in higher intertidal zones.
Sampling was carried out in the three intertidal shellfish beds on three consecutive days during low tides in the same spring tide period to prevent physiological date-dependent differences between shellfish beds.Sampling dates were November 2019, May 2020 and October 2020 (hereafter autumn 2019, spring 2020 and autumn 2020), as the beginning and end of the reproductive period occur in spring and autumn in R. decussatus and R. philippinarum (Matias et al. 2013;Moura et al. 2018) and the main recruitment period occurs in early autumn (Martínez 2009;Dang et al. 2010).Four 25 × 25 cm quadrats (0.0625 m 2 ) were placed in each of the three areas (i.e.NH, LH and IH) in each shellfish bed, and the top 10 cm of sediment was manually removed using a small shovel and sieved through a 1-mm mesh to enable determination of clam recruitment and structural characteristics of Z. noltei, when present.Simultaneously, four 50 × 50 cm quadrats (0.25 m 2 ) were delimited over the previous ones and the top 20 cm of sediment was collected with a shovel and sieved through a 5-mm mesh (see Fig. S1), for assessment of the macrobenthic assemblages, i.e. commercial clams and associated species.In addition, four surface sediment samples were selected as an example for its clarity and conspicuous sea-land thermal contrasts on regional/local spatial scales.The inset on the left-hand side shows mean sea surface temperatures between March 2012 and June 2019 retrieved from the VIIRS-SNPP satellite, represented for the Iberian coast to show large scale temperature patterns obtained in each area, with small corers (diameter 2.5 cm, depth 5 cm), to determine the particle size.Samples were placed in plastic bags, transported to the laboratory and frozen at -20 °C until further processing.

Sample processing
Zostera noltei samples were defrosted and carefully cleaned.The roots and rhizomes (belowground) were separated from the leaves (aboveground).The shoots were then counted and the length of the longest leaf on 20 shoots per sample was measured with a ruler (to the nearest 0.1 mm).The above and belowground biomass was quantified after drying the samples at 60 °C to constant weight (min.48 h).
The sieved sediment samples were examined in a stereo microscope (Leica MZ125), and small specimens (< 5 mm) of both commercial clams (R. decussatus and R. philippinarum) and of noncommercial species were counted.All small specimens of commercial clams were measured using image analysis software (Leica Application Suite V4).Specimens of other fauna larger than 5 mm were counted, and the commercial clams were also measured with a digital calliper, to the nearest 0.01 mm.The flesh of the adult-size commercial clams (≥ 25 mm) was separated from the shell and both were dried to constant weight (min 48 h at 60 ℃).The dry weights were measured with a digital balance, to an accuracy of 0.0001 g.
Sediment particle size was measured in a particle size analyser (Beckman Coulter, mod.LS13320) at the CACTI facilities (Universidade de Vigo).

Reproductive stage
In spring and autumn 2020, 20 adult individuals of each species were collected in the LH and IH areas in each shellfish bed for analysis of the gonadal development stage.NH areas were not included due to the low abundance of adult clams in these areas.As gonad is a diffuse tissue in bivalves, a piece of tissue of ca. 1 cm 2 was dissected from the foot, fixed in Davison formaldehyde for 24 h and rinsed with running water for 15 min.The sample was then dehydrated in an automatic tissue processor (Leica TP1020) and embedded in paraffin, in a paraffin embedding station (Leica EG1150H).Four sections of thickness 5 µm per slide were stained with haematoxylin and eosin (Leica RM2255 rotary microtome and Leica multistainer ST5020).The gonadal stage of each specimen was assessed under a microscope and determined by the condition of the majority of the follicles of the four sections of the slide.For both species, stages were assigned according to the unitary scale proposed by Vázquez et al. (2021), which defines seven categories (sexual rest, early gametogenesis, late gametogenesis, ripe, spawning, gonadal restoration and resorption) and was based on scales for R. decussatus (Delgado and Pérez-Camacho 2005) and R. philippinarum (Holland and Chew 1974).
The condition index (CI) was also determined in adults of both species, i.e. individuals > 25 mm in shell length (Maia et al. 2015;Moura et al. 2018), following Walne and Mann (1975):

Environmental variables
Temperature, salinity and chlorophyll a (chla) time series were recorded in the three shellfish beds between October 2019 and December 2020.The sediment temperature was recorded every 30 min by loggers (EnvLogger version 2.4, https:// elect ricbl ue.eu/) fixed to PVC poles at depths of 3 and 8 cm.Four poles, each fitted with two loggers, were placed in the sediment in each shellfish bed: two in IH and two in LH areas.Salinity was measured every 30 min by two mini-CTDs (SeaStar, DST Logic CTD v8.17, https:// www.star-oddi.com/ produ cts/ data-logge rs) located in the IH area in each shellfish bed.Two types of instruments, each measuring within a different salinity range, were used: Type I (0-18) and Type II (5-50).Daily estimates of chla concentrations (mg m −3 ) were retrieved from the IFREMER primary productivity model with a spatial resolution of 0.01° (http:// tds1.ifrem er.fr/ thred ds/ MARC-IBI-OC5_ L4-OBS/ MARC-IBI-OC5_ L4-OBS_ FULL_ TIME_ SERIE.html).Because of the limited depth and location of the shellfish beds on the shoreline, the model could not reproduce chla estimates in the study sites.Thus, the closest offshore available pixels were used for each shellfish bed at distances of 11,506 m, 2936 m and 6438 m from Testal, O Sarrido CI = dry weight of soft tissue (g) dry weight of shell (g) × 1000 Vol:. ( 1234567890) and A Seca, respectively (Fig. 1).Nevertheless, these distances in the order of few kms match the spatial resolution of satellite derived products commonly used to characterise nearshore populations (Mazzuco et al. 2015;Muñiz et al 2021;Navarrete et al. 2022).
In addition, the offshore estimates may represent good proxy values for primary productivity levels within the rias, as high levels of cross-shore water exchange flows occur due to Ekman transport (Cruz et al. 2021).To exclude salinity data during emersion and to discern temperature patterns during immersion and emersion, salinity and temperature time series were distinguished by the height at which the loggers were located relative to sea level (see Text S1 for details).

Data analysis
Four size-classes were defined in order to assess the population structure of both commercial species (see Table 1).The size classes of recruits, juveniles and adults were determined for each species on the basis of their biological characteristics, while the commercial size (40 mm and 35 mm for R. decussatus and R. philippinarum, respectively) was based on official data (https:// www.pesca degal icia.gal/).Because of the low abundance of clams of commercial size, the data were graphically plotted but added to the adult size class for statistical analysis.To investigate differences in the abundance of commercial clams within the shellfish beds, a generalized linear model (GLM) with Poisson distribution of errors for count data was constructed for each size class and shellfish bed.R. decussatus and R. philippinarum were analysed separately, as they have different biological characteristics.The models for each species therefore included two orthogonal fixed factors: Date (3 levels-autumn 2019, spring 2020 and autumn 2020) and Area (3 levels-IH, LH and NH) and their interactions.When over-or underdispersion was detected (i.e.dispersion parameter, φ ≥ 2 or ≤ 0.5, respectively), the model was corrected using a quasi-GLM model, in which the variance is given by φ x µ, where µ is the mean (Zuur et al. 2009).Overfitting was observed in two models (R. philippinarum recruits in Testal and O Sarrido), and model selection was used to reduce overfitting by removing the interaction term.The spatial variability in the reproductive stage was modelled using multinomial logistic regression, with the gonadal stage as a response variable (7 levelssexual rest, early gametogenesis, late gametogenesis, ripe, spawning, gonadal restoration and resorption).As previous studies did not reveal any differences in spawning chronology between males and females in R. decussatus (Matias et al. 2013;Genez et al. 2015) or R. philippinarum (Dang et al. 2010;Genez et al. 2015), both sexes were pooled for the analysis.The models included two fixed factors, Shellfish bed (3 levels -Testal, O Sarrido and A Seca) and Area (2 levels-IH and LH), and analyses were performed separately for each sampling date (spring and autumn 2020).Gonadal development of R. decussatus was only analysed graphically, owing to the small number of individuals collected in some cases.In addition, the spatial differences in the CI of adults of both clam species were tested separately for each species and sampling date by using Generalised Least Squares (GLS) models, including the two aforementioned orthogonal fixed factors: Shellfish bed and Area.
Time series of environmental variables were used to produce a set of predictors for gonad ripening.These predictors were the mean and the maximum sediment temperatures during emersion and immersion, the mean and minimum salinity, and the mean chla.Mean, maximum and minimum values were extracted within a period of 90 days prior to each sampling in each shellfish bed.We considered a period of 90 days a good proxy for the time during which maturity will be influenced by environmental conditions, given the length of the gametogenic cycle in these clam species (Vázquez et al. 2021).Considering the interspecific differences in the depth at which both species are found (Macho et al. 2016), the sediment temperatures at 3 and 8 cm were used to extract mean and maximum values for R. philippinarum and R. decussatus, respectively.As the salinity range of the Type 1 CTD is better for capturing reductions in salinity, this type of sensor was used to record minimum salinity, while a Type 2 sensor was used to record mean salinity.Environmental predictors were included together in binomial models with the previously defined Shellfish bed factor, to characterise the stage of maturity of both clam species after pooling samples from both IH and LH areas.For this purpose, the binomial factor Maturity (2 levels-ripe and unripe) included the ripe and spawning stages as "ripe", and sexual rest, early gametogenesis, late gametogenesis, gonad restoration and resorption as "unripe".Linear regression and analysis of variance were first used to determine the degree of collinearity between predictors.If Pearson´s R exceeded a value of 0.6 for a pair of predictors, these were not included in the same model, in order to prevent collinearity.The most complex possible models were then constructed separately for both species, including all possible interactions.A backward removal technique with a maximum of 100 steps was applied to these models, and the Bayesian Information Criterion (BIC) was used for model selection.Low BIC values indicate simpler models that explain larger amounts of variability.Unlike the Akaike´s Information Criterion, the BIC accounts for the number of observations in the model (Venables and Ripley 2002).To be considered valid, all of the terms in the selected models must be significant (see Table S1).
Morphometric characteristics of Z. noltei, i.e. shoot abundance (shoot m −2 ), shoot length (mm) and above and belowground biomass (g DW −1 m −2 ), have been widely used to characterise seagrass meadows (e.g.Attrill et al. 2000;González-Ortiz et al. 2014b;Gullström et al. 2008).We therefore integrated the four metrics using a principal component analysis (PCA) to determine the structural characteristics of each meadow and to establish the role of each variable.
The Shannon-Wiener diversity index (H') was computed to determine the diversity of the macrobenthic assemblage in each area within the shellfish beds, excluding the two commercial clam species (i.e.H' for associated species).
Although the statistical analysis was applied to the raw data, the response variables abundance and biomass were transformed according to the International System (ind m −2 and g DW −1 m −2 ) to yield values that could be compared with those obtained in other studies.The normality of the residuals and the homogeneity of variances were checked with Shapiro-Wilk tests and residuals plots, respectively.In GLS models with no homogeneous variance, the heteroscedasticity structure was specified using the varIdent function.To explore significant interaction terms or main effects included in models, post-hoc pairwise comparisons were conducted, and the p-values were adjusted with the conservative Bonferroni correction method.All data were processed using R software, version 4.0.5 (R Development Core Team 2021), and reported as means ± standard deviation (SD).The statistical significance was established at p = 0.05.Data processing and graphic representations were conducted using the tidyverse package, while the statistical analyses were performed using the nlme, vegan, nnet, emmeans, MASS packages and the default stats package.

Patterns of abundance of commercial clams
The abundance of both R. decussatus and R. philippinarum differed considerably between replicates (max SD ± 111.23 and ± 341.91 for R. decussatus and R. philippinarum populations, respectively) in the three areas characterised by different harvesting pressure, indicating enormous spatial variability at a scale of few metres (see error bars in Figs. 2 and 3).
The abundance of R. decussatus differed significantly between areas for the three size classes and shellfish beds (Table 2) except for juveniles in O Sarrido (p = 0.501).In addition, for adults in O Sarrido and the three size-classes in A Seca, the differences between areas depended on the sampling date (i.e.significant Area x Date interaction, p < 0.05, Table 2).Most of the differences indicated greater abundances in the IH than in LH and NH areas (Fig. 2); however, in some cases the abundance was greater in LH and NH areas than in the IH area, e.g. the mean abundance of juveniles in A Seca in autumn 2019 (Fig. 2).As expected, the abundance of recruits differed significantly between sampling dates in Testal and A Seca (p < 0.001, Table 2), with lower values in spring 2020 than in autumn of both years.Remarkably, only a few recruits were found in O Sarrido on the three samplings dates.By contrast, Vol:.( 1234567890) a huge peak in recruitment was observed in autumn 2020 in A Seca, where the mean abundance of recruits was 344 ± 111.23,48 ± 47.10 and 24 ± 20.65 ind m −2 in IH, LH and NH areas respectively, which revealed inter-annual differences as recruitment was almost null in autumn 2019.
All size classes of R. philippinarum were significantly more abundant in the IH areas in all three shellfish beds, except for recruits and juveniles in A Seca (Fig. 3, Table 2).In contrast to R. decussatus, R. philippinarum was scarce in the LH and NH areas in O Sarrido and A Seca.A notable abundance of adults was observed in IH areas, with mean abundances higher than 50 ind m −2 in most cases.However, the abundance differed in LH and NH areas between dates in O Sarrido and A Seca (i.e.significant Area x Date interaction, p < 0.05, Table 2).For example, the abundance was greater in the NH area (26 ± 7.58 ind m −2 ) than in the LH area (18 ± 5.97 ind m −2 ) in A Seca in autumn 2019, while it was similar in both areas in autumn 2020.As expected, the recruitment differed significantly between dates in the three shellfish beds (p < 0.05, Table 2) and was scarce in spring 2020 (i.e. 8 ± 8, 4 ± 4 and 4 ± 4 ind m −2 in Testal, O Sarrido and A Seca respectively).A high degree of inter-annual variability was observed as the maximum recruitment of R. philippinarum in A Seca was much greater in autumn 2020 (668 ± 46.96 ind m −2 ) than in autumn 2019 (64 ± 18.47 ind m −2 ).
The stock of commercially sized individuals of both species seems to be stable because the number of adults above and below the minimum harvesting size was similar, although the variability in abundance depended on the sampling date and harvesting pressure (Figs. 2 and 3).

Reproduction
In May 2020 the reproductive cycle was at an advanced stage in both R. decussatus and R. philippinarum.In all three shellfish beds, the R. decussatus individuals were at a late stage of gametogenesis or were mature, and some were spawning (Fig. 4A).Although statistical analysis was not applied to the R. decussatus data due to the small number of individuals in some samples, the graph appears to indicate spatial synchrony across the areas.
Almost the whole population of R. philippinarum in Testal was already spent, and half of the population was undergoing gonad restoration in preparation for another spawning event, while in A Seca only half and in O Sarrido only one third of the population had released gametes (Fig. 4B).Thus, the gonadal development of R. philippinarum was asynchronous across shellfish beds (LR χ 2 = 32.677,df 8, p < 0.001), mainly in A Seca and O Sarrido relative to Testal (post hoc df 12, p < 0.01), while harvesting intensity did not affect gonadal maturity (LR χ 2 = 6.053, df 4, p = 0.641).
In October 2020, the reproductive cycle of R. decussatus was already completed as the population was in the resting stage, with no or scarce follicles.In A Seca the results were uncertain because of the small number of clams collected (Fig. 4A).By contrast, around 20% of the population of R. philippinarum were still releasing gametes (Fig. 4B).In fact, although (as in May) gonadal development did not differ between areas within each shellfish bed (LR χ 2 = 0.210, df 3, p = 0.971), the shellfish bed had an important effect on the reproduction of this species (LR χ 2 = 61.996,df 6, p < 0.001).Thus, in A Seca more than 80% of the population and in Testal two thirds of the population were mature or spawning, while in O Sarrido the majority of the population Despite the small number of adult individuals of R. decussatus in the LH and NH areas (see N in bars, Fig. 5), because of which these comparisons must be considered with caution, the CI values for R. decussatus and R. philippinarum differed significantly between areas, but depended on the shellfish bed (i.e.significant interaction term Shellfish bed x Area, p < 0.001 in all cases except for autumn 2019 in R. philippinarum, Table 3).Ruditapes decussatus generally exhibited higher CI values in the three shellfish beds, reaching maxima of 130. 75 ± 47.53, 74.69 ± 12.43 and 115.54 ± 15.16 in Testal, O Sarrido and A Seca respectively, while the corresponding values for R. philippinarum were 69.81 ± 35.11, 87.82 ± 18.89 and 90.86 ± 18.21 (Fig. 5).The harvesting intensity did not have a clear influence on CI values, as the Area factor was only significant in spring 2020 for R. decussatus (p = 0.046, Table 3).

Environmental factors affecting reproduction
The environmental conditions differed between shellfish beds during 2019 and 2020.The Chla reached peaks of 20 mg cm −3 in O Sarrido and A Seca, but did not exceed 16 mg cm −3 in Testal.The sediment temperature during both emersion and immersion showed a clear seasonal pattern, and it was much lower in Testal than in O Sarrido and A Seca. Daily variations in salinity ranged between 5 and 30 and were approximately linked with tides, except in summer when the salinity was more constant (see detailed explanations in Fig. S2 and Text S2).According to the binomial analysis, the amount of variability in gonad ripening explained by the environment varied between species and depended on the specific predictors used.The best model for R. decussatus consisted of the additive combination of the maximum temperature of the sediment at 8 cm depth in emersion and shellfish bed (Figs.6A  and B), which explained ~ 36% of the total variability in gamete maturity (Table 4).As the model indicated, when the sediment temperature at a depth of 8 cm was higher than 24 °C, the capacity of R. decussatus gametes to mature decreased.By contrast, temperature was not included in the best model for R. philippinarum.The factorial combination Mean Salinity x Mean Chla was selected as the best model, although it explained only ~ 12% of the total variability in maturity; five different models explained a similar percentage of the variability (Table 4).Given the modest amounts of variability explained, the ripening stage of R. philippinarum did not seem to be determined by the environmental conditions in the different shellfish beds.Nevertheless, the ripening stage was more frequent at intermediate levels of the Mean Salinity x Mean Chla interaction (Fig. 6C, Table 4).

Characterisation of harvesting areas
The effects of different levels of harvesting pressure modified the intertidal system.In all three shellfish beds the sediment particle size was smaller in areas with no harvesting pressure (287.37 ± 31.9,Fig. 4 Percentage of A R. decussatus, and B R. philippinarum individuals at each gonadal stage, for each sampling date, shellfish bed and area.The number of clams in each category is shown inside the bars.Gonadal stages: sexual rest (0), early gametogenesis (1), late gametogenesis (2), ripe (3), spawning (4a), gonadal restoration (4b) and gonadal resorption ( 5) 470.64 ± 44.55 and 280.66 ± 39.89 µm in Testal,O Sarrido and A Seca,respectively), than in intensely harvested zones, i.e.IH areas (511.84 ± 31.8,800.84 ± 53.34 and 536.03 ± 67.62 µm in Testal,O Sarrido and A Seca,respectively).Only bare sediment was observed in IH areas, while the seagrass  Z. noltei was generally present in areas with lower harvesting pressure in all three shellfish beds.These meadows were obviously affected by the harvesting activities because even only a few days of harvesting per year reduced the shoot density from ~ 7995 shoots m −2 in NW to ~ 3405 shoots m −2 in LH areas.
PC1 explained up to 84% of the variation in the four metrics used to characterise the LH and NH areas, with the highest values corresponding to dense, stable Z. noltei meadows (Fig. 7).Aboveground biomass contributed most to PC1 (0.53), while the other metrics measured were of lower and similar importance, i.e. shoot abundance (0.49), shoot length (0.48) and belowground biomass (0.49).High values of PC1 always corresponded to unharvested areas and thus denser, more stable and continuous Z. noltei meadows were found in unharvested areas (NH); the mean aboveground biomass was 66. 92 ± 19.44, 50.74 ± 17.19 and 94.57 ± 59.83 g DW −1 m −2 in NH areas in the Testal, O Sarrido and A Seca shellfish beds, respectively.In addition, the dispersion of the points in the PCA revealed a greater heterogeneity in NH, while those of LH were closer indicating smaller variability within the areas affected by harvesting activities.For instance, in NH the shoot length varied among samples, from 195.98 ± 76.72 to 116.66 ± 42.11 mm, and the shoot abundance varied from 13,152 ± 5644.44 to 3733 ± 176.96 shoots m −2 .Belowground biomass was not the best predictor explaining the structural characteristics of meadows, but followed the same pattern as the other metrics, as it was on average 71.31, 72.46 and 69.21% lower in LH than in NH areas in Testal, O Sarrido and A Seca, respectively (Table 5).

Patterns of abundance and diversity of associated fauna
In general, the Shannon-Wiener diversity (H') values were higher in the unharvested areas, i.e. in wellestablished meadows of Z. noltei.In the Testal and O Sarrido shellfish beds, greater diversity was observed in the NH (1.15 and 1.37, respectively) than in the LH (0.85 and 0.53, respectively) and IH areas (0.76 and 0.42, respectively).By contrast, in A Seca the highest H' values occurred in IH (1.35), although the diversity was similar across the shellfish bed (Fig. 8A).
The Testal shellfish bed was dominated by the cockle C. edule, with the abundance reaching up to 232.33 ± 167.59 ind m −2 in the LH area.Likewise, there was a notable presence of the clam Scrobicularia plana (da Costa 1778) in the three areas and of the polychaete Lanice conchilega (Pallas 1766) in the LH and NH areas (Fig. 8B).Cerastoderma edule also dominated in the IH and LH areas in O Sarrido, whereas five species of bivalves and three species of polychaetes were found in the NH area, with remarkably high abundances of the bivalve Dosinia exoleta (Linnaeus 1758) (129.67 ± 24.58).However, no species dominated in A Seca, as three bivalve species coexisted in high abundance in the IH as well as in the LH and NH areas (S. plana, D. exoleta and C. edule, Fig. 8B).

Discussion
Coastal areas worldwide are highly dependent on SSFs, which are deep-rooted in society and provide food and jobs for local communities.However, most of the target resources are endangered due to several environmental and anthropogenic stressors (Selig et al. 2017;Domínguez et al. 2020Domínguez et al. , 2021)).This situation occurs in intertidal shellfish beds, estuarine habitats which are not only challenged by inherent biotic (i.e.predation, competition and grazing) and abiotic factors (i.e.warming conditions, storms, decreased salinity and tidal currents), but also by Table 4 Final selected binomial models investigating variation in the gonad ripening of R. decussatus and R. philippinarum ordered by the Bayesian Information Criterion (BIC).human harvesting pressure (Dumbauld et al. 2009).
In these areas, short-term management plans have traditionally been implanted without the approval of shellfish harvesters and only considering commercial resources, rather than the whole habitat; most of such enterprises have therefore failed (Selig et al. 2017;Ruiz-Díaz et al. 2020).This paper provides evidence of the current population structure of two commercial clams and includes an ecological overview of their habitats, which are influenced by different levels of harvesting pressure.
The sampling surveys showed that specimens of the introduced clam R. philippinarum of marketable size were much more abundant (up to 149 ind m −2 in O Sarrido) than those of the native clam R. decussatus (up to 20 ind m −2 in A Seca) in the three shellfish beds, as indicated by official landing figures (www.pesca degal icia.gal).However, the results also revealed a notable variability in a scale of few metres, as previously observed for C. edule in NW Spain (Martínez 2009), i.e. large differences between replicates in all size classes, which could be addressed in future studies with high resolution sampling focused on a single shellfish bed.The variability and population structure may be a consequence of multiple, non-exclusive factors.For example, decades of harvesting activities or artificial seeding of high densities of R. philippinarum in some of the fishing beds may account for the increase in populations of the introduced clam.We also observed that the abundance of adults of both species was similar within Z. noltei meadows (NH and LH), which are less modified by seeding and harvesting pressure, despite the generally greater abundance of R. philippinarum.Moreover, the intertidal shellfish beds are exposed to high aerial temperatures during summer emersion periods (Macho et al. 2016;Román et al. 2022b), and the effect on physiological performance, growth and reproduction is more acute in R. decussatus than in R. philippinarum (Domínguez et al. 2021;Vázquez et al.  2021; Román et al. 2022b).Thus, the thermal buffering provided by complex seagrass meadows (Crespo et al. 2017;Román et al. 2022b) may be crucial in maintaining the populations of native clams, particularly in the current context of climate change.
Seagrass meadows also act as nursery and recruitment areas for commercially important species (Herrera et al. 2022), including those of both larger-scale fisheries (LSF) and small-scale fisheries (SSF) (Unsworth et al. 2018).However, our findings indicated a generally higher abundance of recruits and juveniles of both commercial clam species (R. decussatus and R. philippinarum) in bare sediment (IH) than in areas where Z. noltei is present (LH and NH), although recruitment of R. decussatus was higher in areas characterised by low harvesting pressure in A Seca.This finding highlights the importance of carrying out studies at shellfish bed scale to characterise each area rather than making management plans based on general information.As seagrass meadows clearly have role as nursery areas for commercially important species in subtidal habitats (Bertelli and Unsworth 2014; Leslie et al. 2017) and as recruitment can be negatively affected by long emersion periods (Lamb et al. 2014), the location of unharvested areas high up in the intertidal zone may explain the low level of recruitment.Unfortunately, we could not investigate this effect in detail because of the relationship between intertidal height and harvesting effort in the shellfish beds under study.However, further studies disentangling this confounding effect could shed some light on differential recruitment rates in these habitats.
Another key aspect explaining the population structure is the reproductive cycle and modifying factors.Individuals of the introduced clam R. philippinarum were still releasing gametes in October, whereas most R. decussatus individuals were at the resting stage in all three shellfish beds.Thus, our results suggest that the introduced species had a longer reproductive period, which implies a longer recruitment period, as also indicated by Genez et al. (2015).The presence of Z. noltei did not affect the maturation of any clam species, as found for the infaunal amphipod Corophium volutator (Pallas 1766) (Boström et al. 2006).However, the role of the seagrass meadows on the reproductive cycle of invertebrates has not yet been addressed and requires further study.
Food availability is known to cause differences in the condition index (CI) of bivalves (Orban et al. 2008;Tsai et al. 2010).The CI of bivalves is usually associated with reproductive activity, although it can also refer to the flesh content relative to shell.This parameter is therefore of economic relevance for comparing the commercial quality of different species (Orban et al. 2008).In this respect, our results highlight the commercial quality of the native clam, as the CI was higher than that of the introduced species in all cases.Although previous studies clearly indicated a seasonal pattern in CI linked to reproduction (e.g.Dang et al. 2010;Hernández-Otero et al. 2014), we observed more constant values.For instance, in autumn 2020 most R. decussatus individuals were at the resting stage, but the CI was similar to that in spring; a lack of correspondence between the CI and reproduction was also observed in the Ria de Aveiro population of this species (Matias et al. 2013).In line with Dang et al. 2010, we observed differences in CI values from those in nearby shellfish beds, although the differences were not constant across areas.In our study, the presence of Z. noltei did not affect the CI, which was higher inside the meadows in some cases and outside in others.Contrasting results regarding the relationship between CI and the presence of seagrass have been reported in diverse studies.For instance, the body and gonadal mass of the mixotrophic Loripes lucinalis (Lamarck, 1818) increased when the bivalve fed on the suspended particulate organic matter retained inside seagrass meadows (van der Geest et al. 2014).By contrast, the presence of the introduced seagrass Zostera japonica Ascherson and Graebner led to a reduction in the CI of R. philippinarum in northeast Pacific Ocean (Tsai et al. 2010).
The results clearly demonstrated asynchrony in reproduction across shellfish beds separated by only ~ 40 kms, which may be a consequence of intermediate-scale shifts in environmental conditions, as previously observed for stalked barnacles along rocky shores in Galicia (Román et al. 2022c).Small-scale variability associated with shifts in environmental factors (i.e.sea water temperature and salinity) was also observed in the gametogenesis and maturation stages of the invasive mussel Xenostrobus securis (Lamarck, 1819) and the commercial razor clam Ensis magnus (Shumacher, 1817) inhabiting the Galician Rías Baixas (Hernández-Otero et al. 2014; Vol.: (0123456789) Montes et al. 2020).In this study, high sediment temperatures, which varied among shellfish beds, seemed to have a negative effect on the ripening of R. decussatus, while ripening of R. philippinarum was not notably affected by any of the factors studied.Although this finding is consistent with those of several previous laboratory studies, it must be considered with caution due to the small number of samplings.For instance, the increase in surface sediment temperature to 32 °C led to an increase in burrowing activity in R. decussatus, with an associated energetic cost (Macho et al. 2016), and to the presence of more abnormal oocytes in postspawning follicles (Vázquez et al. 2021).By contrast, He et al. (2022) observed a mortality rate of only 2% in R. philippinarum exposed to a temperature of 40 °C, confirming the ability of this species to cope with extreme thermal events as well as extreme salinity (Domínguez et al. 2020).The reduced impact of high temperature or low salinity on the ripening of the introduced species may explain the long reproductive period and the ability of this clam to colonize habitats with different environmental characteristics (Harris 2016), as occurs with other non-native invertebrate species (Gosling 2003;Montes et al. 2020).Sediment temperature was the most important variable determining reproduction of the native clam, which highlights the value of this parameter (rather than the more traditionally used seawater temperature) in explaining physiological mechanisms of intertidal species exposed to temperatures close to their thermal tolerance limits (Olabarria et al. 2016).
Regarding the effects of harvesting activity on the intertidal system, this study indicated an inverse relationship between harvesting pressure and the presence and density of seagrass meadows in the three shellfish beds under study.This obvious pattern is aligned with previous findings of reduced shoot density, biomass and total seagrass area caused by the continuous digging and trampling associated with shellfishing activities (Cabaço et al. 2005;Román et al. 2020;Garmendia et al. 2021;Barañano et al. 2022).However, this negative effect does not always hold true (e.g. Park et al. 2011) and, in some cases, the results of different studies are not comparable as they may involve different species (Herrera et al. 2022) or use distinct metrics.To unify the criteria, we used a PCA in which aboveground biomass was the most important measure for defining meadow structure.Thus, we propose using this as a proxy measure in Z. noltei meadows to simplify further research and facilitate comparison of different studies.Determining the characteristics of the Z. noltei meadows is important as it could indicate the capacity to provide habitat and refuges from predation for both commercial and non-commercial benthic species (Bartholomew et al. 2000).
The infaunal biodiversity associated with intertidal soft bottom habitats has been addressed in several locations worldwide (e.g.Irlandi 1997;Boström and Bonsdorff 2000;Boström et al. 2006;González-Ortiz et al. 2014a), but recent literature concerning the distribution of the fauna on shellfish beds in relation to harvesting pressure is scarce (Vélez-Henao et al. 2021).The present study findings indicate greater biodiversity in areas with low fishing/harvesting pressure (i.e.established -NH-or less established -LH-seagrass meadows) than in areas with more continuous pressure (i.e.bare sediment).The lower biodiversity may be a consequence of the continual physical disturbance caused by clam harvesting, although the impact appears to be site dependent, as several studies have revealed that, in other locations, clam harvesting does not have a direct, significant impact on the associated benthic assemblage (Peterson et al. 1988;Boese 2002).Nonetheless, significant negative effects on particular areas, such as seagrass meadows, are known to occur worldwide, and we suggest that strong harvesting pressure leads to less complex habitats and, consequently, to a decrease in the associated infaunal biodiversity.Both laboratory and field studies have shown that dense and continuous meadows appear to enhance the biodiversity of infaunal species cohabiting with commercial species (González-Ortiz et al. 2014a;Crespo et al. 2017;Brun et al. 2021).
Although the sampling dates were chosen on the basis of previous knowledge about the reproduction and recruitment of clams (Martínez 2009;Matias et al. 2013;Moura et al. 2018), seasonal and/or monthly sampling for a longer period of 2-3 years would be desirable to confirm the patterns found in the present study.Our findings suggest that the continual clam harvesting pressure exerted in the most productive areas of shellfish beds throughout decades has had two clear impacts; (1) a notable increase in the abundance of the introduced clam R.
Vol:. ( 1234567890) philippinarum relative to that of the native clam R. decussatus and, (2) a reduction in establishment of seagrass meadows linked to a decrease in the associated biodiversity.The former may be due to shellfishers seeding the beds with R. philippinarum rather than with the native clam (Royo et al. 2002;Mantovani et al. 2006); this is reinforced by the fact that the introduced species has an extraordinary capacity to cope with environmental stressors such as elevated sediment temperatures (Domínguez et al. 2021).Seagrass recovery after harvesting disturbance is dependent on the presence of adjacent stable meadows (Qin et al. 2016) and its role as habitat and thermal buffer is crucial for both commercial and noncommercial species (e.g.Crespo et al. 2017;Román et al. 2022b).Management plans based on ecological knowledge of each bed and fruitful discussions with local stakeholders are therefore needed to protect seagrass meadows in productive areas of shellfish beds.Moreover, we propose integrating seagrass meadows in a long-term recovery programme for R. decussatus to enhance the abundance of the clam by periodic and controlled seeding to yield higher levels of recruitment and thus maintain a sustainable population structure in the near future.Otherwise, as the abundance of recruits was found to be very small in this study, the population of the native clam could collapse in the NW Iberian Peninsula in the next few decades, with negative ecological and socioeconomic consequences for SSFs and associated jobs.

Fig. 1
Fig. 1 Study area including the locations where the different environmental variables were measured and biological sampling was conducted.Sampling in areas of intensive, low and no harvesting (IH, LH and NH) pressure are shown.Sea surface temperatures in the main map and the 3 insets on the right-hand side correspond to an ECOSTRESS image obtained on 15-08-2021 09:23 UTC, approximately 2 h after high tide,

Fig. 2
Fig. 2 Mean (+ SD) abundance of R. decussatus populations for each size class, sampling date, shellfish bed and area.Raw abundances are shown above the error bars and the number of

Fig. 3
Fig. 3 Mean abundance (+ SD) of R. philippinarum populations for each size class, sampling date, shellfish bed and area.Raw abundances are shown above the error bars and the num-

Fig. 5
Fig. 5 Mean (+ SD) condition index of adults of R. decussatus and R. philippinarum for each sampling date, shellfish bed and area

Fig. 6
Fig. 6 Effects of the significant predictors for the best binomial model.A Observed (Obs) and model-predicted values (Pred) of R. decussatus abundance (number of individuals) at mature or immature stages in each shellfish bed.B Effects of the maximum temperature in emersion on R. decussatus maturity.C Effects on the interaction between mean salinity and mean chla on R. philippinarum maturity.Red dots and lines represent model-predicted values, and black dots represent observed values.Numbers refer to the number of observations within each stage of maturity for each value of maximum temperature in emersion in B and for each interaction value in C)

Fig. 7
Fig. 7 PCA analysis of Z. noltei metrics measured in the three shellfish beds.Percentage of variance explained and the linear coefficients of the most important variables are shown for PC1 and PC2.AGB: aboveground biomass, SL: shoot length

Fig. 8
Fig. 8 Mean (+ SD) A Shannon-Wiener diversity index (H') and B abundance of species coexisting with the commercial clams R. decussatus and R. philippinarum (i.e.associated species) in the different shellfish beds and areas . The conflict continues to occur because clam catches decrease Vol.: (0123456789)

Table 5
Mean (± SD) values of morphometric traits of Z. noltei populations in each shellfish bed and area.AGB: aboveground biomass, BGB: belowground biomass