The manipulated stone substrate experiment was conducted at the littoral area in Upper Lake Constance (N: 47°41′31.0″, E: 9°12′07.0″), one of the two major basins in Lake Constance, Germany. The littoral area has predominantly stony gravel and pebble substrate. It is surrounded by a semi-natural habitat, with scarce suburban housing and with hardly any extensive artificial land use (Gergs et al., 2011). The stone substrates were deployed at a depth of 0.80 m.

The experimental design followed the former study by Mörtl and Rothhaupt (2003). The standardized stone substrates (commercially available concrete cobblestones) had an upper surface area of 121 cm². The upper surface area of each stone substrate was manipulated to obtain four treatments with eight replicates each. The first control treatment (B treatment) consisted of blank and clean stones with no shells and no mussels on them. Treatment 2 (S treatment) consisted of stones with eviscerated shells of mussels that were glued together to mimic intact dead mussels (Botts et al., 1996). Shells were glued to the substrate at their posterior side and positioned in a similar direction, simulating field observations of adult dead mussels at the lakeshore. Treatment 3 (Q treatment) and treatment 4 (Z treatment) consisted of stones with upper surfaces covered by living adult quagga and zebra mussels, respectively. Living mussels of individual species were glued to the substrate at their posterior side. Care was taken to ensure that individual mussels could open and close their shells as in natural cases to ensure survival. The abundance of shells or living mussels on the upper surfaces of stones was equivalent to 3,719 ± 413 ind./m². The substrate appearance resembled that of a natural mussel bed or druse, simulating field observations of living mussels. Caution was taken to ensure that glued individual mussels or shells were approximately equidistant.

To prepare the stone treatments, adult quagga and zebra mussels of ca. 20 mm shell length had been collected by hand at a littoral site of the lake where the two dreissenids coexisted by then (i.e., June 2018). Based on the morphological traits, mussels had been identified and separated as quagga or zebra mussels (Pathy and Mackie, 1993; Kerambrun et al., 2018). To prepare the shell treatments, shells of eviscerated mussels were washed and dried in the laboratory at room temperature for 72 h. Stones were deployed at the study site for the experiment between June and July 2018 (28 days). Elsewhere, artificial substrates have been reported to reach stabilized newly settled macroinvertebrate communities within 20–60 days (Roby et al., 1978; Wise and Molles, 1979). Therefore, through the experimental period of 28 days, we envisioned to capture the invasion processes and evaluate their effects on the density and taxonomic richness of newly settled macroinvertebrates within the first 4 weeks (short time duration).

On the last day of the experiment, all the stones were collected from the lake, thoroughly washed, and carefully inspected for attached (newly settled) macroinvertebrates. The newly settled macroinvertebrates were gently removed with a brush, washed, and collected from each stone after filtering through a 200 μm mesh. We took all the newly settled macroinvertebrates in the entire area of the respective stone into account. Afterward, living adult quagga and zebra mussels that had been glued on the stones were gently removed from their respective stones, and their survival was assessed by observing shell gape and soft tissue. If no soft tissue was present inside the shells, then mortality was recorded. One stone in the Z treatment was upside down at the end of the experiment, and this sample was omitted from further analysis.

In the laboratory, all newly settled macroinvertebrate samples were sorted and counted under a dissecting microscope (× 10 magnification, Zeiss Stemi 2000-C). Identification was made to species or genus level, except for Diptera, which mainly included Chironomidae. Chironomidae were identified at the subfamily level due to morphological ambiguity at their smallest size ranges. The size of macroinvertebrates was measured using a microscope equipped with an ocular micrometer. Newly settled dreissenids larger than 1 mm in length were separated into quagga or zebra mussels based on morphological traits. We pooled the small dreissenids (≤ 1 mm in length) into one group.

The five common newly settled macroinvertebrate species recorded from the experiment, which included the two mussel species (quagga mussel D. rostriformis bugensis and zebra mussel D. polymorpha), two snail species (New Zealand mud snail Potamopyrgus antipodarum and faucet snail Bithynia tentaculata), and one gammarid species (killer shrimp Dikerogammarus villosus), were used for stable isotope measurement. The soft tissues of the two mussel species and two snail species were eviscerated for stable isotope analysis. Individuals of the gammarid species were used for stable isotope analysis.

Lipid content was removed from dried and pulverized samples (Post et al., 2007; Rothhaupt et al., 2014) and weighed (ca. 0.85 mg) in tin cups to the nearest 0.001 mg, using a micro-analytical balance. The ¹³C/¹²C or ¹⁵N/¹⁴N ratios of samples were determined using two instruments. Half of the samples were measured at the stable isotope lab of the University of Konstanz in Germany. Here, samples were combusted in a vario PYRO cube elemental analyzer (Elementar Analysensysteme, Germany) at 1,120°C. Resulting CO₂ and N₂ were separated by gas chromatography and passed into an IsoPrime isotope ratio mass spectrometer (IRMS; Isoprime Ltd., Manchester, United Kingdom). The rest of the samples were measured at the stable isotope laboratory of the University of Erlangen-Nürnberg, Germany. Here samples were combusted in a Costech Elemental Analyzer (ECS 4010; Costech International, Pioltello, Italy; now NC Technologies, Bussero, Italy). The resulting gases were separated and measured using Thermo Scientific Delta V plus IRMS (Thermo Fisher Scientific, Bremen, Germany). The repeatability and comparability of results obtained from the two laboratories were first ensured by measuring multiple macroinvertebrate samples (“dummy” samples) in both laboratories before proceeding with the samples used for this experiment.

All stable isotope values (δ¹³C or δ¹⁵N) were reported in the δ notation (per mill) where δ = [1,000 × (Rsample/Rstandard) –1] ‰, relative to the Pee Dee Belemnite (PDB) for carbon and atmospheric N₂ for nitrogen in parts per thousand deviations (‰).

Non-parametric Kruskal–Wallis test was used to analyze the differences in the abundances of newly settled macroinvertebrates among four different treatments, as the normality of the whole dataset was not reached (Shapiro–Wilk test). Dunn’s test was used for post hoc pairwise comparisons. Non-parametric Mann–Whitney Wilcoxon test was used to compare the abundances of quagga and zebra mussels in each treatment.

To identify differences in community composition across all treatments, we examined dissimilarities among macroinvertebrate assemblages by applying non-metric multidimensional scaling (NMDS) using Bray–Curtis dissimilarities. Mean abundance data (ind./m² of newly settled macroinvertebrate per taxon) were used in the analysis, and data were standardized by the decostand function using the method “total.” We then assessed whether the composition of the macroinvertebrate community differed among treatments using analysis of similarity (ANOSIM). In order to reduce distortion of macroinvertebrate assemblage in the community, taxa accounting for less than 0.5% of the total mean abundance per treatment were excluded from the analysis.

Graphical summaries of the mean abundance of abundant newly settled macroinvertebrates (ind./m² + SE) per treatment were made for taxa that can be identified up to the genus or species level, which included D. rostriformis bugensis, D. polymorpha, Dreissena (≤ 1 mm), P. antipodarum, B. tentaculata, D. villosus, Katamysis warpachowskyi, and Caenis (Figure 1). A summation of abundances for all the newly settled macroinvertebrates present in each sample is expressed as the total abundance.

The stable isotope values of newly settled macroinvertebrates were not normally distributed (Shapiro–Wilk test). Kruskal–Wallis test was used to analyze the differences in δ¹³C and δ¹⁵N values for the five main newly settled macroinvertebrate species among four treatments, and Dunn’s test was used for post hoc pairwise comparisons. Mann–Whitney Wilcoxon test was used to compare the δ¹³C and δ¹⁵N values of quagga and zebra mussels in each treatment. We used Bayesian standard ellipse area (SEA) and size-corrected standard ellipse area for small samples (SEAc) to describe species and community isotopic niche width of newly settled macroinvertebrates among different treatments. We also calculated and reported the 95% confidence intervals of SEA using 10,000 posterior Markov chain Monte Carlo draws (Jackson et al., 2011; Catry et al., 2016).

All calculations and figures were conducted using R (R Core Team, 2021) with the package readxl, tidyverse, DescTools, ggplot2, ggpubr, vegan, and SIBER, in order to import data, reshape data, conduct Dunn’s test, make figures, arranges figures, conduct NMDS analysis, and analyze Bayesian standard ellipse areas, respectively.

The most abundant taxa in all the four treatments were mussels [D. rostriformis bugensis, D. polymorpha, and Dreissena spp. (≤ 1 mm)], snails (P. antipodarum and B. tentaculata), gammarids (D. villosus), mysids (K. warpachowskyi), mayfly nymphs (Caenis spp.), and chironomid larvae (Chironominae and Orthocladiinae) (Figures 1, 2). Other taxa, such as Trichoptera, that occur infrequently were used only for total abundance analysis.

Differences in the abundances of newly settled macroinvertebrates among the four treatments are presented in Table 1. The abundances of newly settled quagga mussels were not significantly different. The abundances of newly settled zebra mussels were also not significantly different. The abundances of newly settled snail P. antipodarum and mayfly nymphs Caenis were significantly different (p < 0.01), with the B treatment having the lowest numbers. The abundances of newly settled gammarid D. villosus were less in B treatments than in Q and Z treatments (p < 0.05). The abundances of newly settled chironomid larvae of subfamily Chironominae were higher in Z and B treatments than in S and Q treatments (p < 0.05). Total abundances of all newly settled macroinvertebrates were significantly higher in Z treatments than in B treatments (p = 0.002).

The total abundance of quagga mussels in the entire experiment area was higher than the total abundance of zebra mussels. In Z treatments, the abundance of newly settled quagga mussels was significantly higher than that of zebra mussels (p = 0.021), and additional information is provided in Supplementary Table 1.

In the NMDS plot (stress 0.09), the average composition of colonizing macroinvertebrate community across treatments appeared to be most segregated by mussel treatments (ANOSIM global test R = 0.473, p = 0.001) (Figure 2). For this 28-day experiment, Q and S treatments were closely grouped and separated from the Z treatment, while the B treatment was more isolated on the plot.

The mortality rate of glued adult zebra mussels (24.00 ± 2.28%) was significantly higher than that of glued adult quagga mussels (7.33 ± 1.33%) (p < 0.01; Figure 3). The abundance of newly settled chironomid larvae (subfamily Chironominae) showed a positive correlation with the mortality rate of adult dreissenids (Figure 3).

The stable isotope values of five main macroinvertebrate species are shown in the stable isotope (δ¹³C and δ¹⁵N) bi-plots (Figure 4). Confidence intervals of the standard ellipse areas for different treatments and the whole experimental area are shown in the density plot (Figure 5 and Supplementary Figure 1).

Among four treatments, the δ¹³C values of newly settled quagga mussels were significantly higher in the B treatment than in the Z treatment (p = 0.016; Table 2). The δ¹³C values of newly settled snail P. antipodarum were significantly higher in Z treatment than in S treatment (p = 0.014; Table 2). The other three species showed no significant differences in δ¹³C values among the four treatments. For all five main species, there were no significant differences in δ¹⁵N values among the four treatments.

The SEAc results indicated that the newly settled quagga mussels showed a smaller isotopic niche in Z treatment than in B, S, and Q treatments (Figure 5 and Table 3). The newly settled zebra mussels also exhibited a smaller isotopic niche in Z treatment. The remaining three newly settled (frequent) macroinvertebrate species showed equivalent isotopic niches among the different treatments. Taken together, SEAc for the whole community was higher in B treatment.

For the two dreissenids, the δ¹³C values of newly settled quagga mussels were higher than the δ¹³C values of newly settled zebra mussels. This pattern was evident in all four treatments (p ≤ 0.01; Supplementary Table 2). In each treatment, the δ¹⁵N values of all newly settled quagga mussels showed no significant differences from zebra mussels. Based on SEAc results, there were no overlaps in the isotopic niches between quagga and zebra mussels in each treatment.

The reaction of newly settled macroinvertebrates toward the stone substrate treatments complimented the hypothesis of facilitative associations between invasive dreissenids (Ward and Ricciardi, 2007; Ozersky et al., 2011), wherein one species has a positive effect on the establishment of another species. For example, the highest mean density of macroinvertebrates was consistently found in substrate treatments with shells or living mussels. By implication, the presence of mussels or shells was driving significant variation across experiments. Given the uniform in situ physicochemical conditions for stone substrate, this variation could be attributed to differences in the presence or absence of zebra and quagga mussels and their respective shells. Indeed, such effects were anticipated, since macroinvertebrates have shown a positive reaction to the physical structure of mussel shells (Ricciardi et al., 1997; Burlakova et al., 2012).

The most frequent newly settled macroinvertebrate taxa in our study, P. antipodarum, D. villosus, and Caenis spp., were mainly present at comparatively higher densities on stone substrate shells (S) and living mussel treatments (Q and Z). The rarity of settling macroinvertebrates on blank stone treatments (B) clearly indicates the feasible benefit of macroinvertebrates from the effects provided. In all three cases, each taxon could indeed benefit from the increased complexity of the substratum and surface area associated with natural dreissenid beds (Stewart et al., 1998; Ward and Ricciardi, 2007). Briefly, for P. antipodarum, shells may provide a bed substratum on which tubular refuges with protrusion could be constructed for grazing and scraping on detritus and sediments (Broekhuizen et al., 2001). In particular, the zebra mussel has shown mutualistic interactions with killer shrimp D. villosus and can benefit them in various ways. As shown during ex situ laboratory experiments on dreissenids, killer shrimp evidently utilizes mussel beds as protection from fish and allopatric predators more efficiently than other gammarids (Kobak et al., 2014). Similarly, juvenile zebra mussels can attach to the hard cover of the killer shrimp which could facilitate their dispersal (Yohannes et al., 2017), while the production of feces and pseudo-feces by zebra mussels may provide food (Klerks et al., 1996; Gergs and Rothhaupt, 2008). For Caenis, the taxon’s size might have allowed the use of the physical structure of dreissenids as nearby refuges through the exploitation of small interstices that are typical characteristics of dreissenid beds (Werner and Rothhaupt, 2007).

Like most of the three numerically dominant species, the abundances of newly settled chironomid larvae of subfamily Chironominae were higher in response to the presence of living zebra mussels (Z treatments). However, this species was one of the few species that showed higher density in the absence of mussels (B treatment). Almost all the newly settled chironomid larvae in B treatment were rather small with the average width of head measuring less than 180 μm. The observed pattern is presumably linked to the natural temporal or seasonal life history patterns (in growth, size, and abundance of the species), and perhaps this response was evident due to the 4 weeks of methodological effect. The overall pattern of chironomid larvae abundance and richness from our experimental study reinforces previous findings that zebra mussels facilitate the growth of benthic macroinvertebrate communities. This is in complement with previous studies that showed that the abundance of chironomid larvae increased after the colonization of zebra mussels (Botts et al., 1996; Gergs and Rothhaupt, 2008). Nevertheless, the abundance of chironomid larvae in Z treatments was higher than in Q treatments (Table 1). The high mortality rate of glued adult zebra mussels than quagga mussels might contribute to the differences, as the animal remains are part of the diet for chironomid larvae (Johnson, 1987).

The results of our analysis using analog shell and live quagga mussels should well be noted again, wherein no significant difference in the abundance of all taxa between S treatments and Q treatments was observed. This finding suggests that macroinvertebrates might primarily respond to the physical structure of quagga mussel shells, as discussed in previous investigations (Botts et al., 1996; Burlakova et al., 2012).

Although zebra mussels invaded the lake first, the rapid expansion in the range of the invasive quagga mussel in Lake Constance has caused a dominance shift from zebra to quagga mussels. However, specific effects of live mussels could offer additive benefits for multiple taxa, including newly settled conspecific zebra and quagga mussels. In fact, taking the whole experimental setup into consideration, the entire abundance of newly settled quagga mussels was higher than the entire abundance of newly settled zebra mussels. These data coupled with the current field observation in Lake Constance indicate that this and/or other similar facilitative invasions do cause major changes in biodiversity and ecosystem. The higher mortality rate of adult zebra mussels than adult quagga mussels may also contribute to the dominance shift. We also observed, based on a laboratory study (own unpublished data), that the quagga mussel has a faster growth rate than the zebra mussel.

Our experiment was conducted in 2018, only 2 years after the initial record of quagga mussels in the lake. Nevertheless, as shown by prior studies in other invaded systems, such initial reactions of macroinvertebrate communities toward dreissenids might reduce after several years (Strayer and Malcom, 2006; Karatayev et al., 2015). At the current particular population levels (5,708 ind./m²), we may conclude that quagga mussels would significantly impact the structure of the benthic community in this system. Indeed, this initial establishment of the species at this littoral site could be considered surprising; however, the highest densities of invasive populations have been recorded in higher depths of North American lakes, such as 342,000, 75,000, and 16,000 ind./m² in Lakes Erie, Huron, and Michigan, respectively (Nalepa, 1995; Howell et al., 1996; Nalepa et al., 2009).

We were able to note significant differences in the food resources of newly settled macroinvertebrates using mainly δ¹³C values (Layman et al., 2012). Among the four treatments applied, the newly settled quagga mussels in B treatments showed higher δ¹³C values than in Z treatments, which was probably related to the smaller size of newly settled quagga mussels in B treatments (p < 0.001). The newly settled snail P. antipodarum is known as a grazer and feeds on periphyton (Larson and Ross Black, 2016), and as expected, showed higher δ¹³C values in Z treatments than in S treatments, thus the higher mortality of the adult zebra mussels in Z treatments may have generated a different nutrient resource. Comparing Q treatments and the other treatments, there were no significant differences in δ¹³C and δ¹⁵N values of newly settled macroinvertebrates. Previous studies by Gergs et al. (2011) found that the stable isotope values of gammarids were significantly and positively related to the biodeposition of zebra mussels. Here, the presence of the adult quagga mussels did not show a strong influence on the isotope values of the newly settled macroinvertebrates, which might be related to the smaller study area (experiment) and shorter study (invasion) period than the previous study.

Nonetheless, the newly settled quagga mussels had significantly higher δ¹³C values than the newly settled zebra mussels in each treatment. Quagga and zebra mussels attained from shallow depth of Lake Constance in 2018 also showed higher δ¹³C values of quagga mussels than that of zebra mussels (own unpublished data). The diet of dreissenids can be very flexible, and contains phytoplankton and suspended detritus (Garton et al., 2005). The benthic primary producers normally have a higher δ¹³C value compared to pelagic primary producers (more negative) (France, 1995; Mitchell et al., 1996). In our study period, the newly settled quagga mussels might consume a higher percentage of benthic food or a lower percentage of pelagic food, in comparison to newly settled zebra mussels. A similar result showed that after the invasion by smallmouth bass and rock bass, lake trout switched to consuming a lower percentage of littoral prey and had a lower δ¹³C value than the reference lakes (vander Zanden et al., 1999). One study found that the δ¹³C values of quagga mussels are more negative than those of zebra mussels (Verhofstad et al., 2013), and the explanation was that quagga mussels consumed more chemoautotrophs which have highly reduced δ¹³C values. The difference between our study and the study of Verhofstad et al. (2013) may be due to the differences in the substrates, depth ranges, or measurements. Some other studies did not find any difference in the stable isotope values between quagga and zebra mussels (Garton et al., 2005).

No overlaps in the isotopic niche of quagga and zebra mussels were observed in our field experiment. The lengths of newly settled quagga and zebra mussels in S, Q, and Z treatments did not show significant differences, thus the differences in stable isotope values of quagga and zebra mussels were not caused by the differences in the size of dreissenids. At our study site where quagga and zebra mussels coexisted, the dietary segregation or partition of colonized space may be the mechanism for the co-existence of quagga and zebra mussels (Diggins et al., 2004; Rothhaupt et al., 2014). Our experiment was conducted in 2018 in Lake Constance, only 2 years after the first record of quagga mussels. During the early colonization period, quagga mussels might occupy a different niche from zebra mussels to reduce resource competition with zebra mussels (Jackson and Britton, 2014).