The study area includes the coastline of Mediterranean Sea, a semi-enclosed basin connected to the Atlantic Ocean by the Strait of Gibraltar. The average depth is around 1500 meters, ranging from 0 to 5270 m, with the deepest point found in the eastern Mediterranean Sea, in the Calypso Deep, in the Ionian Sea. The Mediterranean Sea borders Europe, Africa, and Asia with a Mediterranean climate with mild to cool, rainy winters and warm to hot, dry summers. A wide variety of climatic and anthropogenic pressures affect the area (Durrieu de Madron et al., 2011). For example, thermal anomalies due to rising water temperature caused mass mortality events among marine species (Crisci et al., 2011). Human activities exert significant stress on Mediterranean marine ecosystems, including (i) coastal development, which alter natural habitats, (ii) increased input of nutrients due to farming, commercial and industrial activities, (iii) non-native species introduction through shipping, aquaculture and the Suez Canal opening, and (iv) overfishing, which depletes fish stocks and alters community structures (Katsanevakis et al., 2014; Korpinen et al., 2021).

Occurrences of Atlantic blue crab were obtained by extracting data from the Global Biodiversity Information Facility (GBIF.org, 2024). The whole dataset of blue crab in the GBIF was filtered to include only data with georeferenced information obtained from human observations of living organisms recorded between 2021 and 2023 along the Mediterranean coasts. These data were contributed by citizen science volunteers who recorded sightings in the citizen science platform iNaturalist (https://www.inaturalist.org/). Data from citizen science plays a critical role in monitoring biodiversity, especially at large-scales (Cardoso et al., 2017; Lucy et al., 2016) and in monitoring biological invasions (Encarnação et al., 2021; Pocock et al., 2024). To ensure data accuracy, observations were cross-checked to avoid duplicate species occurrences. While comprehensive datasets on the occurrence of Atlantic blue crab in Mediterranean Sea are available covering its historical detection in the early 1900s (Castriota et al., 2024; Mancinelli et al., 2021), we chose to extract data directly from the GBIF, because it provided updated occurrence data at the scale of the entire Mediterranean Sea for the recent period 2021-2023. This choice of data aligned the blue crab occurrence data with environmental and human activities data over the same period. Consequently, recent blue crab occurrence data inevitably lack information on the status of older occurrences (which in Mediterranean started from 1937 according to Mancinelli et al., 2021, as well as possible crab locations not captured in our dataset. However, the overlap between our data and the more comprehensive datasets provided by Castriota et al. (2024) and Mancinelli et al. (2021) showed that our data cover a significant proportion of the areas where blue crabs have been historically recorded ( Supplementary Figure 1 ).

NetCDF environmental data of water temperatures (°C), salinity (P.S.U.), nitrate (mg·m^-3), orthophosphate (mmol·m^-3) and chlorophyll-a (mg·m^-3) for the whole Mediterranean Sea was retrieved from the database of CMEMS (Copernicus Marine Environment Monitoring Service, available at https://data.marine.copernicus.eu/products, Supplementary Figures 1A–F ). Monthly water temperatures and salinity data were acquired from 2022 to 2023 with resolution of 0.083°x0.083° (E.U. Copernicus Marine Service Information (CMEMS), 2024b). Monthly data of nitrate, orthophosphate and chlorophyll-a were acquired from 2022 to 2023 with resolution of 0.25°x0.25° (E.U. Copernicus Marine Service Information (CMEMS), 2024a). Summer and winter water temperature means annual salinity, phosphate, chlorophyll-a and nitrate means ( Supplementary Figures 2A–F , respectively) were calculated for the surface water using QGIS software (QGIS.org, 2024).

Data on water depth ( Supplementary Figure 2G ) and the density of commercial ships (hereafter, ship density) in the Mediterranean Sea were derived from the European Marine Observation and Data Network (EMODnet, 2014). Ship density data for 2022 was extracted in a GeoTIFF format, representing the annual number of ship route crossings in 1x1km cells of a grid covering all Mediterranean Sea waters ( Supplementary Figure 2H ). Data on coastal wetlands ( Supplementary Figure 2I ) was derived in a GeoTIFF format with a resolution of 0.008°x 0.008°from Lehner and Döll, (2004).

Human population data were collected from the Socioeconomic Data and Applications Center (CIESIN, C. University, 2018) in a GeoTIFF format and represent the population size with resolution of 0.0083°x0.0083° ( Supplementary Figure 2J ). The most recently available data on human population size was for 2020. They are included in the analysis to investigate both the potential linkage to cultural eutrophication and whether observations in blue crab occurrence data might be biased by human abundance (i.e., more people more observations).

For each 50 km of coastline, the total number of observed occurrences of blue crabs was calculated. We assumed that higher numbers of observed occurrences along 50 km of coast are expected in areas with higher abundance of blue crab. This assumption is necessary since density or abundance data are not available in GBIF for C. sapidus, and studies with quantitative data on this species are not available for the entire Mediterranean Sea. A good approximation between the number of observations in GBIF and the abundance of a species in the world has already been used and verified for bird species (Callaghan et al., 2023). However, to account for potential bias due to the higher numbers of crab observations occurring in higher densities of humans, we also calculated the proportion of blue crab occurrences relative to local human population size (hereafter “proportional occurrence”).

Water physical variables (i.e, water temperature, salinity), water quality variables (i.e., nitrate, orthophosphate, chlorophyll-a), water depth and human activities (i.e, shipping density and population count) data were paired with each coastal stretch of 50 km where blue crabs occurred using QGIS software (QGIS.org, 2024). To do this, for each 50 km coastal stretch, the average values of environmental and shipping density data were calculated up to 5 km offshore of the coastline. Similarly, for each 50 km coastal stretch, the total population and the extension of coastal wetlands present in 25 km onshore (on the land side) were calculated. At the end of this process, for each coastal stretch of 50 km of length, the total number of blue crab occurrences, the averages of summer and winter water temperatures (°C), annual salinity (P.S.U.), annual nitrate (mg·m^-3), annual orthophosphate (mmol·m^-3) and annual chlorophyll-a (mg·m^-3), water depth, ship density and total population were obtained ( Figure 1 ; Table 1 ).

All analyses were performed in RStudio software (R Studio Team, 2024). Linear mixed-effect (LME) models were used to test the association of geomorphic, water physical variables, water quality variables and human activities with blue crab occurrences and proportional occurrence. To avoid strongly correlated variables in the models (r_s>0.7), the correlations between variables were investigated through Spearman correlation using corrplot R package (Wei and Simko, 2017). The country of the observation was included as a random factor to account for nested observations. The explanatory variables were standardized by centering and scaling the values (Bates et al., 2015). Response variables (i.e., blue crab occurrences and proportional occurrence) were log-transformed with log10(x+1). Models were validated by checking residual patterns in the “DHARMa” R package (Zuur and Ieno, 2016; Hartig, 2021). The lme model was fitted using the R package ‘lme4’ (Bates et al., 2015). The best model was selected based on the Akaike Information Criterion (AIC; Akaike, 1974) corrected for small sample sizes (AICc; Hurvich and Tsai, 1993) using AICcmodavg R package (Mazerolle, 2019). The selection of the best model was based on Akaike weights (models with large Akaike weights have strong support) and low AICc values (Snipes and Taylor, 2014). Explanatory variables with p-values<0.157 were retained in the best models as they still hold explanatory power and can describe true relationships (Sutherland et al., 2023). The Variance Inflation Factor (VIF<4) of the explanatory variables used in the selected model allowed assessment of the presence of collinear variables (Zuur et al., 2009), using “car” R package (Fox & Weisberg, 2020).

To estimate the variance explained by each of the fixed and random effects of the best models selected, the marginal and conditional R2 values were calculated for each linear mixed model (Stoffel et al., 2021). The marginal R2 gives an estimate of the variance explained by each fixed effect relative to the total variance in the response, whereas the conditional R2 gives an estimate of the variance explained by fixed effects and random effects together, which better reflects the heterogeneity of the variables. The 95% confidence intervals were estimated for the marginal and conditional R2 using 1000 parametric bootstrap iterations (Stoffel et al., 2021). Marginal and conditional R2 were calculated with “partR2” R package (Stoffel et al., 2021).

Among geomorphic variables, the water depth of the coastal strip extending 5 km from the coastline where the blue crab was reported ranged from 4.9 to 202.3 m, with a mean depth of 45.3 ± 41.6 m. Onshore, coastal wetlands within 25 km of the coastline extended from 0 to 708 km², with a mean area of 31.2 ± 127.2 km² ( Table 1 ). Among water physical variables, summer water temperatures ranged from 20.8 to 28.2°C, with an average value of 25.7 ± 1.2°C. Winter water temperatures ranged from 9.1 to 19.3°C, with an average value of 15.0 ± 2.1°C. Water salinity ranged from 30.5 to 39.4 PSU, with an average value of 37.4 ± 1.8 PSU ( Table 1 ). Water quality variables showed a wide range of values with highest values in the northern Adriatic Sea ( Table 1 ; Supplementary Figure 2 ). The shipping density (i.e., annual shipping routes per km²) ranged from 0 to 2719.18 annual shipping routes per km², with an average of 62.1 annual shipping routes per km². Human population size around 25 km onshore of the 50 km stretch ranged from 2695 to 5,302,038 people with an average value of 530,509 people ( Table 1 ). Means and medians were similar for 4 variables: summer and winter water temperatures, water salinity and chlorophyll-a. The medians for the other 8, including crab occurrence variables, were less than means, indicating skewed distributions to the right.

From 2021 to 2023, the blue crab occurrences were reported along the Mediterranean coast from 13 countries, in a total of 102 stretches of 50 km each where crab occurrences ranged from 1 to 42, with highest values in northern Italy, Albania, Greece and Spain (>26 observations·50 km^-1; Figure 2A ). Proportional blue crab occurrences exhibited comparable patterns, with a more pronounced invasion in the eastern Adriatic Sea and Greece ( Figure 2B ). The number of blue crab occurrences is higher along the coasts of Europe than along the coasts of North Africa and Eurasia ( Figure 2 ).

There was correlation among some of the independent variables. Water quality variables (i.e., Chla, NO3 and PO4) were negatively correlated between WinterTemp, Salinity and Depth, whereas they were positively correlated with PeopleSum and CoastWet ( Supplementary Figure 2 ). CoastWet resulted negatively related to Salinity and Depth (r_s=-0.28 and r_s =-0.30 P-values<0.05, respectively; Supplementary Figure 3 ). Shipping density showed a positive correlation with PeopleSum (r_s=0.19, P-values<0.05; Supplementary Figure 3 ).

The best LME model retained from the AIC model selection process (Mod3) indicated that blue crab occurrence (response variable) increased significantly with decreasing depth, and increasing nitrate concentrations ( Table 2A , Figure 3A ; Supplementary Table 1A ). The partitioning of R² showed that NO3 was the variable with the highest value of conditional and marginal R² explanatory ability (conditional R²: 0.32, IC: 0.19–0.46, marginal R²: 0.21, IC: 0.07-0.36) followed by water depth (conditional R²: 0.16, IC: 0.00–0.20, marginal R²: 0.03, IC: 0.00–0.22; Figure 4A ).

Models with the combination of nitrate, chlorophyll and water depth were less plausible (Mod4 and Mod5; ΔAICc=5.2 and ΔAICc=1.07, respectively, Table 2A , followed by models with only environmental variables (Mod2; ΔAICc>18.5; Table 2A ), with all variables (Mod 1, ΔAICc=29.90; Table 2A ) and with only physical parameters (Mod 6; ΔAICc=29.0; Table 2A ). The model including only human activities (Mod7) was the least plausible (ΔAICc=32.8; Table 2A ).

AIC model selection analysis using proportional occurrence of blue crabs as response variable produced similar results, retaining the model with depth and nitrate concentration as the best model (ModS3; Table 2B ). According to ModS3, the proportional occurrence of blue crabs showed a positive relationship with both variables, although depth resulted in lower explanatory ability than nitrates ( Figures 3B , 4B ; Supplementary Table 1B ). The partitioning of R² showed that NO3 was the variable with the highest value of conditional and marginal R² in explaining the proportional occurrence of blue crab (conditional R²: 0.20, IC: 0.08–0.37, marginal R²: 0.16, IC: 0.06–0.31) followed by water depth (conditional R²: 0.05, IC: 0.00–0.24, marginal R²: 0.01, IC: 0.00–0.17; Figure 4B ).

Blue crab occurrences (also corrected by human population size) from 2021 to 2023 were predominantly concentrated along the European coastline with few reported occurrences along African shores, mainly located in Marocco, Algeria and Tunisia. The highest blue crab occurrences in Spain, Greece and Italy coincided with areas with large rivers that flow into the sea or create estuarine lagoons. These large river outlets and coastal wetlands provide brackish environments which appear to be relevant in the spatial variation of blue crabs. These habitats provide nutrient loads which promote high primary production and more food availability for blue crab, as discussed below. Also, shallow waters and salinity gradients are suitable for the blue crab reproduction (Hines, 2007). First occurrences and large spreading of blue crab populations were in fact reported in such Mediterranean brackish environments such as the Ebro Delta in Spain (Castejón and Guerao, 2013; Clavero et al., 2022), Venice Lagoon (Nehring, 2011) and Po River Delta in Italy (Manfrin et al., 2015), Gulf of Thessaloniki in Greece (Kampouris et al., 2020; Nehring, 2011), estuary of Oued Z’hor in Algeria (Benabdi et al., 2019) and Bizerte Lagoon in Tunisia (Shaiek et al., 2021).

Compared to extensive, historical occurrence data provided by Mancinelli et al. (2021) and Castriota et al. (2024), we found a lack of occurrences from 2021 to 2023 in regions where blue crabs have been previously documented but where they are still likely. For example, blue crab occurrences were reported in Egypt (Abdel Razek et al., 2016; Mehanna et al., 2019), Gulf of Lions in France (Labrune et al., 2019), western part of Sardinia in Italy (Culurgioni et al., 2020) and Malta (Vella et al., 2023). This suggests a lack of citizen reporting rather than absence of blue crab populations in these areas, and represents a limitation of this and similar data. Nevertheless, this discrepancy was limited to few areas, and our analysis covered almost all regions where blue crabs have been reported previously. Consequently, we inferred that while all specific regions may not be represented, our results were generally robust for interrelationships among variables.

Although the blue crab life cycle (e.g., larval recruitment, juvenile growth, reaching of maturity and reproduction) is known to be influenced by seasonal temperature variations, as demonstrated by studies on the U.S. Atlantic coast (Hines, 2007), neither winter nor summer water temperature was associated with higher blue crab occurrences in this study, suggesting that other factors primarily influence large-scale crab distribution. On the other hand, the blue crab exhibits a large thermal adaptability ranging between 0-40°C with an optimum value at 24°C along a Thermal Performance Curve (Marchessaux et al., 2022), and dormancy induced below 9°C (Brylawski and Miller, 2006). During the summer season, optimal water temperatures for blue crab metabolism were reached throughout the Mediterranean, but maximum tolerance levels were not exceeded, at least at the Mediterranean scale in this study. This suggests that warm water temperature alone is unlikely to limit or prevent population growth in this region. Similarly, also lower temperature reached in winter appear to not be an obstacle for blue crab populations growth as global warming will likely reduce severe winter in the region (Shaltout and Omstedt, 2014). Furthermore, increasing water temperature under global change, rather than playing a direct role in limiting the expansion of the blue crab, may enhance eutrophication-related phenomena such as bottom hypoxia and anoxia. In fact, on July 28, 2024, in the northern part of the Po Delta, in the locality of Boccasette, municipality of Porto Tolle, province of Rovigo in Italy, a massive die-off of blue crabs was registered, linked to the deoxygenation of the coastal bottom waters (ANSA, 2023). Such events, although sporadic and localized, could have a significant impact on the survival and reproduction of the crabs and requires further investigation.

Another important factor influencing the life cycle of blue crab is salinity. This species is adapted to a wide range of salinities, with its abundances varying along salinity gradients depending on sex and developmental stage (Hines, 2007). Salinity plays a key role in the blue crab reproduction, as the brackish conditions of estuaries and lagoons are essential for courtship and mating (Epifanio, 2019). However, our large-scale dataset limits our ability to explore finer-scale relationships that might explain local blue crab distribution and abundance patterns, such as the dynamics of estuarine systems. Local outlet geomorphology, river discharge and salinity distributions affect crab biology (e.g., abundances and distribution of sexes and reproductive stages). Unfortunately, this detail is unavailable from GBIF data.

The importance of estuaries and brackish lagoons in determining the presence of blue crabs is indirectly supported by our analysis, which shows an increase in blue crab populations associated with nutrient loads and shallow waters, conditions typical of these estuarine environments (Mancinelli et al., 2017a). We found that blue crab occurrences were generally higher in areas with shallower depths, such as the northern Adriatic Sea. An exception was observed in terms of proportional occurrences, as the depth was retained in the best model with, a weak and positive association. This was probably due to the lower human population density in deeper areas of the Mediterranean Sea, such as around islands and close to the Calypso Deep, its deepest point, rather than an actual increase in blue crab populations in deeper waters. To our knowledge, the only blue crab found in deeper Mediterranean waters was a female collected at 220 meters depth off the northern part of Imbroz Island, in the eastern Mediterranean Sea (Daban et al., 2016). In their native range, blue crabs exhibit depth-specific habitat uses, with females migrating to deeper waters for egg development and juveniles preferring shallower, vegetated areas (Hines, 2007; Eggleston et al., 2015; Ogburn and Habegger 2015). While our study did not differentiate between sex or life stage, it is evident that shallow waters, resulting from a gentle slope of bottom, can support blue crab populations in the Mediterranean by providing suitable habitat for various life stage. Shallow waters contain seagrass beds, oyster reefs, or intertidal sand flats that support dense populations of shellfish which, in turn, provide food and refuge habitat for blue crabs at various life stages (Cheng et al., 2022). Moreover, depth may play another important role in determining trophic availability. For example, shallow areas support relatively high benthic primary and secondary production in estuarine systems (e.g., Kemp et al., 2005), thereby providing food resources for blue crabs. Interestingly, nitrate concentration and depth were the most important predictors of blue crab occurrence. In general, phytoplanktonic chlorophyll is widely regarded as an indicator of coastal primary production, given its role as the foundation of the pelagic food web and its direct and indirect influence on benthic primary production. However, the density of blue crabs is more likely to be determined by benthic primary and secondary production, particularly in shallow waters where a consistent amount of solar radiation reaches the bottom. In these conditions, water column nitrate is the most important nutrient source for benthic diatoms and other micro- and macroalgae. This contribution to benthic primary production is not retained directly by chlorophyll a level in the water column; however, it is of paramount importance to the richness of the benthic consumers and thus to the abundance of the blue crab. Under this hypothesis, it was not surprising that the best model selected by AIC was explained by depth and nitrate concentration which, in turn, were related to the presence of river estuaries and lagoons, essential ecosystems for blue crab growth and reproduction. Excessive nitrate concentrations, which can indicate anthropogenic activities like nitrogen loads from diffuse agricultural sources and wastewater point sources (e.g., Viaroli et al., 2018), could provide less favorable conditions for blue crab populations, particularly when dissolved oxygen concentrations suffer.

Human population and shipping density were not main factors in explaining variation in blue crab occurrences. Thus, within the scales of observation, the blue crab occurrences were independent of the human population size and shipping in each area. The negligible influence of human population density on blue crab occurrences was further supported by similar results obtained using population size-corrected data (i.e., proportional occurrence). On the other hand, the low genetic diversity within Mediterranean blue crab populations suggests that the Mediterranean invasion was not a gradual process, but it likely resulted from a few isolated colonization events involving a small number of individuals (Schubart et al., 2023).

Shipping and ballast waters might have played a role in the initial phases of the blue crab invasion, providing continuous introduction of blue crab larvae within and across the Mediterranean Sea, potentially determining the pattern of distribution (Marchessaux et al., 2022). However, the subsequent spread in the Mediterranean has likely been influenced by other factors such as water mass dynamics and habitat connectivity, and suitable environmental conditions (Jeschke and Strayer, 2005).