Time dynamics of the mollusk communities
In this section, the presence-absence of the species recorded in the three periods (T1, T2, T3) are analyzed in relation to time, habitat, and human impact. The list of the 28 species of freshwater mollusks (17 gastropods and 11 bivalves) in T1–T3, their codes, and origins are given in Table 1.
The number of mollusk species has increased in time as the river has shifted from lotic to a mixture of flowing and stagnant sectors due to the building of reservoirs. T1 was characterized mainly by rheophilic elements and prosobranchs. Some species became extinct during the hydro-technical works (before or during T2) and are unlikely to recover, such as the rheophilic Theodoxus transversalis and Lithoglyphus naticoides. Other rheophilic species disappeared between T1 and T2 but managed to survive in tributaries and repopulated some sectors during the last years. The most remarkable recovery is that of the thick-shelled river mussel Unio crassus, a species protected by EU legislation. T2 was characterized by some extinctions but also colonization by lentiphilic pulmonates and tolerant, resistant species such as some clams. A few lotic species also survived in the river sectors between the dams. In T3, we encountered a rich and diverse community, including some newly established populations of AIS and the discontinuous presence of both lentic and lotic communities. Overall, the present-day fauna is richer than in former periods, consisting of 15 species of gastropods and 8 bivalves. The AIS included the gastropods Physa acuta and Ferrissia californica, which arrived in the area most likely during the XXth century, Viviparus acerosus, which is native to the Danube, but unknown until after 2000 in the upper-middle Olt River basin, the bivalves Dreissena polymorpha, also native in the Danube but an invader in the middle Olt since 2008–2010, Sinanodonta woodiana, first found in 2015, and Corbicula fluminea, which was first discovered in the Olt (and also in Transylvania) during our survey in February 2020. The mean number of native species per river’s sector increases almost linearly (2.8 species per sector in T1, 3.3 in T2, and 4.6 in T3), while the corresponding values for AIS increase non-linearly (no AIS in T1, 0.6 species per sector in T2 and 3.2 in T3).
In the CCA of freshwater mollusk community changes through time (Period as predictor), the adjusted explained variation was 23.6% (test on all axes, pseudo-F = 5.9, p = 0.001). The polygons delimiting the positions of the sites during the three periods of time show no overlap, and they were distinct and separated in the ordination space (Fig. 1a). T2 (the period with maximum human impact) is distinctly placed and separated from the period without impact (T1) along both ordination axes. Meantime, T3 is closer to T1, having an intermediate position between the other two periods, showing a trend of recovery, such as the return of some species. In the CCA of T1–T3 species presence-absence predicted by the selected environmental descriptors (Period, Habitat, and Impact) (Fig. 1b), the adjusted explained variation was 28.36% (test on all axes, pseudo-F = 4.2, p = 0.001). FD(Rao) computed on all FT was plotted as isolines by GAM on the ordination space (model AIC = -17.19, model test F = 5.1, p = 0.003; tests of non-linearity in predictor effects: F = 3.9, p = 0.03). The functional diversity decreased from T1 to T2, then increased sharply to T3; it also decreased from rivers (R) to lakes (L) and along the human impact gradient (Impact).
Canonical correspondence analysis (CCA) of mollusk communities: (a) classification diagram of sampling sites based on period (as predictors): T1—XIXth century, T2—1995–2000, T3—2020 (adjusted explained variation 23.6%; first axis accounts for 17.6% the second for 6.0%, both axes are significant, p = 0.001); (b) CCA diagram of species occurrence constrained by environmental predictors (period, habitat: L—lakes, lentic sector in reservoirs, R—river, lotic sectors, and Impact—human impact) with functional diversity expressed as Rao quadratic entropy index (FD (Rao)) isolines plotted by generalized additive models (GAM) on the ordination space (adjusted explained variation 28.36%; first axis accounts for 16.3%, the second for 6.0%, both axes are significant, p = 0.001) .
In the dc-CA with the selected predictors on T1–T3 presence-absence data, the first two axes separate the communities by period, each positioned in a distinct quadrant (Fig. 2). After a decrease in diversity from T1 to T2, in T3, there were more species and higher functional diversity. In time, there was a reduction in body size, a switch from species with separate sexes to hermaphrodites, a transition of oviposition towards ovo-viviparity (in snails), and external fecundation (in bivalves), and a switch of the feeding type. The dc-CA adjusted explained variation was 16.47%; tests based on sectors and species showed significant relationships (combined test for all axes, pseudo-F = 2.6, p = 0.006), the dimensionality test based on case scores was significant for the first axis (pseudo-F = 4.2, p = 0.001) and marginally significant for the second one (pseudo-F = 1.1, p = 0.053). In contrast, the dimensionality test based on species scores was significant only for the first axis (pseudo-F = 1.6, p = 0.004). The adjusted variation explained by environmental predictors (Hab, Impact, and Period) was 28.36%, and by the selected functional traits (Sexes, FeedT, SizeM, and Ovipos) was 14.64%.
Double-constrained correspondence analysis (dc-CA) with selected predictors on presence-absence data in T1–T3. The selected functional traits (in blue) are Sexes (circles: H—hermaphrodite, S—separate sexes), Feeding type (squares: SCR—scraper, SS—scraper and sediment, SF—scraper and filter, F—filter, SEDF—suspension and deposit feeder), Oviposition (diamonds: OV—ovo-viviparity, CAP—capsule/eggmass, BE—parental care, juveniles in brood pouches of demibranchs, No—no oviposition, external fecundation), and mean body size (SizeM); the selected environmental predictors (in red) are time (Period, with levels T1—XIXth century, T2—1995–2000, T3—2020), habitat (R—river, lotic sector; L—lake, a lentic sector in reservoirs) and human impact (Impact). Species are coded by the first three letters of the genus and species names. The adjusted explained variation was 16.47%, the first axis accounts for 12.7% and the second for 2.2%. Native species have black labels, while aliens (AIS) are written in green.
We have split the binary data describing communities into two parts: natives and AIS, using the latter as predictors. We partitioned the variation in native species composition explained by the three predictor groups (Period, Environment, and AIS) (Fig. 3), subjecting the explanatory variables to an interactive forward selection procedure. We used RDA with centered response variables (CCA can not be used because the empty rows in some tables hinder the use of a proper hierarchical permutation scheme). The adjusted explained variation was 39.6% (the simple effects: time accounted for 22.33%, habitat and impact 24.73%, and the selected AIS 20.82%). All simple and unique effects were significant (p < 0.01). Among the invasive species, S. woodiana, V. acerosus, and P. acuta had a significant effect, while D. polymorpha had a marginally significant effect (padj = 0.08). The unique effects were similar, with small differences between the environmental predictors (6.4%), AIS (6.3%), and time (5.9%). The highest value was registered for the shared effect of all predictors (7.3%), indicating a certain correlation between them.
Structure of native communities explained by variation partitioning among time (Period: with levels T1—XIXth century, T2—1995–2000, T3—2020), environmental features (human impact—Impact, and habitat—Hab: R—river, lotic sector; L—lake, a lentic sector in reservoirs) and the selected alien invasive species—AIS (Drepol—Dreissena polymorpha, Sinwoo—Sinanodonta woodiana, Vivace—Viviparus acerosus, Phyacu—Physa acuta). RDA with hierarchical permutation scheme was used, interactive forward selection of predictors; both simple and conditional effects are significant (p < 0.01).
Actual (T3) mollusk communities’ response to the environment, space, and alien invasive species
The following results refer only to species counts data from 2020. In the dc-CA on all species (Fig. 4) with selected environmental predictors (Dis_dam, Flow, Impact, Hab) and functional traits (SizeM, Nooff, Ovipos, FeedT, Sexmat), adjusted R2 = 42.64% for functional traits, adjusted R2 = 33.95% for the environmental predictors, and dc-CA adjusted R2 = 36.6% (combined test on sectors and species on all axes pseudo-F = 5.3, p = 0.001, on first axis pseudo-F = 0.5, p = 0.026). The effects of all selected predictors were significant (p < 0.05). The dimensionality tests based on case scores and on species scores were significant for the first four axes (p < 0.04). The human pressure rises in reservoirs toward the dams. It is associated with species having reduced number of offspring, reduced sizes, and oviposition in egg capsules, grouping the lentiphilic, pulmonates, and phytophilic species on the left side of the figure. The right side groups are mainly bivalves, rheophilic species, breeders, and species with an increased number of offspring and large body size. Some clams (Pisiidae) prefer canals, while naiads (Unionidae) occur in larger habitats, especially the riverbed.
Double-constrained correspondence analysis (dc-CA) on species counts in 2020 (T3), with selected environmental predictors (in red: Dis_dam—distance to the nearest dam downstream; Flow—average water flowing conditions; Impact—human impact intensity; Habitat: C—canals, R—river, lotic sector, L—lake, lentic sector in reservoirs) and functional traits (in blue: SizeM—mean body size; Nooff—number of offspring; Sexmat—sexual maturity; Feeding type, as squares: SCR—scraper, SS—scraper and sediment, SF—scraper and filter, F—filter, SEDF—suspension and deposit feeder; Oviposition, as diamonds: OV—ovo-viviparity, CAP—capsule/eggmass, BE—parental care, juveniles in brood pouches of demibranchs, No—no oviposition, external fecundation). Species are coded by the first three letters of the genus and species names. Selected 18 species are depicted. The adjusted explained variation was 36.6%, the first axis accounts for 11.8% and the second for 10.4%. Native species have black labels, while the aliens (AIS) are written in green.
We split the response data table into two components: natives and AIS. In the CCA with native species, the conditional effects of habitats Hab.C (canals) and Impact were significant (Table S1). For AIS, we used RDA since their gradient is short (2.8 standard deviation units), requiring a linear instead of a unimodal model. In the RDA with AIS, the conditional effects of Dis_dam and Flow were significant (according to the p-adj values, they were marginally significant, i.e. p-adj < 0.096). Thus, natives and AIS share no common predictors, indicating that the species composition of native and alien communities is linked to different environmental drivers (Table S1).
Analyzing relations between functional diversity and human pressure (Fig. S1), we have learned that AIS communities show a monotonic increase of functional diversity with impact, while the autochthonous communities respond by a sigmoid curve, with an accelerated increase at mid-values of impact, followed at higher values by an asymptotical reduction and then a decrease.
Since dams (and reservoirs) are the most prevailing sources of human impact in the research area, and their effect is a major question the present paper addresses, we explored the species response curves related to the distance from the nearest dam placed downstream (Dis_dam), using GLM (Fig. 5, Table S2). Natives and AIS show mostly independent segregation of abundance curves in response to continuously modified conditions along this gradient (flow, sediment deposits, presence of plants along the banks), resulting in the coexistence of ecologically heterogeneous elements. There is a continuous range of responses, adaptations, and tolerances to the conditions associated with the flow and distance to the dam. There are no apparent groups of species but rather overlapping distributions of counts related to the gradient.
Species response curves to the Dis_dam: distance (in m) to dams (count data, generalized linear models, only a selection of significant models illustrated; the model parameters are given in Table S2; alien species are depicted in green).
We explored relationships between FD of the native community and AIS counts. FD(Natives) computed on all traits was found to be negatively related to the abundance of only two invasive bivalves. The relationship was significant for C. fluminea (pseudo-F = 4.6, p = 0.045, r2-adj = 16.0%), and marginally significant for S. woodiana (pseudo-F = 3.7, p = 0.062, r2-adj = 12.62%) (Fig. 6a,b). In contrast, taxonomic diversity measures of native community (Fig. 6c,d) were positively and significantly related only to two gastropods: P. acuta (SNat and H; pseudo-F = 12.8, p = 0.001, r2-adj = 38.23%) and F. californica (SNat; pseudo-F = 8.1, p = 0.006, r2-adj = 27.17%). All the other relations with AIS were not significant.
The t-value biplot (Van Dobben circles) for the response of Rao functional diversity of the native communities (FD(Natives)) to alien invasive species counts (AIS), with (a) significant negative effect of Corbicula fluminea (Corflu) and (b) marginally significant effect of Sinanodonta woodiana (Sinwoo), and for the response of taxonomic diversity measures (SNat—number of native species, H—Shannon index, N2—Hill’s measure) of native communities in relation to AIS counts; two species of gastropods, (c) Physa acuta (Phyacu) and (d) Ferrissia californica (Fercal) had a significant positive effect. All the other relations of response variables with AIS were not significant.
We also performed a variation partitioning of H and A, using E and S as explanatory variables, and a reverse analysis by inversing the roles of these matrices (Table S3). The null hypothesis states that the response variables are similarly distributed against the classes of variance explained by the predictors.
(1) A and H as response variables: the test statistic χ2 is 13.18 if A is the expected distribution and 16.17 if H is the expected distribution. Both are higher than the critical value for df = 3 and p = 0.01, thus, the two distributions are significantly different. Both A and H are better explained by the unique effect of S, followed by the unique effect of E, and then by the shared effect of S and E, but the difference is given by the residual variance, which is higher, almost double, in A (26.9%) than in H (14.1%).
(2) E and S as response variables: the χ2 test statistics (0.62 in both analyses) is less than the critical value for df = 3 and p = 0.05, thus, the two distributions are not significantly different. Natives and AIS predict the distribution of the variability in environment and space in the same manner: unexplained variation is highest in E and S, followed by the unique effect of H, the shared effect of H and A, the smallest values being estimated for the unique effect of A.
The complete extended novel algorithm for the cumulative variation partitioning (CUVARP) is given in the Supplementary Information. The algorithm for partitioning the variation in the native communities (H) explained by alien communities (A), environment (E), and space (S) is synthesized in Table S4, and intermediary results are given in Table S5 and illustrated in Fig. S2 to S5. The percentages (CUVARP diagram) are represented in Fig. 7 and Fig. S6.
CUVARP diagram (Cumulative Variation Partitioning with Multiple Response and Predictor Matrices). All datasets are synthesized through selected PCoA axes scores of H—native species, A—alien invasives, E—environmental parameters, and S—space. No overlapping colored circles depict unexplained (residual) variation. The percentage of variation (explained and residual) from the total variation in all matrices is shown (the null value of AES overlap is not represented). Graphic alternatives are given in the Supplementary Information.
The total variation in our data (PCA of all matrices—H, A, E, and S) was 160.0, of which the unexplained variation (residual variation of selected PCoA axes of H, A, E, and S) was 56.442 (35.28%). The explained variation was 103.558 (64.72%). The unexplained variation in the AIS community was higher than the corresponding part of the natives, indicating relative independence of AIS from the measured spatial and environmental predictors. The unique variations of A related to E and S were also lower than those corresponding to H.
A variant of this method also includes the FT (for the natives and the aliens), and it is introduced and depicted in Fig. S7 as the CWM-AVARP (Community Weighted Means—Average Variation Partitioning with Multiple Response and Predictor Matrices) diagram.
Source: Ecology - nature.com
