© V. Gulin Beljak et al., Published by EDP Sciences 2023

1 Introduction

Ecological restoration refers to the process of managing and/or supporting the recovery of a degraded or destroyed ecosystem as a means of maintaining ecosystem resilience and conserving biodiversity (Convention on Biological Diversity, 2016). Restoration efforts (including river and stream restoration) are increasing daily worldwide (Palmer et al., 2014). The year 2021 was designated as the beginning of the United Nations (UN) Decade on Ecosystem Restoration (2021–2030) − a global movement initiated by the UN General Assembly to support ecosystem restoration and reverse global ecosystem degradation. In these frames, the European Union (EU) Biodiversity Strategy aims to restore 25,000 km of European rivers by 2030 (at least 20% of land and water areas in the EU), and by 2050 to cover all ecosystems in need of restoration (Publications Office of the European Union, 2022).

The European Mediterranean region has been identified as one of the most climate-vulnerable regions and a climate change "hotspot" (Salvia et al., 2021). The increasing intensity of droughts, driven by longer dry periods and shorter wet periods have made this region one of the most susceptible to soil degradation and desertification in Europe (Ferreira et al., 2020). Carbonate formations and karst features characteristic of this area are subjected to the global problem of climate change, as many karst springs, rivers, and streams face significant hydrological changes such as increased flow intermittence characterised by seasonal loss of flow and surface drying (Frollini et al., 2022; Patekar et al., 2021; Sivelle et al., 2021; Stubbington et al., 2018).

Croatia hosts many karst features, as almost half of the country is karst territory (Patekar et al., 2021). Most recognized karst features are the tufa barriers, which are formed by interactions between precipitation of carbonate minerals and resident organisms at specific ambient conditions (Golubić et al., 2008). This process is largely determined by geochemical conditions/environmental parameters (i.e., oversaturation of water with calcium carbonate, slightly alkaline pH, low concentration of dissolved organic matter) and is therefore very sensitive to environmental changes (Pentecost, 2005). Many tufa-depositing forms are highly affected by local anthropogenic interference with the environment, the rapid expansion of tourism, and the increasing threat of climate change (Šiljeg et al., 2020).

Although not many restoration activities have been carried out so far in tufa-depositing systems in Croatia, there have been some successfully finished attempts. In the tufa-depositing system of the Plitvice Lakes National Park, macrophyte overgrowth was identified as one of the main factors affecting the tufa-depositing process. Restoration methods that focused on macrophyte removal resulted in improved tufa deposition due to better aeration and reduced organic matter loading (Miliša et al., 2016;Pavlus and Novosel, 2004). In another Croatian national park (Krka National Park), parts of the Skradinski buk tufa barrier dried up due to uncontrolled growth of the invasive plant species Ailanthus altissima (Mill.) Swinge, whose strong root system caused overgrowth of the barrier, interruption of water flow, and eventual drying up of the streams (Gulin et al., 2021). Restoration efforts to remove the invasive plant from the Skradinski buk tufa barrier consisted of a comprehensive aerial survey (Phantom 4) and detailed vegetation mapping of a 1-ha experimental area at the barrier, which helped to locate the invasive plant specimens. Once the invasive plant individuals were identified, they were mechanically removed on several occasions during August 2017 with permission from the Croatian Agency for Environment and Nature. Within two months of removal, five streams in the experimental area, that had previously dried up, were reactivated (restored in terms of water connectivity).

The success of these restoration efforts was explored and confirmed by the results of Gulin et al. (2021, 2022), who used periphytic phagotrophic protists (e.g., ciliates) as bioindicators of restoration success. These results showed that the protist community structure and its functional metrics at restored tufa-depositing streams were responsive to the environmental conditions such as changing hydrology (i.e., occasional high runoff or drought), soil drainage, and extensive inorganic matter (i.e., tufa) deposition. The close connection between periphytic protists and the process of tufa deposition is an important element of the dynamics of the tufa-depositing ecosystems such as Krka National Park (Matoničkin and Pavletić, 1962; Primc-Habdija et al., 2001). The periphyton matrix retains and accumulates both organic matter particles and tufa particles, increasing microhabitat complexity/structural heterogeneity (Matoničkin Kepčija et al., 2011; Risse-Buhl et al., 2015; Dzubakova et al., 2018). Because periphyton is a substrate-bound community, it has shown promise in several studies for restoring lotic ecosystems. Timoner et al. (2012, 2014) have found changes in functional composition in response to disturbance and environmental changes such as desiccation, flooding, or eutrophication, and more recently Atristain et al. (2023) have shown that periphyton functional composition responds to restoration activities such as dam removal, opening up room for further implementation of periphyton in future studies. The usefulness of periphytic protists in detecting the ecological success of restoration efforts is due to the fact that they have several characteristics of reliable ecological indicators: (1) they are ubiquitously distributed throughout ecosystems, (2) they respond rapidly to environmental changes, (3) they are quantifiable at multiple levels of biological organization (species and community), and (4) have effects at levels above and below their position in the food web (Wu, 2016). Ciliates are commonly found in aquatic and terrestrial ecosystems (Esteban and Fenchel 2021; Singer et al., 2021) and, with over 8000 species described, play an essential role in microbial food webs as mediators of matter and energy transfer across trophic levels (Singer et al., 2021).

The results of the previous study by Gulin et al. (2022) showed that ciliate taxonomic metrics (abundance, species richness, Shannon and Simpson derived True diversity indices) were higher in restored vs. unrestored (i.e., control) streams for two different immersion periods (one- and two-months). On the other hand, functional metrics such as functional dispersion (FDis) in restored streams were lower during the one-month immersion period, but higher during the two-months period. Thus, we hypothesized that there should be differences in species and functional trait composition between the two (restored vs. unrestored/control) stream types. The first objective of the present study was (i) to examine the ciliate species occurrence in restored and control streams, and recommend species with indicator potential for restored streams. The second objective was (ii) to compare restored vs. control streams using community weighted means (CWM), to determine if there are differences in certain functional traits between the two stream types at the two different immersion periods (one- vs. two- months). We assumed that the lower functional diversity in restored streams as observed by Gulin et al. (2022) was a consequence of functional redundancy, i.e., that restored streams were dominated by species occupying the same functional niche. Furthermore, we hypothesized that ciliate communities of restored streams would consist primarily of bacterivorous motile species readily mobilized during occasional high flow or drought, and that these species are also characteristic representatives of terrestrial ecosystems as restored streams drain the local forest floor. Our third objective was (iii) to investigate the effects of environmental parameters (i.e., water physicochemical parameters, organic and inorganic (i.e., tufa) matter content, chlorophyll a concentration) on ciliate species composition and functional traits in restored streams. We assumed that organic matter content and nitrite concentration would affect ciliate occurrence and functional trait composition the most, as these parameters were identified as the most important for community-level metrics in the previous study (Gulin et al., 2022).

The results of the present study will help to understand how ciliates respond to stream restoration at the species and functional trait levels. Furthermore, these results could serve as guidance in identifying species that can serve as indicators of stream restoration processes, not only in Croatian tufa-depositing systems, but also at a larger spatial scale (e.g., regional/Mediterranean).

2 Materials and methods

2.1 Study area

The Krka River is a karst river in the Dinaric region of Croatia, located in the central part of the north-eastern Adriatic coast (Fig. 1). Its headwaters are located near the Dinara Mountain and consist of several independent springs from which the river enters a canyon and is characterized by several lotic and lentic (e.g., Lake Visovac) areas intersected by tufa barriers (e.g., Roški slap, Skradinski buk) before it flows into the Adriatic Sea near the town of Šibenik (Fig. 1). With a total catchment area of about 2427 km² (Bonacci and Ljubenkov, 2005) and an average annual discharge of 47.4 m³ s⁻¹ (1990–2009), the Krka River is considered a medium-sized river in Croatia (Čanjevac and Orešić, 2015).

This study was conducted at the longest and last tufa barrier (Skradinski buk) in the Krka River watercourse (Fig. 1), part of which has been protected by the National Park category since 1985. A 1 ha trial area was delineated, from which the invasive plant Ailanthus altissima (Mill.) Swinge (tree of heaven) was mechanically removed in August 2017 with the permission of the Croatian Agency for Environment and Nature, after being classified as dominant by aerial survey and vegetation mapping. Within two months of removal, five streams in the experimental area, that had previously dried up, were completely reactivated. The sampling design included seven sampling sites: two located in unrestored, i.e., control streams (C sites) where water was present before and after plant removal, with the C1 site having typical, well-developed moss cover and the C2 site having no moss cover; and five located in restored streams, i.e., newly reactivated streams (N sites) selected after the removal of the invasive plant species. Growth of mosses on surfaces can provide substrates for calcite nucleation and trap detrital calcite, accelerating tufa deposition; thus, moss-covered substrates represent typical microhabitats of tufa barriers in Krka National Park (Matoničkin and Pavletić, 1962; Primc-Habdija et al., 2001). Aquatic mosses in lotic ecosystems are also a suitable substrate for colonization of numerous microscopic and macroscopic metazoans, as they provide protection from water flow and wash-off from the habitat and have a large capacity for retention and accumulation of organic matter, which is a food source for benthic organisms (Traunspurger, 2002; Dražina et al., 2013, 2014). Details of the sampling sites and their descriptions are provided in Table 1.

Fig. 1 Location of Croatia within the Mediterranean region (A) and the Krka River catchment with the location of Skradinski buk barrier (B).

2.2 Sampling design and ciliate identification

Due to the high seasonality of the Krka River (Schöll et al., 2012), the research included four sampling campaigns covering four different seasons, in period between May 2019 and June 2020. At the beginning of every sampling event, at each sampling site, two 3D-printed Plexiglass carriers with three glass slides (7.6 × 2.6 cm) in each one (Gulin et al., 2022), were immersed in the water in the middle section of each stream, one for a period of one month, and second for a period of two months (Tab. 2). Since glass slides were partially covered by Plexiglass, the total surface area of each slide available for periphyton development was 17.18 cm². Prior to immersion into the stream, the slides were cleaned with detergent, 1 M hydrochloric acid and distilled water. During each sampling event (after one-month and two-months immersion period in each season) three slides were collected per sampling site. In total, 168 slides were placed into the streams, 84 for the one-month immersion period and another 84 for the two-months immersion period. However, only 55 slides were collected for the one-month immersion period and 52 for the two-months immersion period, due to seasonal drying of streams and disturbance by visitors of the National Park (sampling sites were located in the close vicinity of visitor trails and probably aroused curiosity). Collected slides were placed in the plastic containers filled with a small amount of ambient water and stored at 4 °C in the dark. The slides were examined within a maximum of 48 h from collection, using Zeiss Axioimager A2 with DIC objectives and Axiocam 305 digital camera. Ciliates were identified at species level using Zen 2.4 imaging software and relevant literature (Foissner et al., 1991, 1992, 1994, 1995; Foissner and Berger, 1996). For each taxon, 10-15 photomicrographs were taken and subjected to morphometric measurement. Additional video clips were used to record movements and distinguishing features.

2.3 Environmental parameters

On each sampling event, at all sampling sites (C1, C2, N1-N5), in situ measurements of the following environmental parameters were made using the appropriate portable field meters: temperature (T) and dissolved oxygen concentration (DO) (oximeter OXI 96, WTW GmbH, Weilheim, Germany), pH (pH meter 330i, WTW GmbH, Weilheim, Germany), conductivity (Cond) (conductometer Sension 5, Hach, Loveland, Colorado, USA), and flow velocity (FV) (flow velocity meter P600, Dostmann electronic GmbH, Wertheim- Reicholzheim, Germany). An additional water sample (1 L) was collected and transported at 4 °C for subsequent laboratory analysis of the following parameters: alkalinity (Alk), total water hardness (TWH), and concentrations of nitrite (N-NO₂⁻), nitrate (N-NO₃⁻), and orthophosphate (P-PO₄³⁻) (APHA, 1985), and total chemical oxygen demand (COD) using the standardized acidic potassium permanganate titrimetric method (Deutsches Institut für Normung, 1986).

After microscopic identification of ciliates, periphyton sample from each slide was separated into two equal parts. One part was used to determine the content of organic and inorganic matter in the sample. For this purpose, the samples were dried at 104 °C to constant weight, then ashed at 400 °C for 4 h and reweighed. The mass difference between the dried and the ashed sample represents the organic matter (OM) content, while the mass difference between the ashed sample and the glass slide represents the content of inorganic matter (IM), i.e., deposited tufa. The values were expressed as mass (mg) of organic or inorganic matter content per cm² of surface area. The other part of the sample was used to determine chlorophyll a concentration by ethanol extraction method (Nusch, 1980). Values were expressed as mass (µg) of chlorophyll a per cm² surface area.

2.4 Data preparation

The functional categorization of ciliates used for calculation of community weighted mean (CWM) values was conducted according to Gulin et al. (2022). Calculating CWM values is a useful way to calculate community trait values weighted by abundance of species in that community (Ricotta and Moretti, 2011). Ciliates were assigned to categories for the following functional traits: food source (“algae”; “algae, bacteria”; “algae, diatoms”; “diatoms”; “bacteria”; “cyanobacteria”; “omnivorous”; “phagotrophic protists, small metazoans”), feeding strategy (“filtration”; “predation”), ecosystem preference (“active sludge”; “lentic (freshwater)”; “lentic and lotic (freshwater)”; “soil”), habitat preference (“benthos”; “benthos, periphyton”; “periphyton”; “planktonic”), motility (“motile”; “semi sessile”;“sessile”), mode of locomotion (“crawling”; “free-swimming”; “gliding”; “jumping”, “rotating”; “unknown”) and life form (“colonial”; “solitary”). The detailed table of allocated functional traits can be found in Gulin et al. (2022), and it has been used as a convenient reference in other recent studies concerning ciliates in the Krka River (e.g., Gulin et al., 2022; Gulin Beljak et al., 2022).

2.5 Data analysis

2.5.1 Indicator species analysis

Indicator species analysis (IndVal), which combines both abundance and frequency of a species's occurrence, was used to determine potential indicator species for C and N sites. IndVal values were calculated in R v. 4.2.2. (R Core Team, 2022) using the package indicspecies v. 1.7.12. Correlations between indicator species and environmental parameters were calculated with Spearman's correlations using the “rcorr” function from the R Hmisc package v. 4.7-2 and visualized with the “corrplot” function from the R corrplot package v. 0.92.

2.5.2 Functional composition of periphytic ciliates

The CWM values were calculated for each trait category at each site to quantify shifts in mean functional trait values within communities resulting from environmental selection for certain functional trait categories (Ricotta and Moretti, 2011). The CWMs were calculated using the “functcomp” function from the FD package v. 1.0-12.1. (Laliberté and Legendre, 2010; Laliberté et al., 2014). To determine if there were significant differences in CWM values for each trait category between C and N sites for both immersion periods, generalized linear mixed models (GLMMs) were constructed using SPSS Statistics v. 28.0. (IBM Corp, 2021). The variable “site” was considered a fixed effect, while the variables “replicate” and “season” were included as random effects in both model variations, as recommended by Jost (2007). First-order autoregressive (AR1) covariance structure was assumed in all models, since data were collected repeatedly over time (in the course of several months) (Field, 2009). The differences in CWM values between C and N sites were visualized using R ggplot2 package v. 3.4.0. (Wickham, 2016).

2.5.3 Effects of environmental parameters on species and functional trait composition

The differences in the values of environmental parameters between C and N sites were tested using the analysis of similarity (ANOSIM) and have already been reported in Gulin et al. (2022). Correlations between CWM values that differed significantly between C and N sites and environmental parameters, were calculated with Spearman's correlations using the “rcorr” function from the R Hmisc package v. 4.7-2 and visualized with the “corrplot” function from the R corrplot package v. 0.92.

Redundancy analysis (RDA) in R was used to assess the influence of environmental parameters on the spatial distribution of ciliate species (only those represented by more than 5% in total abundance per site) and CWM values (only those that differed significantly between C and N sites) at C and N sites. RDA is designed for correlated response variables (which is often the case for CWMs) and therefore allows many response variables to be analyzed simultaneously (by forming RDA axes, which are linear combinations of response variables) (Kleyer et al., 2012; Šmilauer and Lepš, 2014). Values for ciliate abundance and CWMs were transformed by the Hellinger's transformation prior to species-based RDA and CWM-based RDA (Legendre and Gallagher, 2001). Multicollinearity between environmental parameters was checked using the “vif.cca()” function from the R vegan package v. 2.6-4, and if the value was greater than 2, multicollinearity was considered high. The function “forward.sel()” from the R adespatial package v. 0.3-20. was used to select the best explanatory variables. Monte Carlo test with 999 unrestricted permutations was used to determine the significance of the model, axis, and explanatory variables (environmental parameters) at the P < 0.05 level.

3 Results

3.1 Species occurrence and potential indicator species

Based on 1837 recorded ciliate individuals, 78 species were identified. The most abundant species at C sites during the one-month immersion period was Pleurotricha grandis (14.49% of the total abundance at C sites). Those with more than 10% of total abundance at C sites were Trithigmostoma cucullulus (12.08%,) and Cyclidium sp. (11.11%) followed by Holophrya ovum (9.18%), Vorticella aquadulcis-complex (8.21%), Aspidisca cicada (7.73%) and Pseudochilodonopsis algivora (7.25%). Other species were represented with less than 5%.

The most abundant species at N sites during the one-month immersion period was Platyophrya vorax, accounting for 72.58% of the total ciliate abundance at N sites, while additionaly only Cinetochilum margaritaceum had a proportion greater than 10% (10.27%). Chilodonella uncinata had a share of more than 5% in the total abundance (5.62%), while the other species were represented with less than 5%.

During the two-months immersion period, the species with highest proportion in total ciliate abundance at C sites were Holophrya discolor (35.71%) and Holophrya sp. (10.71%). The remaining species were represented by less than 5% in the total abundance. At N sites, C. uncinata (25.12%), Holophrya sp. (21.40%), P. vorax (16.98%) and Chlamydonella alpestris (6.28%) had the highest proportion in the total ciliate abundance, respectively. The remaining species were represented with less than 5%.

Indicator species analysis (IndVal) for the one-month period showed that P. vorax was the only species with indicator potential for N sites (stat = 0.275, p = 0.034). For C sites, the following species were detected as potential indicator species: V. aquadulcis-complex (stat = 0.318, p = 0.003), H. ovum (stat = 0.269, p = 0.011), and A. cicada (stat = 0.188, p = 0.037).

Indicator species analysis for the two-months period did not reveal any indicator species neither for N sites nor for C sites.

3.2 Functional trait composition

3.2.1 One-month immersion period

GLMMs for the one-month immersion period showed statistically significant differences between C and N sites for CWM values related to food source, ecosystem preference, habitat preference, and mode of locomotion, while CWM values related to motility, feeding strategy, and life form did not differ significantly between C and N sites. CWM values associated to feeding strategy, motility and life form can be found in Supplementary Figure 1.

The CWM values of food source at N sites differed significantly from those at C sites only for the functional categories “bacteria”, and “phagotrophic protists, small metazoans” (Tab. 3). N sites had significantly higher CWM values for the category “bacteria” and significantly lower values for the category “phagotrophic protists, small metazoans” than C sites (Tab. 3, Fig. 2).

The CWM values of ecosystem preference differed significantly between C and N sites for the functional categories of “lentic (freshwater)”, “lentic and lotic (freshwater)”, and “soil”. In comparison to C sites, N sites had significantly higher CWM values for the category “soil”, but significantly lower values for the categories “lentic (freshwater)” and “lentic and lotic (freshwater)” (Tab. 3, Fig. 3).

The CWM values of habitat preference differed significantly between C and N sites for the functional categories of “benthos, periphyton”, “benthos” and “planktonic”. The CWM value of category “benthos, periphyton” was significantly lower at N sites compared to C sites, whereas the opposite was found for the categories “benthos” and “planktonic” (Tab. 3, Fig. 4).

The CWM values of mode of locomotion differed significantly between C and N sites only for the functional categories “crawling” and “free swimming”. The CWM value of category “crawling” was significantly lower at N sites than at C sites, while the opposite was found for the category “free swimming” (Tab. 3, Fig. 5).

Fig. 2 Differences in community-weighted mean (CWM) values (± standard error) of food source functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.

Fig. 3 Differences in community weighted mean (CWM) values (± standard error) of ecosystem preference functional trait categories in ciliate communities between control (C) and restored streams (N) during the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.

Fig. 4 Differences in community weighted mean (CWM) values (± standard error) of habitat preference functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.

Fig. 5 Differences in community weighted mean (CWM) values (± standard error) of mode of locomotion functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.

3.2.2 Two-months immersion period

For the two-months immersion period, GLMMs showed statistically significant differences between C and N sites for CWM values of ecosystem and habitat preference, while the other CWM values did not differ significantly between restored and control streams. CWM values of food source, feeding strategy, motility, mode of locomotion and life form can be found in Supplementary Figure 2.

The CWM values of ecosystem preference differed significantly between C and N sites only for the functional category "soil", whereby higher value was found at N than at C sites (Tab. 3, Fig. 6).

The CWM values of habitat preference differed significantly between C and N sites only for the functional categories of “benthos, periphyton” and “benthos”. The CWM values of category “benthos, periphyton” was significantly higher at N than at C sites, whereas the opposite was found for the category “benthos” (Tab. 3, Fig. 7).

Fig. 6 Differences in community weighted mean (CWM) values (± standard error) of ecosystem preference functional trait categories in ciliate communities between control (C) and restored streams (N) during the two-months immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.

Fig. 7 Differences in community weighted mean (CWM) values (± standard error) of habitat preference functional trait categories in ciliate communities between control (C) and restored streams (N) for the two-months immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.

3.3 Environmental parameters

For the one-month immersion period, N sites had significantly lower COD values but significantly higher nitrite values than C sites. For the two-months immersion period, N sites had significantly lower COD, pH and orthophosphate values, while conductivity values were higher compared to C sites. A detailed overview of all environmental parameters has been reported in Gulin et al. (2022) and can be found in Supplementary Figure 3.

Organic matter content was lower at the N sites than at C sites in both periods, but with occasional extremely high values, while there were no such extreme outliers at C sites. Inorganic matter content was lower at N sites than at the C sites during the one-month period, and the reverse was true during the two-months period. Chlorophyll a concentration was higher at N sites compared to C sites during the one-month period and lower during the two-months period, but with extremely high chlorophyll a concentrations during both periods in contrast to N sites. A detailed overview of organic, inorganic matter content and chlorophyll a concentration has been reported in Gulin et al. (2022) and can be found in Supplementary Figure 4.

3.4 Effects of environmental parameters on species and functional trait composition

Spearman's correlation coefficient between ciliate species with indicator potential and environmental parameters was performed only for the one-month immersion period because IndVal analysis did not identify any indicator species for the two-months immersion period. P. vorax was found to be significantly positively correlated with nitrite concentration (R = 0.409, p = 0.013). Detailed Spearman's correlation results for the one-month immersion period can be found in Supplementary Figure 5.

RDA ordination for the one-month immersion period was statistically significant (F = 2.856, p = 0.001). The first constrained axis (RDA1) described 13.6% of the variance (F = 6.288, p = 0.001), the second (RDA2) explained 11% of the variance (F = 5.069, p = 0.003), while the third axis (RDA3) explained 8.5% (F = 3.897, p = 0.028) and was therefore not visualized (Fig. 8). Three explanatory variables (environmental parameters) explained 33.1% of the variance in the data: nitrites (F = 5.041, p = 0.001), COD (F = 2.140, p = 0.030), inorganic matter (F = 2.358, p = 0.032). RDA showed that species representative for N sites, such as P. vorax and C. margaritaceum, were highly correlated with nitrites and organic matter. C sites exhibited stronger grouping compared to N sites. RDA ordination for the two-months immersion period was not statistically significant.

Spearman's correlation coefficient between CWM values and environmental parameters was calculated for both immersion periods. During the one-month immersion period, CWM food source values for the category “bacteria” were significantly negatively correlated with inorganic matter (R = –0.333, p = 0.013), while CWM mode of locomotion values for the category “crawling” were significantly negatively correlated with COD (R = –0.279, p = 0.039). The detailed Spearman correlation results for the one-month period can be found in Supplementary Figure 6.

As for the two-month immersion period, the CWM ecosystem preference values for the category “soil” were found to be significantly negatively correlated with COD (R = –0.447, p = 0.001) and organic matter (R = –0.303, p = 0.029), while they were significantly positively correlated with pH (R = 0.314, p = 0.024). CWM habitat preference values for the category “benthos, periphyton” were significantly negatively correlated with pH (R = – 0.341, p = 0.013), while the opposite was found for the category “benthos” (R = 0.381, p = 0.005). A significant negative correlation between the category “benthos” and COD was also found (R = –0.276, p = 0.048). The detailed Spearman's correlation results for the two-months period can be found in Supplementary Figure 7.

Only significantly different CWM values between C and N sites were used for RDA ordination for both immersion periods. RDA ordination of CWM values and environmental parameters for the one-month immersion period was found to be statistically significant (F = 2.067, p = 0.001). The first constrained axis (RDA1) described 11.4% of the variance (F = 6.891, p = 0.001).

It was found that the first unconstrained axis (PC1) represented 23.4% of the total variance. None of the explanatory variables (environmental parameters) had a significant impact on the CWM values (Fig. 9).

The RDA ordination of CWM values and environmental parameters for the two-months immersion period proved to be statistically significant (F = 3.521, p = 0.006). The first constrained axis (RDA1) described 20.8% of the variance (F = 12.979, p = 0.004) while the first unconstrained axis (PC1) represented 56.8% of the total variance. Significant correlations were found between CWM habitat preference values for the category "benthos, periphyton" and COD (F = 3.776, p = 0.041) (Fig. 10).

Fig. 8 Redundancy analysis (RDA) showing the relationships between ciliate species (red arrows) and environmental parameters (blue arrows) among control (C) and restored (N) sites for the one-month immersion period. Abbrevations: N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; AC = Aspidisca cicada; CU = Chilodonella uncinata; CM = Cinetochilum margaritaceum; CYC = Cyclidium sp.; HO = Holophrya ovum; PV = Platyophrya vorax; PG = Pleurotricha grandis; PA = Pseudochilodonopsis algivora; TC = Trithigmostoma cucullulus; VAC = Vorticella aquadulcis-complex.

Fig. 9 Redundancy analysis (RDA) showing the relationships between CWM values (red arrows) and environmental parameters (blue arrows) among control (C) and restored (N) sites for the one-month immersion period. Abbrevations: N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; FSB = food source “bacteria”; FSPPSM = food source “phagotrophic protists, small metazoans”; EPLF = ecosystem preference “lentic (freshwater)”; EPLLF = ecosystem preference “lentic and lotic (freshwater)”; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”; HPP = habitat preference “planktonic”; MLFS = mode of locomotion “free-swimming”; MLC = mode of locomotion “crawling”.

Fig. 10 Redundancy analysis (RDA) showing the relationships between CWM values (red arrows) and environmental parameters (blue arrows) among control (C) and restored (N) sites for the two-months immersion period. Abbrevations: COD = chemical oxygen demand; OM = organic matter; CHL = chlorophyll a concentration; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”.

4 Discussion

4.1 Species with indicator potential for restored tufa-depositing streams

Previous study by Gulin et al. (2022) showed that several taxonomic metrics of ciliates were higher in restored streams than in control streams at the Skradinski buk tufa barrier for both immersion periods, while functional metrics such as functional dispersion (FDis) were lower during the one-month immersion period but higher during the two-months period. This study investigated further the possible reasons for these results by examining the differences in ciliate species and functional trait composition between restored and control streams; highlighting species that are characteristic and distinctive of restored and control streams (i.e. species with indicator potential); and understanding the environmental parameters that influence the composition of ciliate species and functional traits.

The results of this study revealed that although restored streams had higher species richness than the control streams during both immersion periods, their communities were strongly dominated by few species in terms of relative abundance. Contrary to that pattern, control streams had comparable share of 6 to 7 species. This was more pronounced during the one-month immersion period, when the colpodid P. vorax accounted for nearly 73% of the total ciliate abundance in restored streams, while the other species such as C. margaritaceum and C. uncinata accounted for only 5–10%. P. vorax was also identified as the only one with indicator potential for restored tufa-depositing streams, suggesting that it was not only abundant but also common in the periphyton of restored streams during the entire study period, which included repeated sampling. However, during the two-months immersion period, the community structure of the periphyton of restored streams showed less pronounced dominance of certain species, and proportions of species with the highest relative abundance, such as P. vorax, C. uncinata, Holophrya sp. and C. alpestris, were more or less equal.

4.2 Differences in ciliate functional traits between control and restored streams − restored tufa-depositing streams as a dynamic interface between terrestrial and aquatic ecosystems

The strong dominance of a few species in terms of relative abundance at restored streams was reflected in functional diversity during the one-month immersion period most likely because these species performed the same functional roles, resulting in higher competition between species and within species (Hooper et al., 2005). Restored streams favoured better competitors such as P. vorax. Species from the genus Platyophrya have several adaptations to highly competitive habitats such as being free-swimming (vs. periphyton crawling) and having bacteriophagous microstomes (Foissner and Wolf, 2009). This led to lower functional diversity at restored streams compared to control streams and was also evident regarding the species feeding on phagotrophic protists and small metazoans, which were significantly less abundant in restored streams than in the control streams as they are selective feeders. On the other hand, a longer immersion period, i.e. two-months immersion period, allowed the separation of functional niches, i.e., niche partitioning (Loreau and Hector, 2001) in restored streams, which facilitated species coexistence by reducing competition (DeLong and Vasseur, 2012). These results highlight the importance of study design in periphyton-based studies, as sampling periphyton that has only been immersed in water for a short period of time can lead to potentially incorrect conclusions about stream restoration outcomes and success evaluation.

The most abundant ciliate species of restored streams indicate a close link between terrestrial and aquatic ecosystems, both in terms of taxonomic and functional trait composition: P. vorax, C. margaritaceum, C. uncinata, and C. alpestris are frequently found not only in aquatic (freshwater) ecosystems, but also in soils and nutrient-rich man-made systems such as activated sludge (Foissner, 2016) in contrast to usually oligotrophic tufa-depositing systems. Restored streams in our study drained the forest floor that was previously completely dry, and it appears that this was the predominant factor affecting species and functional trait composition of ciliates in restored tufa depositing streams, suggesting that soil drainage in restored streams highly affects periphyton structure and ecosystem function.

Soil ciliate species are generally bacterivores, as was the case with the most species characteristic for restored streams in our study (C. margaritaceum, C. uncinata and C. alpestris), with the exception of P. vorax, which is an omnivorous species meaning that it feeds on other food sources additionally to bacteria (Esteban and Fenchel, 2021). By grazing, ciliates control the size and composition of bacterial communities, indirectly regulating decomposition processes (Jia et al., 2021). The functional category of soil species showed a significant positive correlation with pH and a negative correlation with COD and organic matter in our study. Because restored streams often had higher nitrite concentrations than control streams due to increased drainage favouring reduced organic forms (e.g., nitrites), their concentration may have influenced the organic matter decomposition process and bacterial communities, indirectly affecting the abundance of soil ciliate species. Dissolved forms of nitrogen such as nitrites have the potential to affect organic matter decomposition depending on their concentration and pH (Riggs and Hobbie, 2016). Higher concentrations of dissolved forms of nitrogen can promote microbial growth and thus the decomposition process (Riggs and Hobbie, 2016), resulting in more food for ciliates, as found in the present study where the abundance of P. vorax was significantly positively correlated with nitrite concentration. However, when the concentration of dissolved nitrogen forms exceeds a threshold, it can lead to a reduction in microbial biomass (Jing et al., 2021) and consequently to a reduction in food for most soil ciliate species. This may have been the case in our study, but more data on carbon and nitrogen fluxes in restored streams are needed to better understand and verify these conclusions. P. vorax has recently been associated with the degradation of refractory carbon sources in terrestrial ecosystems (Pastorelli et al., 2022) and may therefore play an important role in carbon fluxes of restored streams as they seem to represent an interface between aquatic and terrestrial ecosystems.

Unlike control streams, restored streams in our study supported more motile, free-swimming ciliate species. These findings are consistent with the observed trends of Risse-Buhl and Kusel (2009), who showed that periphyton in streams with high flow rates (in our case restored streams) is dominated by highly motile, flattened, gulper feeding ciliates such as P. vorax, C. uncinata, and C. alpestris in our study. Occasional events of high flow in restored streams reported by Gulin et al. (2022) also seems to have benefited species characterized as planktonic such as C. margaritaceum but more during the one-month immersion period than the two-months immersion period. Planktonic species probably appeared with seston from the upstream Visovac Lake, which can be a source of planktonic organisms along lotic stretches in the study area as it has already been reported in Gulin et al. (2021). Risse-Buhl and Küsel,(2009) also showed that flattened, filter-feeding ciliates dominate at slow flow rates, as was the case for species found in control streams, such as crawling A. cicada and semi-sessile V. aquadulcis-complex. Control sites favoured crawling species commonly found in stream benthos and periphyton of lentic and lotic parts along the Krka River (Kulaš et al., 2021; Gulin et al., 2021; Gulin Beljak et al., 2022).

Besides of flow velocity, the predominance of motile ciliates in restored streams could be a consequence of microhabitat complexity, a parameter we did not quantify in our study and which could be the missing explanatory variable in our RDA (evident from the higher representation of unconstrained axis in total variance compared to constrained). The combination of soil particles drained from the local soil, deposited tufa and retention of organic matter can increase periphyton heterogeneity in restored streams and provide a greater range of resource niches. These niches might be better exploited by motile, free- swimming than crawling ciliates, as they are likely better able to avoid being buried by the intense tufa deposition. The urostilids recorded in our study, such as the genera Oxytricha, Tachysoma, and Urosomoida, seem to support this idea. These ciliates possess morphological adaptations that allow them to hide or hunt among or within particle aggregates − most are short to long elliptical in outline, ventrally flattened, and dorsally distinctly vaulted (Berger, 1999). Although their relative abundance was not comparable to that of other species in the community and therefore they were not recognized as potential indicators of restored streams, they were more abundant in restored streams than in control streams (Gulin et al., 2022).

Heterogeneity, i.e., habitat complexity has often been considered in successful stream restoration efforts (Baattrup-Pedersen et al., 2022), but usually at the reach level (i.e., channel morphological features) rather than at the microhabitat level. Higher heterogeneity achieved through stream restoration has been associated with shifts in the functional composition of macroinvertebrate communities (Frainer et al., 2017; Hasselquist et al., 2018), but periphyton-level responses have been little studied (Huang et al., 2018) and require further attention. In their most recent comprehensive literature review, Baattrup-Pedersen et al. (2022) call for more studies that use functional parameters to evaluate stream restoration efforts of all ecosystem types at different spatial and temporal scales to allow for more robust comparisons between natural and restored ecosystems. Phagotrophic protists such as ciliates provide excellent model systems for studying some fundamental effects of ecological interactions on functional traits, that may also apply to larger organisms (Montagnes et al., 2012), so the results of this study may reveal potential responses of macroinvertebrates or fishes in the long term. Given the ubiquity of most ciliate species, as well as community structure in similar habitats (Fenchel and Finlay, 2004), a comparable functional response to restoration might be expected in freshwater ecosystems with high levels of sedimentation or in other karstic waters that have close links to terrestrial ecosystems.

5 Conclusions

Our results emphasize the potential value of functional traits as a tool for assessing the success of ecological restoration at the microscale. Future evaluation of the success of stream restoration in tufa-depositing environments should incorporate data on traditional taxonomic metrics, but also data on functional traits of substrate-related biota (periphyton, benthos) and quantifiable measures of microhabitat heterogeneity. In addition, it would be useful to include measures of organic matter decomposition (carbon and nitrogen fluxes) to enable quantification of microbially mediated decomposition for a comprehensive assessment of ecological responses to stream restoration.

Acknowledgments

We would like to thank Gordana Goreta who supervised the project and Nikolina Smolić who facilitated with the selection of sampling sites from Krka National Park; Antun Beljak for his assistance in the field; Sunčica Bosak for loaning microscope equipment and Ivan Tekić for his kind help in preparing the map of the study area.

Supplementary Material

The Supplementary Material is available at https://www.kmae-journal.org/10.1051/kmae/2023012/olm.

Fig. S1: Differences in community-weighted mean (CWM) values (± standard error) of feeding strategy; motility; and life form functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period.

Fig. S2: Differences in community weighted mean (CWM) values (± standard error) of food source; feeding strategy; motility; mode of locomotion, life form functional trait categories in ciliate communities between control (C) and restored streams (N) during the two-month immersion period.

Fig. S3: Box plots showing selected enivironmental parameteres at control (C) and restored sites (N) for the one -month and two-months immersion period, Asterisk symbol (*) indicates statistically signification differences amnog estimated means (ANOSIM p<0.05), Upper and lower edges of the boxes are the first and third quartiles; thi line inside the box represents the median; individual dots are outliers, Abbreviation; COD- total chemical oxygen demand, N-NO₂ -nitries, P-PO³₄ -orthophosphates.

Fig. S4: Box plots showing organic matters and inorganic matter content; and chlorophyll σ concentrations at control (C) and restored sites (N) for the one -month and two-months immersion period. Asterisk symbol (*) indicates statistically signification differences amnog estimated means (ANOSIM p<0.05), Upper and lower edges of the boxes are the first and third quartiles; thi line inside the box represents the median; individual dots are outliers,

Fig. S5: Pairwise Spearman's rank correlation matrix of ciliate taxa abundance and environmental parameters for the one-month immersion period. Only taxa identified by Indicator species analysis (IndVal) were used for the purpose of correlation. Circle size reflects the magnitude of the correlation coefficient. Color represents the level of Spearman's correlations (blue means positive correlation and red means negative correlation). X represents statistically insignificant correlations (p>0.05). Abbrevations: AC- Aspidisca cicada, HO = Holophrya ovum, PV = Platyophrya vorax, VAC = Vorticella aquadulcis-complex; N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration.

Fig. S6: Pairwise Spearman's rank correlation matrix of community weighted mean (CWM) values of functional trait categories and environmental parameters for the one-month immersion period. Only categories that were significantly different according to GLMMs were used for correlation. Circle size reflects the magnitude of the correlation coefficient. Color represents the level of Spearman's correlations (blue means positive correlation and red means negative correlation). X represents statistically insignificant correlations (p>0.05). Abbrevations: N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; FSB = food source “bacteria”; FSPPSM = food source “phagotrophic protists, small metazoans”; EPLF = ecosystem preference “lentic (freshwater)”; EPLLF = ecosystem preference “lentic and lotic (freshwater)”; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”; HPP = habitat preference “planktonic”; MLFS = mode of locomotion “free swimming”; MLC = mode of locomotion “crawling”.

Fig. S7: Pairwise Spearman's rank correlation matrix of community weighted mean (CWM) values of functional trait categories and environmental parameters for the two-months immersion period. Only categories that were significantly different according to GLMMs were used for correlation. Circle size reflects the magnitude of the correlation coefficient. Color represents the level of Spearman's correlations (blue means positive correlation and red means negative correlation). X represents statistically insignificant correlations (p>0.05). Abbrevations: COD = chemical oxygen demand; COND= conductivity; ORTH= orthophosphates; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”.

References

Supplementary Material

The Supplementary Material is available at https://www.kmae-journal.org/10.1051/kmae/2023012/olm.

Fig. S1: Differences in community-weighted mean (CWM) values (± standard error) of feeding strategy; motility; and life form functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period.

Fig. S2: Differences in community weighted mean (CWM) values (± standard error) of food source; feeding strategy; motility; mode of locomotion, life form functional trait categories in ciliate communities between control (C) and restored streams (N) during the two-month immersion period.

Fig. S3: Box plots showing selected enivironmental parameteres at control (C) and restored sites (N) for the one -month and two-months immersion period, Asterisk symbol (*) indicates statistically signification differences amnog estimated means (ANOSIM p<0.05), Upper and lower edges of the boxes are the first and third quartiles; thi line inside the box represents the median; individual dots are outliers, Abbreviation; COD- total chemical oxygen demand, N-NO₂ -nitries, P-PO³₄ -orthophosphates.

Fig. S4: Box plots showing organic matters and inorganic matter content; and chlorophyll σ concentrations at control (C) and restored sites (N) for the one -month and two-months immersion period. Asterisk symbol (*) indicates statistically signification differences amnog estimated means (ANOSIM p<0.05), Upper and lower edges of the boxes are the first and third quartiles; thi line inside the box represents the median; individual dots are outliers,

Fig. S5: Pairwise Spearman's rank correlation matrix of ciliate taxa abundance and environmental parameters for the one-month immersion period. Only taxa identified by Indicator species analysis (IndVal) were used for the purpose of correlation. Circle size reflects the magnitude of the correlation coefficient. Color represents the level of Spearman's correlations (blue means positive correlation and red means negative correlation). X represents statistically insignificant correlations (p>0.05). Abbrevations: AC- Aspidisca cicada, HO = Holophrya ovum, PV = Platyophrya vorax, VAC = Vorticella aquadulcis-complex; N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration.

Fig. S6: Pairwise Spearman's rank correlation matrix of community weighted mean (CWM) values of functional trait categories and environmental parameters for the one-month immersion period. Only categories that were significantly different according to GLMMs were used for correlation. Circle size reflects the magnitude of the correlation coefficient. Color represents the level of Spearman's correlations (blue means positive correlation and red means negative correlation). X represents statistically insignificant correlations (p>0.05). Abbrevations: N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; FSB = food source “bacteria”; FSPPSM = food source “phagotrophic protists, small metazoans”; EPLF = ecosystem preference “lentic (freshwater)”; EPLLF = ecosystem preference “lentic and lotic (freshwater)”; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”; HPP = habitat preference “planktonic”; MLFS = mode of locomotion “free swimming”; MLC = mode of locomotion “crawling”.

Fig. S7: Pairwise Spearman's rank correlation matrix of community weighted mean (CWM) values of functional trait categories and environmental parameters for the two-months immersion period. Only categories that were significantly different according to GLMMs were used for correlation. Circle size reflects the magnitude of the correlation coefficient. Color represents the level of Spearman's correlations (blue means positive correlation and red means negative correlation). X represents statistically insignificant correlations (p>0.05). Abbrevations: COD = chemical oxygen demand; COND= conductivity; ORTH= orthophosphates; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”.

All Tables

All Figures

	Fig. 1 Location of Croatia within the Mediterranean region (A) and the Krka River catchment with the location of Skradinski buk barrier (B).
In the text

	Fig. 2 Differences in community-weighted mean (CWM) values (± standard error) of food source functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.
In the text

	Fig. 3 Differences in community weighted mean (CWM) values (± standard error) of ecosystem preference functional trait categories in ciliate communities between control (C) and restored streams (N) during the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.
In the text

	Fig. 4 Differences in community weighted mean (CWM) values (± standard error) of habitat preference functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.
In the text

	Fig. 5 Differences in community weighted mean (CWM) values (± standard error) of mode of locomotion functional trait categories in ciliate communities between control (C) and restored streams (N) for the one-month immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.
In the text

	Fig. 6 Differences in community weighted mean (CWM) values (± standard error) of ecosystem preference functional trait categories in ciliate communities between control (C) and restored streams (N) during the two-months immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.
In the text

	Fig. 7 Differences in community weighted mean (CWM) values (± standard error) of habitat preference functional trait categories in ciliate communities between control (C) and restored streams (N) for the two-months immersion period. Functional trait categories showing significant differences between C and N sites are indicated by asterisk (*) in the legend.
In the text

Fig. 8 Redundancy analysis (RDA) showing the relationships between ciliate species (red arrows) and environmental parameters (blue arrows) among control (C) and restored (N) sites for the one-month immersion period. Abbrevations: N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; AC = Aspidisca cicada; CU = Chilodonella uncinata; CM = Cinetochilum margaritaceum; CYC = Cyclidium sp.; HO = Holophrya ovum; PV = Platyophrya vorax; PG = Pleurotricha grandis; PA = Pseudochilodonopsis algivora; TC = Trithigmostoma cucullulus; VAC = Vorticella aquadulcis-complex.

In the text

Fig. 9 Redundancy analysis (RDA) showing the relationships between CWM values (red arrows) and environmental parameters (blue arrows) among control (C) and restored (N) sites for the one-month immersion period. Abbrevations: N = nitrites; COD = chemical oxygen demand; OM = organic matter; IM = inorganic matter; CHL = chlorophyll a concentration; FSB = food source “bacteria”; FSPPSM = food source “phagotrophic protists, small metazoans”; EPLF = ecosystem preference “lentic (freshwater)”; EPLLF = ecosystem preference “lentic and lotic (freshwater)”; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”; HPP = habitat preference “planktonic”; MLFS = mode of locomotion “free-swimming”; MLC = mode of locomotion “crawling”.

In the text

Fig. 10 Redundancy analysis (RDA) showing the relationships between CWM values (red arrows) and environmental parameters (blue arrows) among control (C) and restored (N) sites for the two-months immersion period. Abbrevations: COD = chemical oxygen demand; OM = organic matter; CHL = chlorophyll a concentration; EPS = ecosystem preference “soil”; HPBP = habitat preference “benthos, periphyton”; HPB = habitat preference “benthos”.

In the text