Our study included more than 1.5 million km of lotic systems (rivers, creeks, channels, ditches) and over 50,000 of lentic ecosystems (lakes, ponds, reservoirs) in continental Europe spanning from 15°W- 45°E and 34°N -70°N (Figure 2).

Four approaches were used to obtain occurrence data of S. woodiana in Europe. The first approach was aimed at on retrieving peer reviewed and grey literature with known occurrences of S. woodiana in Europe. The search consisted of three searches using the ‘Web of Science’ and ‘Google Scholar’ search engines. The first search was performed using Web of Science and focussed on published scientific articles. A total of 57 hits were retrieved using the search term ‘Sinanodonta woodiana’. The second search was performed using Google Scholar to also obtain the grey literature and reports with data on its distribution. The search yielded 4080 hits of which the first 100 were checked for relevance. To ensure that specific local publications and/or reports were not missed a third search was performed using Google Scholar with a search term consisting of the species name with each European country’s name in its native language. Fifty-three separate searches were performed to cover the multitude of countries and regions on the European continent of which the first 10 hits were scanned on relevance resulting in 467 hits that were considered. In total, 622 hits including year of occurrence were retrieved and all data were georeferenced. All sightings and relevant information were subsequently entered into a database. The second approach was regularly updating the database with new publications published in the period 2018–2021. The third approach focussed on contacting malacological experts throughout Europe with the aim of retrieving unpublished locations with S. woodiana occurrences. Also, verified data sent by malacological experts from citizen science were included. The fourth approach consisted of using a backwards snowballing technique to search for any additional known occurrences based on acquired scientific literature using the 1^st, 2^nd and 3^rd data acquisition approach. Using the four approaches a total of 1322 georeferenced occurrences of S. woodiana for the European continent were compiled for the period between 1979–2021.

Nineteen bioclimatic data for three time periods (recent: 1970–2000, future: 2041–2060, and 2081–2100) and three shared socio-economic pathways (hereafter SSP, represented as the radiative forcing measured as watt/m²: 2.6, 4.5, 8.5) were downloaded from WordClim global climate database (www.wordlclime.org) with a spatial resolution of ~ 3 km (Table 1). SSP’s represent scenarios of projected socio-economic global changes and are used to derive greenhouse gas emissions with different climate policies (IPCC 2021). To account for multicollinearity between variables a correlation matrix was used to exclude variables with a Pearson correlation coefficient r ≥ 0.8 (Elith et al. 2010). Additionally, 11 environmental variables (Table 1) that likely have an influence on the distribution of S. woodiana were derived from online sources. For instance, geologic strata likely affect the chemistry of the groundwater, while the type of landcover may affect runoff, sedimentation, surface water chemistry and temperature (Liberoff et al. 2019; Hamid et al. 2020).

We have analyzed the relationships between the distribution of S. woodiana across Europe, defined by spatial coordinates (latitude and longitude) and time, given as the year of sampling (as reported by samplers or deduced from the published literature). Reported occurrences and cumulative values have been related to time by scatterplots and regression analyses. Relationships between the time and coordinates have been characterized first with trend-surface polynomials using Redundancy Analysis (RDA) on time constrained by centered coordinates (denoted Xc for longitude and Yc for centered latitude). The centering was done by subtracting the mean value from each term for reducing dependency among polynomial terms, as recommended by Šmilauer and Lepš (2014). Since polynomial terms reveal only monotonous or smooth non-linear changes in dependent variables, we have also used the method of distance-based Moran’s Eigenvector Maps (db-MEM) as described by Legendre and Legendre (2012), formerly known as Principal Coordinates of Neighbor Matrices (PCNM) as described in Šmilauer and Lepš (2014) in its use according to Ter Braak and Šmilauer (2018), because this method might reveal more types of spatial heterogeneity patterns occurring at multiple spatial scales. For this, the table of coordinates was subject to a db-MEM analysis, using the Euclidean distance for calculating sample distances, the cut-off threshold value being determined for the nearest neighbor. From the Principal Coordinates Analysis (PCoA) results, all axes (eigenvectors) with positive eigenvalues were used. The axes were saved and used as explanatory variables in a new RDA with time (Year) as the response variable, with an interactive forward selection procedure, using the values of probability adjusted (p-adj) by false discovery rate (FDR) as criteria of significance. The results of all constrained analyses included the adjusted explained variation, and the significances were tested by the Monte Carlo permutation test with 999 unrestricted permutations. The selected PCoA axes were subject to a reverse analysis, using time (year) as predictor within a new RDA, searching for patterns through T-value biplot (also known as Van Dobben circles) and by regression analyses between the PCoA axes and time, using generalized additive models (GAM). Then, the same PCoA axes were subject first to a Principal Components Analysis (PCA), and then to a RDA constrained by time. The scores from the first axes (unconstrained and constrained) were subject to a GAM analysis with Gaussian response distribution and model testing and selection by Akaike Information Criterion (AIC). The results were depicted in a combined contour plot against the latitude and longitude. For analysing and modeling relations between space and time we have used the software Canoco 5.15 (Ter Braak and Šmilauer 2018).

The maximum entropy (MaxEnt) approach was used because it is particularly efficient when handling complex interactions between response and predictor variables (MaxEnt, version 3.3.3.; Elith et al. 2006; Phillips et al. 2006). MaxEnt uses presence-only data and environmental predictor variables over a defined area of concern. The algorithm compares the presence locations to many background points, i.e., locations where the presence of a species is unknown. It predicts species occurrence by finding a probability function that is closest to uniform or most spread out considering the limits of the environmental variables at the inputted presence locations (Elith et al. 2011). We used 70% and 30% of the S. woodiana presence data for model training (calibrating) and model testing (validating), respectively. Five replicates were done in each model run with subsampling as the replicate run type. The random seed function was used to select different test/training samples in each replicate. To account for sampling bias, we limited the spatial extent to extract background points to those countries where S. woodiana findings were confidently verified. The predictive performance of the models was tested using the Area Under the Receiving Operator Curve (AUC) and the True Skill Statistics (TSS). AUC was used to and is defined as the area under a plot of the proportion of true positives versus the proportion of false positives. A minimum AUC value of 0.5 indicates that the model did not perform better than random while values close to 1 indicates perfect prediction (Swets 1988). TSS compares the number of correct predictions, minus those attributable to random guessing and accounts for correct predictions of a species presence and absence (Allouche et al. 2006); values close to -1 represent a model that is no different than random, while values close to 1 indicate excellent model performance. The 10^th-percentile threshold rule (excludes 10 per cent of the locations having the lowest predicted value) was applied as it is commonly used for data sets collected over a long time period and by different observers and methods (Radosavljevic and Anderson 2014; Bosso et al. 2017). Permutation importance was calculated by randomly permuting each variable among the training points (presence and background) to determine the decrease in training AUC. Jackknife tests were used to determine a heuristic estimate of each predictor variable on the overall distribution by iteratively excluding each variable and by using each variable in isolation (Phillips and Dudik 2008). To account for sampling bias (i.e., S. woodiana absence in some areas is because these areas were not yet or infrequently surveyed and not because environmental and climatic conditions in these areas prevent the occurrence of that species) a bias file was created that puts more weight on background data that were not or less frequently surveyed (Phillips et al. 2009). The seven model outputs were imported into ArcGIS 10.7.1 (ESRI, 2008) to create maps representing a predicted habitat suitability ranging from 0 to 1 for each grid cell. Different colors were used with red indicating a high probability of suitable conditions, green indicating conditions typical of those where the species is found, and lighter shades of blue indicating a low predicted probability of suitable conditions. We regrouped habitat suitability into three categories: least (< 0.35), moderate (0.35–0.65), and high (> 0.65).

The relationships between the spatial distribution of S. woodiana and time are depicted in Figure 3. After the year of the first report (1979), there were hardly any new locations and findings, but after a while, the reported occurrences increased unevenly and exponentially to a maximum of about 100 new locations (in 2018), followed by a tendency of decline during the last few years (Figure 3A). We have found significant non-linear regressions between reported occurrences (as well as cumulative occurrences) and time, given as years of sampling, ranging between 1979 and 2021. In the RDA (not shown here) with the year as response variable and polynomial terms of centered coordinates as predictors, seven selected polynomial terms (in decreasing order of explained variation: Xc, Xc³, Yc²Xc,Yc², Xc², YcXc, and Yc³) explained 6.99% in time variability (adjusted explained variation), all their effects being highly significant (p-adj < 0.05). This result shows that the relationship between the spatial expansion of the invasive species and time is significant but weak, explained mostly by longitude and less by latitude, considering the relative importance of the terms in the polynomial predictors. When considering their individual separate effects, time was independent of latitude (Yc, pseudo-F = 0.4, p = 0.5), while the longitude (Xc) had a significant and negative relationship with time (pseudo-F = 0.32, p = 0.001), the adjusted explained variation being low (2.29%).

The coordinates have been subject to a db-MEM analysis and the resulted PCoA axes scores (spatial eigenvectors) have been saved in a new data table. The PCoA axes (as predictors) relations to time (as response variable) have been investigated by RDA using an interactive forward selection procedure, resulting in a selection of 17 significant explanatory variables. Then, we did a reverse analysis using time as predictor and all these selected PCoA axes as response variables, for testing and illustrating their relations, by means of a T-value biplot (Van Dobben circles). Most small value PCoA axes (denoted in the figure as PCO 1, 5, 6, 8, etc., representing coarse spatial scales) have a negative relationship with time, meaning a decrease from 1979 to 2021, while the larger order (representing finer spatial scales) are mostly positively related to time (Figure 3B). The relations (GAM response curves) between the selected PCoA axes’ scores against time are depicted in Figure 3C. The axes characterized by a larger coefficient of determination (r²) are mostly monotonically decreasing in time. Some axes show a unimodal response, increasing slowly up to a point and decreasing afterward (such as PCO 3 and 7), while the rest are characterized by a very gentle slope, increasing with a low rate (such as PCO 162). Their coefficients of determination are also rapidly decreasing, from PCO1 r₁² = 16.1% (the maximum value) to about 7% (PCO5) and 5% (PCO8), then eight axes show values between 1 and 3%, while the rest of the axes do not surpass 1%. The coefficients of determination tend to decrease (but not monotonically) with the increase in the axes rank, from the maximum of PCO1 to the minimum value of r²₁₆₂ = 0.4% for PCO162. Most PCO axes intersect each other, and all with the abscissa, in about the same short sector, and the response value of zero is near 2010. PCoA selected axes scores unconstrained by a PCA analysis and constrained by Time (year) in an RDA analysis are depicted (Figure 3D) as contour plots on an XY attribute spatial plot (meaning in a space defined by latitude as Y and by longitude on X), by their scores along the first unconstrained (in black isolines) and, respectively, constrained axis (in red), using a GAM analysis (3 degrees of freedom for both X and Y, Gaussian response distribution and stepwise selection using AIC, significant model, both for the response R² = 35.8%, as well as for the explanatory case scores R² = 9.1%, p < 0.001). The plot resulting from the GAM analysis of the scores from the PCA and RDA first axes in relation to geographical coordinates shows similarities, meaning that there is a symmetrical and correlated change in scores of sampling stations of both analyses from a central core (located around the first introduced location, and broadly a core in the inner Carpathian Basin) radiating unevenly towards all directions. But there are also some dissimilarities, given by the smoother and more predictable distribution towards the east and rougher and more discontinuous towards the west, when considered in relation to time, and by the dispersal along the main axis of continental Europe, the east-west axis, contrasting with almost not related to time along the latitudinal, north-south axis.

For the seven models the following performance value ranges were obtained: training AUC (0.954±0.004 and 0.964±0.005); test AUC (0.917±0.003 and 0.924±0.003); TSS (0.677±0.019 and 0.690±0.012). Under the recent scenario, variables with the highest permutation importance were elevation, distance between water bodies, and mean temperature of the coldest quarter. Jackknife sensitivity analysis of the seven models showed that elevation, followed by mean temperature of the warmest quarter had the highest test and training gain and the highest AUC value when used in isolation. The variable that decreased the training gain most when omitted was either elevation (recent, SSP2.6_2050, SSP2.6_2100, SSP8.5_2050) or land cover (SSP4.5_2050, SSP4.5_2100, SSP8.5_2100). The variable that decreased the test gain most when omitted was either elevation (SSP2.6_2050, SSP2.6_2100) or land cover (recent, SSP4.5_2050, SSP4.5_2100, SSP8.5_2050, SSP8.5_2100). No variable decreased the AUC gain when omitted.

Under recent climate conditions, 2.3% (147,332 km²) of European watershed area is predicted as a highly suitable habitat for S. woodiana (Table 2, Figure 4). More than 75% of this area (113,566 km²) is located in the temperate climate zone between 40°N and 60°N and includes northern Italy, northern Serbia, southern and central Hungary, western, southern, and southeastern Romania, southern and central Poland, eastern Croatia, southern France, and southwestern Ukraine. Areas within a watershed predicted as highly suitable were associated with lowland streams and water bodies (elevation < 200 m), fluviatile deposits from porous sedimentary rocks, agriculture as the dominant landcover, a close distance to urban areas, a mean temperature of the coldest quarter between 2°C and 8°C (Figure 5). Major water bodies in these watersheds include but are not limited to: Po, Etch, and Reno rivers and their tributaries draining into the Adriatic Sea, Tiber and Arno rivers draining into the Mediterranean Sea as well as Alpine lakes including Maggiore, Garda, D’Iseo and Como in Italy; Tisza, Sava, and Danube rivers in Serbia; Danube, Tisza, Rába, Drava, and Mureş/Maros rivers, the Sió Channel, and Lake Balaton in Hungary; Danube, Mureş/Maros rivers, Olt, Someș, Criş/Körös rivers and the Danube Delta in Romania and Ukraine; Prut River in the Republic of Moldova and Romania; Vistula, Oder, and Warta rivers in Poland; Sava and Drava rivers in Croatia, and Rhône and Garonne deltas in France. The predicted suitable S. woodiana habitat in lakes and stream reaches of western Latvia and Estonia (especially lakes Busnieku, Engure, Liepajas and Usma in Latvia and Lake Mullutu-Suurlath at Saaremaa Island in Estonia) would make this the northernmost currently suitable location in Europe. The boreal climate zone above 60°N and the southern Mediterranean climate zone below 40°N were classified as least suitable.

Under future climate conditions, the area of highly suitable habitat expanded slightly in the mid-century (2050, SSP 2.6: 2.4%; SSP 4.5: 2.2%) and then decreased by the end of the century (2100: SSP 2.6: 2.2%; SSP 4.5: 1.7%). For scenario SSP 8.5, areas of highly suitable habitat decreased to 2.2% for both time periods (Table 3). Highly suitable habitat remained along an east-west gradient stretching from northeastern Spain to the Danube Delta in the East between 43°N and 50°N. Our results did not show an expansion of suitable habitat into areas below 40°N and above 60°N in the future. Areas of low habitat suitability contracted under the three climate scenarios and time periods, while areas of medium habitat suitability increased. An exception was scenario SSP 4.5 at the end of the century where areas of low and high habitat suitability expanded and contracted, respectively (Figure 6, Suppl. material 1).

Identifying, recording, and publishing new locations of S. woodiana on the European continent showed an exponential model in time. In the last years, while the minimum number of newly identified locations is still increasing (the bottom edge of the scatterplot in Figure 3A), the upper limit displays a tendency of not surpassing the maximum of about 100 occurrences per year, repeatedly achieved. Following an increased interest in documenting the species’ invasive dispersal, after reaching a peak, the number of new reported occurrences likely will soon start to decrease. All the analyses of relationships between space and time, and their results (some depicted in Figure 3) have the same problem: the impossibility of discriminating between the subjectivity in reporting occurrences (when the species was searched, intentionally and probably most of the time unintentionally, by someone and reported accordingly) and when and where the species has objectively invaded. These two features (human-related and species distribution) are considered together in these analyses, and there is no way to separate them. We can only consider (based on our experience and judgement) that these two are strongly and positively correlated. We are also aware of the problem of relying only on reported locations, after (indetermined time) the species has invaded that particular spot. Most of the time, we don’t have access to reliable information of studies in locations where the species has not invaded and after this has happened. Thus, all the inferred information relies on occurrence (presence) data rather than presence-absence data. All this adds another source of indetermination to our study.

Small ranked axes (lower order axes of PCoA obtained during db-MEM and further analyses) represent eigenfunctions describing continental patterns (broad scale), while the larger the order, the finer the spatial eigenfunction represents (Legendre and Legendre 2012; Šmilauer and Lepš 2014). This means that continental (large scale) and regional scale invasions (Figure 3B, C) are negatively related to time (they have decreased during the invasion timespan). In contrast, the local (small or fine scale) dispersal is still ongoing, but this is unevenly (some higher order axes have positive, some negative responses) and loosely (very small coefficients of determination, and the larger the order, the lower the value) related to time. All these hint at the saturation of continental and regional expansion, while locally the species still invades new sites. The turning point of invasive behavior probably happened around the year of 2010, which is the point where all the PCoA axes collide with each other and the abscissa (Figure 3C). The relative acquisition of new points at the local scale is relatively constant in time.

Since S. woodiana was first encountered in what is usually considered Eastern Europe (an area in the Inner Carpathian Basin - however, from a geographical perspective, this is central-southern Europe), and most studies are conducted in central and western Europe, it is not surprising that the pattern of dispersal (Figure 3D) is developing across the longitudinal axis (and along a latitudinal belt), and from the core to all directions, but especially towards west. Longitude (i.e. reported occurrences related to longitudinal coordinates) is also significantly related to the time, while latitude is not, as our analyses revealed. The longitudinal axis of Europe, along the temperate zone, has historically proved to be a factor favoring the expansion of humans, technology, agriculture, farming, languages, and much more (Diamond 2017), in sharp contrast with continents characterized by north-south dominant geographical axes. Being a (meso-) thermal species with tolerance and possibly also a preference for higher temperature values, S. woodiana showed a predominant expansion along the temperate area and mean latitudes. The unconstrained PCA axis 1 scores (Figure 3D in black) show an asymmetrical rhomboidal shape of waves, with broad width and much narrower height, which resembles the geographical features of Europe’s main axis, but also climatical preferences of the invasive species. The red lines (time constrained RDA axis 1 scores of selected spatial eigenvectors) show the time-related invasive behavior (and its documentation by reports, which bring the subjectivity in this analysis) that is correlated (similar shape and radial waves from the same core) but not identical to the unconstrained contour plot. Limiting factors, and as such, edges of future distribution, are probably different across the latitude. The most prevailing limiting factor towards the north is probably temperature, while certain hydro-geographical features of landscapes and riverscapes might be limiting factors to the south.

Threshold-independent measures of model accuracy, i.e., AUC for recent and future scenarios were above the value of 0.9 indicative for excellent predictive power (Swets 1988), and threshold-dependent measures, i.e., TSS for all models indicated very good model performance (Allouche et al. 2006). One assumption of species distribution models is that the sampling effort is randomly distributed over the area of interest. However, in reality species distribution models rely on occurrence records which are spatially biased as some areas are easier to access and thus more frequently surveyed than others (Phillips et al. 2009). Sampling bias may result in inaccurate models if suitable areas are not included, or some areas are over-represented due to more frequent sampling (Fourcade et al. 2014). To circumvent this issue, the inclusion of a bias file as suggested by previous studies not only increased the predictive performance of the model but also decreased the risk of omission errors (Phillips et al 2009; Kramer-Schadt et al. 2013; Syfert et al. 2013).

Based on our model predictions, S. woodiana can occupy a broad range of lotic and lentic habitats in the European temperate and northern Mediterranean climate zones ranging from northern Turkey and Greece in the South to the Baltic states in the North and from western France to southern Ukraine in the East. The distribution of S. woodiana is controlled by both habitat variables, such as elevation and land cover, and by climate variables such as the mean temperature of the coldest quarter. Highly suitable habitat was associated with lower-elevation watersheds dominated by agriculture and urban areas. This is in line with previous studies, which showed a preference of this species for lowland freshwater bodies, such as lakes and ponds or slow flowing rivers and muddy riverbeds without strong currents (Sárkány-Kiss et al. 2000; Sîrbu and Benedek 2018; Sousa et al. 2021). It is well known that watersheds dominated by agriculture and urban areas are characterized by elevated nutrient levels, thermal pollution and higher sediment loads due to soil erosion. Sîrbu and Benedek (2018) documented that S. woodiana is often associated with anthropogenically altered water bodies and in fact it is sometimes the only Unionidae present due to its higher tolerance for hypoxia (Sîrbu et al. 2005) or thermal (Kraszewski and Zdanowski 2007; Labecka and Domagala 2018) and organic pollution (Bielen et al. 2016). The ‘distance between water bodies’ was ranked a variable with high permutation importance indicating a faster spread of adult mussels or fish infested with glochidia larvae among water bodies close or connected to each other. For instance, fishponds seem to be a ‘hotspot’ for S. woodiana dispersal due to adequate food availability, the high density of host fish, production, and sale of fry infected with glochidia (Urbańska et al. 2019, Figure 7). They are mostly connected to other lotic or lentic systems where ‘escaped’ infested fish enter or juvenile mussels are washed from the bottom sediment. Active fishpond drainage in inland fishery should also be considered as an influence on the distribution of Chinese pond mussels (e.g., transfer of still-living mussels to the water gardens). Temperature variables, such as the mean temperature of the coldest quarter, indicate that S. woodiana would be able to spread almost throughout the entire European continent, except a few areas including southeastern Spain, Norway, central and northern Sweden, Finland eastern Estonia, eastern Latvia, and eastern Lithuania, northern Ukraine, and Belarus. Previous research documented S. woodiana’s occurred in regions with air temperatures ranging from -31.1°C to +38°C (Urbańska et al. 2021). Both larvae and adult mussels have significantly higher thermal tolerances than native unionids (Protasov et al. 1993; Benedict and Geist 2021) and the filtration rate was experimentally shown to be higher for S. woodiana compared to Unio douglasiae under warmer temperatures (Kim et al. 2011). The species’ broad tolerance for a wide range of environmental factors coupled with its high reproductive capacity (Labecka and Domagala 2018, 2019; Labecka and Czarnoleski 2021) provides S. woodiana with competitive advantage over other unionids and usually becomes the dominant species after it establishes new populations (Urbańska and Andrzejewski 2019). Even after an acute human impact (e.g., a wave of pollution) that caused a die-off of all naiads (Unionidae), S. woodiana was the first and sometimes the only recolonizing naiad species (Sîrbu and Benedek 2018).

Under the ‘most optimistic’ climate scenario SSP2.6 there was an expansion of suitable habitat in the middle of the century. This corresponds with previous findings that aquatic alien species will likely benefit from the predicted increase in temperatures especially in northern latitudes fitting better with those in their native habitat (Rahel and Olden 2008). However, suitable habitats declined under scenarios SSP4.5 and SSP8.5 for both time periods. A reason might be climate-change related alterations of temperature and precipitation and therefore lower river discharges by 2050 in Europe during summer and fall season, especially in the temperate continental and in the Mediterranean zone (Schneider et al. 2013; van Vliet et al. 2013). Increase of water withdrawal for irrigation and electricity production could further reduce discharge (Döll et al. 2009). Furthermore, the reduction in discharge may restrict host fish movement and interrupt habitat connectivity (Baldan et al. 2021) which may further reduce the natural dispersal rate of S. woodiana. However, we expected a much larger expansion of S. woodiana into northern latitudes and into higher elevations in the future scenarios due to the predicted increase in temperature.

Our modeling approach comes not without caveats and its critical to acknowledge the limitations of our results. We may have not included variables (especially ecological) which may be important to the distribution of S. woodiana. For instance, previous studies demonstrated the importance of including data on the diversity and distribution of host fish into SDM (Daniel et al. 2017; da Silva et al. 2022). Douda et al. (2012) found that 75% of lotic systems in the Czech Republic would have suitable host fish. These authors also suggested that food availability likely be a limiting factor because S. woodiana was usually found in larger numbers and sizes in eutrophic water bodies. Therefore, our model would likely benefit from including data about the trophic status of lotic and lentic water bodies. Furthermore, the sheer amount of water bodies included in this study prohibited the measurement of in situ water quality parameters. Instead, we used available data sets as proxies. For instance, climate data and catchment geology were used as surrogates for water temperature and instream conditions, respectively. However, air temperature may not accurately track the temperature of larger bodies of water in particular. Although this approach has raised some criticism (Elith and Leathwick 2009), McGarvey et al. (2018) did not find significant differences in the prediction of fish distribution when climate versus instream covariates were used. We suggest that a variation partitioning between effects of fish communities and environmental (and especially habitat) descriptors on distribution and abundance of S. woodiana might be of great importance for understanding its ecology and dispersal strategy. The species’ effects and interactions with other naiads (Urbańska et al. 2021), with the native mollusks (as done by Sîrbu et al. 2022) and other freshwater communities should be also prioritized for a better understanding of its invasive behavior and its consequences.