Arcachon Bay (Fig. 1) is a 182 km² semi-enclosed lagoon located on the south-western Atlantic coast of France (Bay of Biscay) (Bouchet et al. 1997). It is a macro- to meso-tidal lagoon with a maximal tidal range of 4.9 m (Dodet et al. 2019) connected to the Ocean through a permanent 2-km wide channel. Its large intertidal mudflat represents ca. 110 km² and hosts eelgrass beds, primarily consisting of the dwarf eelgrass, Zostera noltei (Hornemann, 1832). Those meadows have been monitored since 1989 and, at that time, covered 70 km² of the lagoon (Auby and Labourg 1996). Then, the meadows strongly regressed during the 2000s resulting in a coverage of ca. 46 km² in 2008 (Plus et al. 2010), and ca. 39 km² in 2019 (unpublished observations). Salinity in the bay is influenced by freshwater inputs, mainly from the Eyre River in the southeast, and to a lesser extent by catchment runoffs, generating an East-West average salinity gradient from 20 to 35 (Auby et al. 1999).

Arcachon Bay is a hotspot for aquaculture activities. The main activity is the extensive farming of the Japanese oyster, Magallana gigas (Thunberg, 1793), combined with professional fishing activities of clams and/or fishing baits. Regarding clams, it mostly targets the Manila clam, Ruditapes philippinarum (A. Adams & Reeve, 1850), cockles (e.g. Cerastoderma edule (Linnaeus, 1758)) and mussels (e.g. Mytilus edulis Linnaeus, 1758). The last two have shown low stock levels in recent years.

The dataset used in this study was obtained during two field campaigns carried out to assess clam-standing stocks (Caill-Milly et al. 2021). These sampling campaigns took place in May-June 2018 and 2021. Samplings followed a stratified sampling strategy that targeted 23 intertidal strata including most tidal flats, except for some areas located in the south-western part of the lagoon (Fig. 1B). The delimitation of the strata was based on expert knowledge (scientists, fisheries and commercial fishers), each stratum representing an area as homogenous as possible in terms of hydrology, sediment characteristics, including particle size, and current patterns(Kermorvant et al. 2017, 2019). Respectively, 533 and 524 stations were sampled in 2018 and 2021 (Suppl. material 1). These stations were randomly distributed within strata using a Generalized Random Tessellation Stratified, GRTS, with a mean of 10 stations per square kilometre, proportionally to the strata’s area (Kermorvant et al. 2019). Samples were collected at high tide during spring tides with a Hamon grab (0.25 m²), which samples sediments down to 10 to 20 cm. Two replicates were made in most strata, except for C, D, I, J, K2, ML and N. Every replicate passed through three successive sieves of decreasing mesh size (20, 10 and 5 mm). All individuals of Arcuatula senhousia sampled by stations were measured and counted to estimate their abundance as a number of individuals by 0.25 m².

In addition, during the 2018 campaign, one sediment sample was collected at each station for grain size analyses. The sediment grain size was measured using a Malvern Mastersizer 2000 laser diffraction particle size analyser. These measures were not repeated during the 2021 field campaign.

Hydrodynamical model: MARS 3D

Environmental variables used in this study were extracted from a coastal environment modelling system: Model for Applications at Regionals Scales, MARS 3D (Lazure and Dumas 2008), designed by IFREMER. This model was adapted to Arcachon Bay by Plus et al. (2009). It takes into account tides, wind and river runoffs to simulate hydrodynamics taking into account the shape and bathymetry of the system. It permits running hydrodynamics simulations which provide salinity, temperature, current speed and erosion potential as well as other model-derived descriptors such as water residence times and age of water bodies (Plus et al. 2009). The erosion potential is a parameter including the sediment type and the current velocity. For this study, MARS 3D computed values for each variable at points of a grid with a resolution of 235 × 235 m, at the scale of the lagoon. For each point, the model computed a mean of a pool of values estimated every 30 minutes over one year. Thus, each point of the grid represents an annual mean of the environmental variables. The obtained points grid was then extrapolated in a raster format at the scale of the bay using the GIS software ArcGIS (ArcMap 10.8.2) by kriging (Ordinary) with the Geostatistical method: Geostatistical Wizard (e.g. Environmental Systems Research Institute, ESRI, Geostatistical Analyst).

Arcuatula senhousia habitat features in the lagoon

To define the main characteristics of the habitat where A. senhousia spread in Arcachon Bay, environmental variables that were proven to drive the distribution and life cycle of the species were chosen according to the literature (e.g., Reusch and Williams 1999; Bachelet et al. 2009; Yamamuro et al. 2010; Hosozawa et al. 2020; Watson et al. 2021). The selected parameters were (1) depth (in meters), (2) average salinity, (3) average water temperature (in °C), (4) average bottom current (in m.s^-1), (5) average erosion potential close to the bottom (day^-1) and (6) Zostera noltei meadow’s cover. The erosion potential represented the sheer force applying on the seabed. It was calculated by the MARS-3D model based on an imposed critical speed of 0.075 m.s^-1 for the boundary condition (Cognat 2019). Average salinity, water temperature and bottom currents were also derived from the output of this model. The meadow cover was estimated based on a mapping survey conducted in 2019 (Rigouin et al. 2022). It varies between zero (no seagrass) and one (fully covered sediment) and can be interpreted in the following way:

A two-way univariate permutation analysis of variance (PERMANOVA) without data transformation was used to test for differences in the presence and average density of A. senhousia between years (2018, 2021) and among strata. The design included 2 crossed factors, namely ‘year’ (2 modalities) and ‘strata’ (23 modalities). The analysis was based on Euclidean distances and used 9999 permutations, thus providing an analogous to traditional ANOVAs but where p-values are obtained by permutation (Anderson et al. 2008). Analyses were performed using the PRIMER software (Primer-E Ltd, 2012, version 6.1.15.). Pairwise tests were then conducted to put in evidence (1) strata showing significant temporal changes, and (2) between strata differences in mussel abundances.

One-way PERMANOVA on untransformed data was also used to test for differences in environmental parameters between stations where A. senhousia was present versus stations where no A. senhousia were found. The design then included a single factor ‘presence’ (two modalities) and the analysis was based on Euclidean distance as well as 9999 permutations.

Species distribution modelling

One of the objectives of the study was to estimate the probable distribution of Arcuatula senhousia in the lagoon. To estimate this potential distribution, in 2018 and 2021, a species distribution model (SDM) implemented in the R software environment using the package Biomod2 (4-1.2 version; Thuiller et al. 2009) was used (Fig. 2). SDMs are defined as models exploring the relationship between geographical occurrences of species and their corresponding environmental variables (Naimi and Araújo 2016). Presence-only records from 2018 and 2021 were used (Anderson and Raza 2010; Sillero et al. 2021).

Two scenarios were explored. The first scenario was based on the presence data of A. senhousia sampled in 2018. The obtained model provided a first estimate of the distribution of the species in the lagoon. This modelling method was then validated by comparing it to the presence data based on the 2021 sampling campaign.

Environmental parameters selection: Pairwise analysis

To characterise the distribution of Arcuatula senhousia in the bay, the environmental parameters were selected from the pool of available environmental data by a pairwise comparison using Pearson’s correlation coefficient. To deal with collinearity issues among variables in the model, a selection of parameters was performed when the correlation level exceeded 0.7 (Dormann et al. 2013). Thereby, the selected parameters were depth, erosion potential of surficial sediment and presence of Zostera noltei seagrass meadow (Suppl. materials 1, 2). These parameters are known to influence the phenology and distribution of A. senhousia within its native habitat range (Reusch and Williams 1999; Takenaka et al. 2018; Hosozawa et al. 2020).

Calculation algorithms

To create a map of the distribution of A. senhousia, four types of algorithms were chosen and applied: Generalised Linear Models, GLM; Gradient boosting machine, GBM; Random Forest, RF; and MAXimum of ENTropy model, MAXENT (Phillips et al. 2006). These algorithms were chosen from a pool of others after showing the best performance on the data and by following recommendations from the relevant literature (Reiss et al. 2011, 2015). More details about these modelling techniques can be found in Thuiller et al. (2009) or Reiss et al. (2015).

Model performance evaluation

To assess the models’ performances, three different comparative methods were used (1) the Receiver Operating Characteristic (ROC), (2) the True Skill Statistic (TSS), and (3) Cohen’s KAPPA coefficient (k). Model discrimination was computed using the area under the receiver operator characteristic curve (Area Under the Curve, AUC). This provides information on the model’s ability to make a distinction between occupied and unoccupied sites. In this study, models with an AUC coefficient over 0.7 (for the three comparative methods) were considered to display good performances. Then, those with coefficients describing the true positive (sensitivity) and true negative (specificity) above 0.8 were also considered correct.

Implementation characteristics

Pseudo-absences were randomly allocated within the study area. The number of pseudo-absences was chosen as equal to the number of presence stations (Barbet-Massin et al. 2012). We ran five iterations and 30 cross-validations for each modelling algorithm, which thus produced a total of 600 models. For each of them, the data was split into 70% for calibration and then 30% for validation. Each model was then tested to evaluate the predictive model and to select the most robust one to build an ensemble model. An ensemble model is a way to synthesise the results of a pool of chosen algorithms to combine forecasts. It was realised with the BIOENSEMBLE option of the Biomod2 package. A threshold of TSS > 0.8 was used to estimate which models would be selected in the BIOENSEMBLE model (Allouche et al. 2006).

Finally, the results were synthesised as a raster, which showed a range of habitat suitability, i.e. the property of the habitat to allow the survival and persistence of the species (Araújo and New 2007). This index varies from 1 when the environmental variables seem suitable for the presence of the species, to 0 when the variables are not suitable for the presence of the species.

A total of 150 and 2,050 individuals of Arcuatula senhousia were sampled during the surveys carried out in 2018 and 2021, respectively. These individuals were characterised by shell length between 10 to 35 mm (23.3 ± 5.5, average ± sd), and between 5 to 35 mm (16.5 ± 6.4, average ± sd) in 2018 and 2021, respectively.

The presence of A. senhousia significantly increased (main test of the two-way PERMANOVA, p(perm) < 0.05, Pseudo-F = 137.5) between 2018 and 2021. A. senhousia was found at 62 sampled stations out of 533 in 2018 and 188 sampled stations out of 524 in 2021. Strata colonised in 2018 remained colonised in 2021, with an increased number of species occurrence in most cases (Fig. 3A, B).

In 2021, A. senhousia was significantly more present (pairwise tests of the two-way PERMANOVA, p(perm) < 0.05) in the northwest strata sheltered from dominant wind waves by an island (Ile aux Oiseaux) and in the southeast strata close to the main freshwater input (i.e. the mouth of Eyre River). In those areas, the species occurred in more than half of the sampled stations. In 2021, A. senhousia was also observed in strata where it had not been monitored in 2018.

Furthermore, densities of A. senhousia significantly increased between 2018 and 2021 (main test of the two-way PERMANOVA, p(perm) < 0.05, Pseudo-F = 65.2), in the same northwest and southeast axis. However, the highest densities were mostly found in the areas shielded by the island. The highest densities observed were 92 and 270 ind.m^-2 in 2018 and 2021, respectively.

The sediment type in the samples was homogeneous. The mean D50 measured was of 50 µm +/- 7.3 µm (mean +/- SD). 81% of the sampled points in 2018 were classified as “sandy clay”, 17% as “muddy sand”, and 2% as “sand”.

No significant differences in temperature and salinity were found when comparing the stations where A. senhousia was retrieved to those where it was absent (one-way PERMANOVA, p(perm) > 0.05, Pseudo-F respectively 9.8 and 3.2 for salinity and temperature). Fig. 4A, B indeed showed little variation for those parameters. On another hand, there were significant differences in current speed, erosion potential, depth and meadows’ coverage (one-way PERMANOVA, p(perm) < 0.05, Pseudo-F respectively 52.3, 26.7, 21.1 and 15.4). A. senhousia occurred more frequently within the lowest values of current speed (< 0.1 m.s^-1, Fig. 4C) and erosion potential (< 35 day^-1, Fig. 4D). Regarding intertidal elevation, the species occurred throughout the mid-littoral zone (ca. -1.5 to -3 m in Arcachon Bay) but its frequency declined toward the high intertidal (Fig. 4E). Moreover, the species was significantly less frequent in the lowest meadows coverage (< 0.3, Fig. 4F).

Maximum densities were found at current and erosion potentials ranging respectively from 0.05 to 0.15 m.s^-1 and 0 to 75 day^-1 (Fig. 4C, D). However, statistical analysis did not support any significant differences in densities between the parameter ranges for all environmental parameters (one-way PERMANOVA, p(perm) > 0.05, Pseudo-F = 13,5).

First estimation of the probability of presence: 2018

Ensemble model based on the 2018 presence

The overall validation scores were high (Table 1). The AUC was above 0.8 for ROC and TSS, which indicated very good performance of the ensemble model. It was above 0.6 for Cohen’s coefficient (KAPPA), which represented a good performance for the model. The sensitivity was above 95 for both ROC and TSS, but was above 75 for the KAPPA, which showed, according to this metric, that the model has a good chance of predicting true positive values. The specificity was near 90, which qualified the model as a good predictor of the true negative. The overall values showed a model with a good estimation of both true positives and true negatives.

The suitability index map (Fig. 5) indicated that the highest index was between 0.8 and 0.9. These values indicated environmental parameters theoretically favourable for the presence of A. senhousia. The areas characterised by this index range mostly followed the same northwest/southeast axis previously described. Furthermore, the remaining intertidal area was identified with indices between 0.5 and 0.8. Although lower, they qualified areas where the species could still be present. Conversely, areas around channels were always characterised by an index lower than 0.2, indicating lower chances of finding the species.

Validation of the model estimates

Comparing the species occurrences recorded in 2021 (Fig. 5) with the ensemble model based on 2018, the modelling results gave a fairly accurate estimation of the potential distribution of A. senhousia: positive occurrences in 2021 were mostly located in areas with a coefficient above 0.6. However, some discrepancies between the results of the model and the presence of A. senhousia were observed at a few stations located near the mouth of the Eyre River, where the values of indices were low (between 0.1 and 0.3).

In the previous observation made based on the dataset collected in 2021 (part 1, p9), the north of the bay was characterised by the highest occurrences and densities of A. senhousia. Based on the model output (Fig. 5), this area was indeed showing a suitability index above 0.8. Furthermore, in 2021, some of the strata were described as newly colonised and the ensemble model based on 2018 data described them as suitable for the establishment of A. senhousia. The model also pointed out the area in the south of the lagoon, which was not sampled in 2021 and 2018, as a suitable habitat.

Probable distribution of the species in the bay: based on 2021 presence

The overall validation coefficient was above 0.8 according to all three methods (Table 2). The sensitivity and the specificity were also above 90. This ensemble model had an excellent performance with a good appreciation of the true positives and the true negatives.

Looking at the map of the suitability index (Fig. 6), the greater part of the intertidal area was characterised by an index above 0.8. Around 90% of the intertidal area was characterised by a homogeneous probability of being colonised by A. senhousia. As predicted in Fig. 5, the south of the bay also showed a suitability index above 0.8. In addition, the southwest of the Ile aux Oiseaux and the strata at the middle of the lagoon were also characterised by these high indices. The area west of the river mouth is still showing values of the indices slightly lower than the rest of the intertidal flats, but still with values above 0.7. This area could still be suitable for the presence of A. senhousia. The channels were the only area characterised by values below 0.2.

Characterising the environmental context associated with the spread of a NIS through time at the scale of an ecosystem requires the use of a dataset able to capture a sufficient amount of the inherent environmental variability within which populations of the species of interest could sustain and extend (Sillero et al. 2021). Moreover, the spatial resolution at which these different types of data are acquired (or aggregated) must be high enough and area-to-point spatial misalignment should be avoided (Mourguiart et al. 2024) to allow for the establishment of meaningful relationships between the presence/absence (or density) of the species of interest and several key environmental abiotic and biotic variables. In this study, we likely fulfilled such requirements to study the spatiotemporal variations and the ecological niche of Arcuatula senhousia in Arcachon Bay by using field survey data describing the Asian mussel distribution combined with environmental data derived from (1) bathymetric survey, (2) airborne habitat mapping, and (3) numerical hydrodynamic model outputs aggregated at the same high spatial resolution (235 × 235 m).

Differences in average salinity and temperature did not explain the spatial distribution of the species in Arcachon Bay. This can be due to the relatively small variation of these parameters in the lagoon (temperate ecosystem). Other authors have found A. senhousia in a wide range of salinity and temperature conditions, ranging from 5 to 30 (e.g. Takenaka et al. 2016; Hosozawa et al. 2020). Knowing the global distribution of A. senhousia, it appears that salinity and temperature measured in a temperate ecosystem would not necessarily influence the distribution of the species within the lagoon. However, that does not mean that those parameters do not affect the species. Indeed, Liang et al. (2009) showed that temperatures below 15 °C and salinities below 20 can impact embryo and spat survival. Consequently, the variability of these environmental factors, instead of their average value should be considered.

Conversely, PERMANOVA results showed that there were significant differences in (1) depth/elevation, (2) current speed, (3) erosion potential and (4) seagrass cover between locations where the Asian date mussel occurred compared to locations where it was absent. Firstly, several authors have shown that A. senhousia was mainly distributed in the intertidal flat (Morton 1974; Creese et al. 1997; Crooks 1998; Como et al. 2018). Intertidal habitats extend over 110 km² in Arcachon Bay. Consequently, it offers favourable conditions for the presence of this predominantly intertidal species. In this dataset, stations located on the highest part of the intertidal flat (higher than 3.5–4 m) were less colonised. This corresponds to the area above the average spring tide watermark. Although this study focuses mainly on intertidal flats, other monitoring surveys in this bay failed to show the presence of A. senhousia in subtidal areas (De Montaudouin et al. 2018). However, in both native and introduced ecosystems, this species was observed in lagoonal or estuarine subtidal areas (as in Japan or a French Mediterranean lagoon; Hoenselaar and Hoenselaar 1989; Hosozawa et al. 2020). Secondly, the probability of the presence of A. senhousia decreases considerably with the increase in current speed and erosion potential (Fig. 4C, D), which potentially regulates the mussel population. This suggests that strong currents negatively affect the species’ ability to stay anchored within the sediment matrix, which may impede its presence in subtidal channels. This relationship is consistent with field observations made after periods of storms, showing high numbers of A. senhousia individuals washed ashore (Massé C., Bernard G., personal observations). For other species of bivalves, it was observed that storms (Nehls and Thiel 1993), waves and currents (Hunt and Scheibling 2001) could dislodge individuals. The observed pattern of high probability of presence and densities also coincided with low current and erosion potential values. Furthermore, although the hydrodynamic model used in this study did not account for the effect of swell and wind-wave events, it can reasonably hypothesise that the hotspots of the presence of A. senhousia we identified north-eastern from the Ile aux Oiseaux island could be linked with relative protection from erosive events connected with the dominant wind-wave regime. Besides their local negative effects on A. senhousia’s presence, strong currents could also be a way of propagation for the species. A. senhousia has a pelagic larval stage that ranges between 4 and 50 days and the larvae can disperse up to 116 kilometres away from spawning grounds (Sgro et al. 2002; Shanks 2009). In addition, the species has been observed in association with macroalgae, such as Gracilaria vermiculophylla in the lagoon (Bernard et al. 2020). This association could promote a secondary dispersal mechanism by the transport of adults fixed to algae at a larger scale, as rafting effects are known to influence the dispersion of individuals of a species and their hitchhikers (Kinlan and Gaines 2003).

Finally, according to our results, A. senhousia was less present in bare sediment compared to sediment covered by seagrass (Z. noltei) meadows. Two non-exclusive hypotheses can support this finding. First, Thomsen (2007) showed that substratum provision was one of the key factors in invasion processes. Accordingly, Reusch and Williams (1998) demonstrated that low-foliage meadows promote the presence of A. senhousia. Indeed, Z. noltei meadows may promote A. senhousia occurrence in the bay, by their roots and rhizomes network, which seem to provide an anchoring point for species individuals and a substratum for settlement. Secondly, seagrass meadows and their canopy are known to substantially attenuate the effects of waves and to reduce current velocity in their vicinity (Fonseca and Fisher 1986; Ganthy et al. 2015). Furthermore, roots and rhizomes networks of Z. noltei meadows also stabilise the sediment matrix, limiting erosive fluxes (Ganthy et al. 2011) as well as biologically induced sediment mixing by burrowing infauna (Bernard et al. 2014). Z. noltei meadows then appear to provide a physically stable shelter to individuals of A. senhousia, lowering the levels of the two main variables previously found to hinder their presence in the lagoon, namely current velocity and associated erosion potential.

Those observations are raising concerns about the probable impact of A. senhousia on the seagrass meadows in the lagoon. Indeed, following seagrass meadows regression in Arcachon Bay (Plus et al. 2010), particular attention should be paid to the interactions between those two species. Reusch and Williams (1999) already showed that high densities of A. senhousia impede the Zostera marina’s leaves and rhizomes growth. Furthermore, A. senhousia forms small patches that generate bulges at the sediment surface. Those bulges, by generating micro topographical changes, could locally increase currents and hence contribute to meadow regression (Cognat et al. 2018). The identified elements emphasise the effects on the mussel and the potential secondary impacts that could cause habitat regression because of the species’ presence.

The comparison of the data from 2018 and 2021 with previous studies (i.e, 2002, (Bachelet et al. 2009); Do 2012 and Dubois 2012; Fig. 7) indicates that the distribution of Arcuatula senhousia has changed drastically over time and space. Indeed, Arcuatula senhousia was first observed in 2002, in the northwest areas of the lagoon (Fig. 7A; Bachelet et al. 2009) in the same sheltered place described before (e.g. north of the island, Ile aux Oiseaux). Then, the Asian date mussel was distributed in various sites throughout the bay (Fig. 7B–D) and was not limited to a particular area. The spread of colonised areas followed a northwest/southeast axis in the lagoon. Furthermore, a consistency in the presence of the species was observed throughout the surveyed period. In most cases, areas colonised in 2002 were still colonised in 2009, 2018 and 2021. The north of the bay seems to be the hotspot of the species’ presence in the bay. This suggests that A. senhousia is strongly established in this area and may be a historical settling area for the species in the lagoon. Although two temporally close sampling surveys are not sufficient to validate the long-term dynamics of this species, the clear increase in abundances between 2018 and 2021 appears to reflect an invasion process. More surveys need to be carried out to confirm this hypothesis.

The species seemed to have spread in two ways. (1) Gradual expansion of the species’ distribution around core areas already colonised. Indeed, newly colonised stations in 2021 were located in the vicinity of the stations where A. senhousia was present in 2018 (Fig. 7). To a lesser degree, this tendency is also observed when comparing 2009/2010 and 2018. Indeed, the settlement and survival rate of many bivalves’ juvenile is most often higher when they are in the vicinity of adults (Wilcox and Jeffs 2017). In this case, adult specimens’ presence could provide an anchor and protection for juveniles’ settlement and therefore increase their abundance. (2) Colonisation of areas far from previously colonised ones by dispersion during the larval stage or by rafting effect. Indeed, between 2018 and 2021, the species appeared in some strata where it had not been seen before. This observation goes in pair with the hypothesis made on the role of current in the distribution of larval and adult individuals of A. senhousia.

Until the observations made in 2018 and 2021, the species was supposed to be in a lag time, as Crooks (2005) described, during which its population did not show any substantial development. However, the recent increase in presence and densities questions a possible species’ dynamic change in the lagoon. According to Boudouresque and Verlaque (2012), two invasion kinetics models could occur. (1) The “boom and bust model” is characterised by a rapid increase in abundance, followed by a collapse with a low but stable number of individuals. This kinetic type has been observed for the same species in other introduced ecosystems, such as Perth, Australia (McDonald and Wells 2010), and Sardinia, Italy (Como et al. 2018). (2) The “natural fluctuation model” is conversely characterised by the population reaching a plateau and is less impacted by fluctuations. To our knowledge, this type of kinetics was not observed in studies carried out on A. senhousia.

However, in 2021, the species’ densities did not exceed 200 ind.m^-2. Furthermore, following another field observation in Arcachon Bay (Benard G. & Massé C. comm. pers.), density in “reefs” was estimated at around 2,000 ind.m^-2. This number is lower than in other areas where A. senhousia was introduced. Indeed, in the Adriatic Sea (Mediterranean Sea) or Mission Bay (USA), densities of more than 10,000 ind.m^-2 were observed (Mistri et al. 2004; Castorani and Hovel 2015). Arcachon Bay appears to be still in an early stage of colonisation. The low densities present in the lagoon could be a clue that the species may still be in its lag time, as Crooks (2005) described. Furthermore, the biomass of A. senhousia does not exceed 5% of the biomass of the congeneric species Ruditapes philippinarum collected during the same sampling campaign.

Limits of the study

Differences between the model results based on 2018 and 2021 highlight the uncertainty of the equilibrium in Arcachon Bay, which, as Barbet-Massin et al. (2018) explained, could lead to an underestimation by the model of the probable area invaded by A. senhousia. Future high-resolution surveys will allow us to confirm the prediction proposed by the 2021 model. Furthermore, including additional sampling points south of the Bay would enhance model validation and improve our understanding of species distribution.

Moreover, the modelling process does not consider additional factors that could influence the presence of the species, such as anthropic activities (Reiss et al. 2015) or biotic interactions such as predations or competitions (Srivastava et al. 2019; Stephenson et al. 2022). It was seen in both native and introduced ecosystems that predation was an important factor in explaining the distribution of A. senhousia. The potential predators of A. senhousia are ducks (Yamamuro et al. 2010), crabs (Cabiddu et al. 2023), shorebirds and fishes (Crooks 2002). There are limited observed cases of predation on A. senhousia in Arcachon Bay, but fragments of A. senhousia shell were found in shorebird faeces (Latry comm. pers.).

Prediction of the species distribution

The two prediction maps gave an idea of the probability of the presence of the species in Arcachon Bay. Results show accordance between areas characterised as suitable by the model based on 2018 and in situ observations of A. senhousia in 2021. Then, the change in the predictions between the two maps confirms that the species distribution has increased and probably could still increase in the future. Furthermore, the model based on 2021 permits to point out areas to monitor in the future where the presence of the species is probable. Both ensemble models pointed out that the south of the lagoon seems to be suitable for the presence of the species. Those areas were not sampled during the campaigns of 2021 and 2018, but field observations confirmed its presence in this area in 2023 (pers. observation). According to the introduction and invasion theories described before, and the recent observations, a scenario of homogenisation of the presence and settlement of A. senhousia in all areas of the bay is possible. Further monitoring in this area of interest is necessary to follow the spread of this species.

The observation made with the modelling results permitted to characterise the lagoon as a typical habitat where the species could be present. This information is of primary interest for management purposes, especially in the case of exotic species. One aim of this work was to test if this modelling method could be used as an indicator to point out areas to monitor. The results obtained for A. senhousia in the lagoon were conclusive and helped to provide additional information on the distribution of the species. Even if we pointed out limits to this method, SDMs are a useful tool to help in the detection of areas to survey. Indeed, non-indigenous species detection is now one of the priorities for managing ecosystems since corridors of the introduction of species have been considerably increasing since 1950 (Hulme 2009). This is especially the case of introduced species such as Arcuatula senhousia, which seems to be introduced in a large variety of ecosystems. Indeed, introductions of this species in England (Watson et al. 2021) or Portugal (Encarnação et al. 2024), as well as a northward expansion of its presence along the French Atlantic coast (Massé et al. 2022) have been increasingly recorded. This underlines that the present study could contribute to a comprehensive knowledge of the main drivers controlling invasion dynamics at a larger scale and therefore to identify future potential areas to be colonised. The next step will be to try this method at a larger scale, maybe coupled with other modelling methods [e.g. COUE (Liu et al. 2020)] and find suitable ecosystems to monitor on the Atlantic coast of France.

This work was supported by the French Biodiversity Agency, Office Français de la Biodiversité (OFB), under Grant Agreement OFB 22-0143/UB.

The data from the Manila clam campaigns (2018 and 2021) used for this work were collected thanks to funding from DLAL-FEAMP Affaires Maritimes et la Pêche du Bassin d’Arcachon-Val de l’Eyre, the Nouvelle-Aquitaine Region, CDPMEM33 (Comité Départemental des Pêches Maritimes et des Élevages Marins de Gironde) and Ifremer. Data on the occurrence of Arcuatula senhousia before 2018 were acquired thanks to the financial support of the Syndicat Intercommunal du Bassin d’Arcachon (SIBA). The latter sampling campaigns were performed with the help of the staff and ships Planula III and IV fitted out by University Bordeaux 1 and CNRS/INSU (Flotte Océanographique Française), respectively.

The authors have no relevant financial or non-financial interests to disclose.

Material preparation and data collection were performed by Florence Sanchez, Nathalie Caill-Milly, Muriel Lissardy and Guillaume Bernard. Environmental data were provided by Florian Ganthy. Analysis were performed by Salomé Coignard. Interpretation of data was made by Salomé Coignard, Marie Fouet, Hugues Blanchet, Guillaume Bernard and Cécile Massé. The first draft of the manuscript was written by Salomé Coignard, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.