Characteristics of the Stone Crayfish Population along a Disturbance Gradient—A Case Study of the Kustošak Stream, Croatia

Abstract

The stone crayfish, Austropotamobius torrentium, is a native European freshwater species sensitive to habitat alterations such as stream channelization and urban development, putting its populations at risk. This study aims to assess how habitat alteration and subsequent water quality changes affect the stone crayfish population in the Kustošak Stream (Croatia) through three selected sites under different levels of anthropogenic disturbance: (1) undisturbed; (2) recently modified, and (3) long-time modified sites. At each site, crayfish were captured, measured, and marked to estimate population size, structure, and crayfish condition. Additionally, we examined whether water quality (18 measured physicochemical parameters) affects relative crayfish abundance. We also used mitochondrial sequencing and microsatellite genotyping to assess species genetic diversity and population connectivity and to compare results among the sites. The results showed that habitat alteration caused an overall reduction in population abundance and changes in size structure; however, it had not yet resulted in detectable genetic differences. Partial least squares regression showed that crayfish abundance was affected by the physicochemical parameters of water, including, among others, oxygen, calcium ions, nitrates, pH, and water temperature. We discuss our findings in the context of the effects of anthropogenic disturbance on the viability and persistence of this EU priority species.

1. Introduction

Austropotamobius torrentium (von Paula Schrank, 1803), the stone crayfish (Figure 1), is a native European crayfish species [1]. It is the smallest among all native European crayfish, and it is considered an ecosystem engineer and a keystone species in freshwater ecosystems [2]. In recent decades, native European crayfish species, including the stone crayfish, have begun facing significant population declines [1,3] caused by water pollution [2], climate change [3], the presence of invasive alien crayfish species and their pathogens [4], and habitat alteration [2]. The stone crayfish is highly vulnerable to all of these pressures, as it is the species with the lowest reproductive output of all native species [5] and is determined to become especially endangered in the future due to climate change and invasive crayfish range expansion scenarios [6]. Its status in Europe is assessed as unfavourable-inadequate (U1) or unfavourable-bad (U2) in all biogeographical regions, according to reporting under Article 17 of the Habitats Directive [7]. Finally, the stone crayfish is listed in Appendix III of the Bern Convention [8], and in Annexes II and V of the EU Habitats Directive (92/43/EEC) [9]. Following the EU Habitats Directive Annex II, the stone crayfish is a priority species for which each EU member state is required to designate specific Natura 2000 sites (special areas of conservation, SACs), while other listings also aim to ensure its conservation and the appropriate management of its populations.

Stone crayfish populations in Croatia, similar to other populations across Europe, are endangered [3]: previous studies estimated that 28% of previously existing populations have disappeared [3,10]. In Croatia, due to sufficient and long-term monitoring data, the stone crayfish is classified as a vulnerable species according to the IUCN [11], and it is strictly protected on a national level by the Nature Protection Act (NN 80/13, NN 144/2013) which should ensure prohibition of its capture and prevent damage or destruction of its habitat.

Finally, the centre of genetic diversity of the stone crayfish is situated in the western part of the Balkan Peninsula, with 7 out of 9 described divergent mitochondrial lineages present in the karstic areas in Croatia [12,13,14]. This cold-adapted species usually inhabits areas with small pristine water bodies at higher altitudes [1]. One of those areas, the Medvednica Nature Park, is a Natura 2000 site situated above Croatia’s capital, the city of Zagreb. The park contains 45 streams which support a high number of the stone crayfish populations [15,16,17]. However, while upper courses of those streams are located within the park, large areas of the streams are located outside of the park’s borders, meaning that an unknown proportion of stone crayfish populations is present in urban areas of the city of Zagreb [18], outside the park and the Natura 2000 network. One of the aforementioned streams is the Kustošak Stream [15,16,17,18], for which 68.85% of the stream is located outside the park’s borders, where it flows through populated areas under elevated anthropogenic pressures. Additionally, the stream has been hydromorphologically changed through channelization, stream bed alterations, and building of dams as a part of municipal stream management and flood protection. Furthermore, part of the Kustošak Stream outside the park’s borders was modified recently during construction work, when the stone substrate and vegetation were removed from the stream bed, resulting in a reduction in the natural shade of the stream. The construction work was completed according to all regulations given by government authorities, but has resulted in potentially damaging habitat alterations for the stone crayfish. Here, we aimed to assess how these recent habitat modifications and subsequent water quality changes affected the stone crayfish population in the Kustošak Stream and whether these effects can be observed at the genetic level. Even though the conducted research encompasses a short period of time, it reflects the effect of anthropogenic pressure on natural stone crayfish populations. In other words, habitat alteration and loss can negatively impact the demographics of natural populations, reducing population size, gene flow, and genetic diversity, thereby decreasing the likelihood of population persistence, which will be more pronounced over a longer period of time. Such reductions lead to inbreeding and make species more vulnerable to extinction because lower diversity decreases adaptability and species long-term persistence [19,20].

2. Materials and Methods

2.1. Study Area and Data Collection

For this study, we selected three spatially close sites (less than 500 m) along the lower part of the Kustošak Stream, outside the Medvednica Nature Park borders, and under different levels of the anthropogenic disturbance: (1) undisturbed, on the border with the Medvednica Nature Park; (2) recently modified, with stone substrate and vegetation removed by construction work in the riverbed; and (3) long-time modified (channelled), with concrete banks (Figure 1). Sites 2 and 3 were situated in a populated area, indicating potentially higher levels of anthropogenic disturbance (Figure 2).

The field work was conducted during a week-long period in September 2021, during a period of increased crayfish activity [2]. The crayfish were trapped using 10 baited handmade traps per site. The traps were made of two plastic bottles (1.5 L) coupled together, with funnel entrances (made from plastic net, mesh size = 2 mm) at each end. The traps were baited with sausages and left in the water overnight for a total of 7 days, regularly checked every day. Upon capture, the sex of each individual was determined, and the samples were measured (weight and total length). Weight (w [g]) was measured using a dynamometer (Pesola 100 g [accuracy ± 0.3 g], Switzerland). The crayfish total body length (TL [mm], from the tip of the rostrum to the end of the telson), was recorded per each crayfish using a digital caliper (Alpha Tools [accuracy ± 0.02 mm], Franklin, NJ, USA). From each individual, one pereopod was sampled and stored in 96% ethanol at 4 °C until further DNA isolation. The survival of crayfish was not compromised because the sampled pereopods will regenerate upon the next moulting. Sampling was conducted in accordance with ethical standards, and the required permission was obtained from the Ministry of Economy and Sustainable Development of the Republic of Croatia (UP/I-612-07/21-48/129). After sampling, the crayfish were marked as described in Section 2.2. and released back into the stream. Finally, at each site, 18 physicochemical water parameters were measured (

Table 1. Physicochemical parameters (18) of water measured in the study.

	Parameter	Abbreviation	Measuring Unit
1	Dissolved oxygen	O2	mg/L
2	Water temperature	w. temp.	°C
3	Water pH	pH	/
4	Oxidation reduction potential	ORP	mV
5	Water turbidity	turbidity	NTU
6	Chlorides	Cl−	mg/L
7	Nitrates	NO3−	mg/L
8	Sulphates	SO42−	mg/L
9	Sodium ions	Na+	mg/L
10	Potassium ions	K+	mg/L
11	Magnesium ions	Mg2+	mg/L
12	Calcium ions	Ca2+	mg/L
13	Ammonium	NH4+	mg/L
14	Chemical oxygen demand	COD	mg/L
15	Nitrate nitrogen	NO3−N	mg/L
16	Nitrites	NO2−	mg/L
17	Total nitrogen	TN	mg/L
18	Total phosphorus	TP	mg/L

). Parameters 1–5 were measured directly in the field [O₂ and water temperature were measured using an oximeter (OXPB-1, Lutron YK-2005WA, Taiwan); pH was measured using a pH meter (PE-03, Lutron YK-2005WA, Taiwan); ORP was measured using a pH and redox meter (SenTix 60, WTW pH 3310, Weilheim, Germany); water turbidity was measured using a turbidimeter (Hach 2100Q, Loveland, CO, USA)]. For analyses of parameters 6–18, water samples were collected in bottles, washed three times before filling them to the top, transported on ice, and stored at −20 °C until the analyses were performed. Concentrations of Cl⁻, NO₃⁻, SO₄²⁻, Na⁺, K⁺, Mg²⁺, and Ca²⁺ were determined in the lab using a two-channel analytical ion chromatograph ICS6000 (Thermo Scientific, Waltham, MA, USA), with suppressed conductometric detection and electrolytic preparation of eluents on both channels. For anion separation, an AS11-HC-4 µm (2 × 250 mm, Thermo Scientific, Waltham, MA, USA) column was used, with a suitable precolumn AG11-HC (Thermo Scientific, Waltham, MA, USA) and KOH as the eluent, while for cation separation, a CS16–4 µm column (2 × 250 mm, Thermo Scientific, Waltham, MA, USA) was used, with a suitable precolumn CG16 (Thermo Scientific, Waltham, MA, USA) and MSA as an eluent; conventional analytical determinations for concentrations of NH₄⁺, COD, NO₃⁻N, NO₂⁻, TN, TP, and COD were conducted with Hach Lange Cuvette Tests (LCK 304, 1414, 339, 341, 338, 350) using a DR3900 spectrophotometer (Hach, Loveland, CO, USA).

2.2. Data Analyses

2.2.1. Population Size

Upon capture, each individual was marked with a number on the carapace using an oil-based marker pen. This approach enabled us to estimate population size per site by two mark-recapture methods: the Schnabel method [21],

$$\mathrm{N}=\frac{\sum \left({\mathrm{M}}_{\mathrm{t}}{\mathrm{C}}_{\mathrm{t}}\right)}{\left(\sum {\mathrm{R}}_{\mathrm{t}}\right)+1},$$

and the Schumacher–Eschmeyer method [22],

$$\mathrm{N}=\frac{\sum \left({\mathrm{C}}_{\mathrm{t}}{\mathrm{M}}_{\mathrm{t}}^2\right)}{\sum \left({\mathrm{M}}_{\mathrm{t}}{\mathrm{R}}_{\mathrm{t}}\right)},$$

where: N = estimation of the number of individuals in a population, M_t = number of individuals marked in sample t, C_t = total number of individuals caught in sample t, R_t = number of individuals already marked when caught in sample t.

For the Schnabel method, the upper and lower confidence limits of the estimated population size were determined,

$$\frac{1}{\mathrm{N}}\pm \mathrm{t}\times \sqrt{\mathrm{var}\left(\frac{1}{\mathrm{N}}\right)},$$

where t = 1.96 (critical value of Student’s t-test, with 95% confidence).

The relative abundance of crayfish was also estimated as the catch per unit effort (CPUE: number of crayfish caught per trap per night).

2.2.2. Size Class Frequency

Collected data regarding total length were used for the analysis of the size structure per site. Based on total length, the crayfish were divided into 5 size classes (0–2.99 cm, 3–4.99 cm, 5–6.99 cm, 7–8.99 cm, and >9 cm). The size classes were arbitrarily set, similar to the method of Maguire and Klobučar [23].

2.2.3. Fulton’s Condition Factor

Data on weight and total length were used for the calculation of the condition of the crayfish. Fulton’s condition factor (FCF [g/mm³]) was calculated as an indicator of crayfish fitness [24],

$$\mathrm{FCF}=\frac{\mathrm{w}}{{\mathrm{TL}}^3}\times 100,$$

where: w = weight of the individual, and TL = total length of the individual.

2.2.4. Analysis of Physicochemical Parameters of the Water

Collected data regarding the physicochemical parameters of the water were used to examine whether a correlation exists between these and the CPUE.

2.3. Statistical Analyses

Crayfish characteristics and population parameters (CPUE, TL, weight, and FCF) were first tested for normality using the Shapiro–Wilk test. Since the data did not meet the assumptions necessary for the use of parametric tests, even after transformations, their non-parametric analogues were used instead. The Mann–Whitney U test was used to determine whether differences in TL, weight, and FCF exist between males and females at each site. If there were no statistically significant differences between the sexes, males and females were pooled in the subsequent analyses. However, only in the case of FCF, analyses were performed separately for males and females due to observed significant differences between sexes (resulting from sexual dimorphism). Differences in all measured crayfish characteristics and population parameters between the sites were analysed using Kruskal–Wallis ANOVA with Dunn’s post hoc test. Analyses were conducted in TIBCO Statistica version 14.0.1.25 (Palo Alto, CA, USA).

The partial least squares regression (PLS-R) method was used to examine the correlation between the physicochemical parameters of the water (explanatory variables, X), and the CPUE at each site (response variable, Y). Initially, the model quality indices were calculated for two components to validate the performance of the model. Then, a correlation radar was generated in order to visually represent the relationship between explanatory and response variables. Lastly, to indicate the importance of the explanatory variables for the response variables, variable importance in the projections (VIPs) was calculated. The analysis was performed using XLSTAT version 2021.3.1.1189 software for data analysis and visualisation (Addinsoft, Paris, France, Microsoft Excel). In all performed analyses, the significance level was set at p < 0.05.

2.4. Genetic Diversity and Structure

The total genomic DNA was extracted from the pereopod muscle tissue with a GenElute Mammalian Genomic DNA Miniprep kit (Sigma-Aldrich, St. Louis, MO, USA), following the manufacturer’s protocol.

The phylogenetic position and haplotype diversity of the studied crayfish were assessed by mitochondrial DNA analysis (cytochrome oxidase subunit I gene), while the fine genetic diversity and population structure were studied by microsatellite DNA analysis.

2.4.1. Mitochondrial DNA Analyses

The mitochondrial COI gene was amplified and sequenced with universal barcode primers LCO-1490 and HCO-2198 [25], allowing for comparison with previously published stone crayfish sequences. PCR, purification, and sequencing were performed according to the methods of Lovrenčić et al. [14]. The sequences were edited and aligned in Bioedit v. 7.2.5 [26]. The final COI alignment did not contain any length variants or ambiguous sites and included all previously published sequences, with the final length of 582 bp. The sequences were subsequently collapsed to unique haplotypes using FaBox [27] in order to associate haplotypes obtained in the present study with the haplotypes obtained in previous research. Newly obtained haplotypes were submitted to GenBank (accession numbers OQ048681 to OQ048683). Additionally, Austropotamobius pallipes (Lereboullet, 1858) was chosen as an outgroup (GenBank Accession numbers: KX369673, KX369674).

The phylogenetic assignment of the studied individuals was estimated using Bayesian analysis (BA) in MrBayes ver.3.2.6 [28]. The optimal model of nucleotide evolution for the COI dataset was HKY + I + G, selected under the Bayesian information criterion (BIC) using the jModelTest 2.1.10 [29]. The BA was carried out with priors set according to the suggested model; two separate runs with four Metropolis-coupled Monte Carlo Markov chains (MMCM) were performed for 10,000,000 generations, and trees were sampled every 1000 generations. After checking the congruence (ESS values > 200) with Tracer v1.7.1 [30], the first 25% of the sampled trees were eliminated as burn-in, and a consensus tree was constructed, with nodal values representing the posterior probabilities (values ≥ 0.95 were considered supported).

2.4.2. Microsatellite DNA Analyses

The population genetic diversity and structure were estimated using seven polymorphic microsatellite loci (ATM78, ATM79, ATM64, AT1, AT37, ATOR37, and Aas3040) and following the PCR protocols according to Lovrenčić et al. [31]. Genotyping was performed by capillary electrophoresis using internal GeneScan 600 LIZ Size Standard in Macrogen, Inc. (Seoul, South Korea), while microsatellite genotypes were scored and double-checked manually using Peak Scanner software v.1.0 (Applied Biosystems, Waltham, MA, USA). Microsatellite loci were assessed for potential presence of genotyping errors due to scoring of null alleles, stuttering, and large allele dropout using MICRO-CHECKER v.2.2 [32].

The within-population genetic diversity was assessed by calculating the proportion of polymorphic loci (P), mean number of alleles (N_A), number of private alleles (A_PR), observed heterozygosity (H_O), unbiased expected heterozygosity (H_E), and deviation from the Hardy–Weinberg equilibrium (HWE) for each site across all loci and using GenAlEx v.6.51 [33]. Allelic richness (A_R) was calculated and corrected for sample size by the rarefaction method using FSTAT v.2.9 [34]. The same software was used for the calculation of the inbreeding coefficient (F_IS) and the levels of genetic connectivity of each site by pairwise comparisons of F_ST values and evaluated using 1000 permutations. Potential signatures of recent bottlenecks were tested under three different mutational models: the infinite allele model, the stepwise mutation model, and the two-phase model in BOTTLENECK v.1.2 [35]. Significant deviations from mutational-drift equilibrium were tested using the Wilcoxon sign rank test with 10,000 simulations.

The genetic structure among the studied sites and the assembling of individuals into genetic clusters was assessed using the Bayesian model-based clustering approach, as implemented in STRUCTURE v.2.3 [36]. The conditions performed were 10 runs for each potential genetic cluster (K) between 1 and 3, with a 100,000 burn-in period followed by 100,000 Markov Chain Monte Carlo iterations, using correlated allele frequencies under a straight admixture model. The best number of clusters within the dataset (K) was determined using the Evanno method (ΔK method) [37] in STRUCTURE HARVESTER [38]. STRUCTURE graphical results were plotted with CLUMPAK [39].

3. Results

3.1. Descriptive Statistics

A total of 204 individuals were caught, 107 at site 1, 56 at site 2, and 41 at site 3. Descriptive statistics of measured crayfish parameters, as well as Fulton’s condition factor and CPUE per site and sex, are presented in

Table 2. Results of descriptive statistical analysis (mean ± standard deviation) of parameters: total length (TL), weight (w), Fulton’s condition factor (FCF), and catch per unit effort (CPUE) of male and female individuals at different sites (site 1 = undisturbed site, site 2 = recently modified site, and site 3 = long-time modified (channelled) site).

Site	Sex	TL [mm]	w [g]	FCF [g/mm3]	CPUE
1	male	49.97585 ± 13.76868	6.57317 ± 6.50822	0.00435 ± 0.00086	2.19756 ± 0.96475
1	female	54.22271 ± 12.92085	6.87879 ± 4.47218	0.00390 ± 0.00074	2.51212 ± 0.99591
2	male	62.50467 ± 9.82911	10.5000 ± 5.34499	0.00402 ± 0.00052	1.13000 ± 0.30529
2	female	65.90153 ± 12.56042	11.30769 ± 5.50468	0.00417 ± 0.00318	1.10370 ± 0.31191
3	male	59.93421 ± 11.54536	10.23684 ± 5.46038	0.00556 ± 0.00659	0.94210 ± 0.24339
3	female	65.24272 ± 8.99621	10.54545 ± 4.51692	0.00358 ± 0.00046	0.95000 ± 0.25023

, while descriptive statistics for the measured physicochemical parameters is shown in Supplementary Tables S1 and S2.

3.2. Estimation of Total and Relative Population Size

The total estimated number of crayfish using both mark-recapture methods (the Schnabel method and the Schumacher–Echmeyer method) followed the recorded disturbance gradient, and the highest number of crayfish was estimated at the undisturbed site (site 1), followed by the recently modified site (site 2), and the lowest number was estimated at the long-time modified site (site 3) (

Table 3. Estimation of total population size per site with two mark-recapture models, the Schumacher–Eschmeyer method and the Schnabel method. For the Schnabel method, the lower and upper confidence limits are shown.

	Schumacher–Eschmeyer Method	Schnabel Method
site	estimation	estimation	lower confidence limit	upper confidence limit
1	328	346	251	555
2	161	98	77	139
3	87	60	46	86

).

These results also reflected differences in recorded CPUE, which exhibited significant differences between sites (Kruskal–Wallis ANOVA: H (2, N = 204) = 121.7805, p < 0.001). The CPUE at the undisturbed site 1 was significantly higher than the CPUE at the two modified sites (Figure 3).

3.3. Size Structure and Crayfish Condition

Crayfish size also differed significantly between the examined sites: total length (Kruskal–Wallis ANOVA: H (2, N = 204) = 35.77619, p < 0.001) and weight (Kruskal–Wallis ANOVA: H (2, N = 204) = 39.11740, p < 0.001), with undisturbed site 1 again showing significant differences in size of caught individuals compared to the other two modified sites (Figure 4).

The size class frequency of crayfish per site is shown in Figure 5. All size classes were present only at the undisturbed site, while the two modified sites lacked smaller size classes.

Crayfish condition, measured as Fulton’s condition factor, exhibited no significant differences between the sites, for either males or females (p > 0.05).

3.4. Effect of Physicochemical Parameters of Water Quality on CPUE

Using the PLS-R method, the relationships between crayfish abundance expressed as CPUE (response variable, Y) and various physicochemical parameters of the water (explanatory variables, X) were analysed. As for the model quality indices R²X and R²Y in component 1, 70.1% of the variance in the set of explanatory variables (X) was used to explain 42.6% of the variance for the response variable (Y). In component 2, 100% in the set of explanatory variables (X) was used to explain 43.6% of the response variable (Y). Q², as a measure of good prediction, showed that component 1 contributed 31.2% and component 2 contributed 30.4% to the model quality. The relationship between the explanatory and response variables is visually represented in the form of a correlation radar (Figure 6), with positively correlated variables close together and negatively correlated examples far apart. PLS-R analysis showed that K⁺, NO₃⁻, O₂, and Ca²⁺ were the most important parameters that positively influenced CPUE, followed by water turbidity, COD, ORP, Cl⁻, TN, and NO₃⁻, N, while the most pronounced negative correlation was found for Mg²⁺, Na⁺, pH, and water temperature, followed by NO₂⁻, NH₄⁺, TP, and SO₄²⁻ (Figure 6, Supplementary Table S3). The majority of the parameters (11 in total) had VIP values > 1, meaning that they are considered highly relevant in explaining the CPUE and contribute significantly to the model (Supplementary Figure S1).

3.5. Genetic Diversity and Structure

3.5.1. Mitochondrial DNA Analyses

We obtained 33 new sequences from the studied stream, belonging to six unique COI haplotypes, three of them (AT154, AT156, AT157) recorded for the first time (Supplementary Table S4, Figure 7). Haplotypes from this study were recovered within the well-supported monophyletic phylogroup Central and South Eastern Europe (CSE) sensu Klobučar et al. [12] (Figure 7). Regarding haplotype diversity at each of the studied sites, we recorded six haplotypes at site 1 (AT9, AT119, AT143, AT154, AT156, AT157), five at site 2 (AT9, AT119, AT143, AT154, AT156), and four at site 3 (AT9, AT143, AT156, AT157) (Supplementary Table S4).

3.5.2. Microsatellite DNA Analyses

A total of 36 alleles were observed across the seven microsatellite loci, ranging from three (loci AT1 and Ator37) to nine (loci ATM64 and AT37). No genotyping error due to stuttering or large allele dropout was apparent across the loci. Null alleles were detected only within locus AT37 at sites 2 and 3. However, since no locus showed null alleles across all populations, a bias in our analyses due to null alleles was unlikely, and all subsequent analyses were conducted using the original dataset.

The summary statistics of the genetic diversity indices for the studied sites are shown in

Table 4. Summary results across seven microsatellite loci of population genetic diversity of the studied stone crayfish sites. P—proportion of polymorphic loci, N—sample size, N_A—average number of alleles/locus, A_R—allelic richness, A_PR—number of private alleles, uH_E—unbiased expected heterozygosity, H_O—observed heterozygosity, F_IS—inbreeding coefficient, and P_HWE—probability of deviation from Hardy–Weinberg equilibrium after Bonferroni adjustments (not significant—ns).

Site	N	P	NA	AR	APR	uHE	HO	FIS	HW
1	15	1	4.00	3.73	5	0.597	0.524	0.126	ns
2	13	1	3.71	3.56	3	0.496	0.465	0.065	ns
3	12	1	3.86	3.73	5	0.569	0.502	0.122	ns

. Overall, the microsatellite markers showed a high level of polymorphism and an acceptable level of genetic diversity (Table 4). Genetic diversity, expressed as P, N_A, A_PR, A_R, and H_O, was the highest at site 1, followed by site 3, while the lowest genetic diversity was indicated for site 2 (Table 4). A total of 13 private alleles were found, ranging from 3 (site 2) to 5 (sites 1 and 3). The inbreeding coefficient (F_IS) per population was low, ranging between 0.07 at site 2 to 0.12 at sites 1 and 3, indicating a slight homozygote excess/heterozygote deficit. Moreover, all sites were in Hardy–Weinberg equilibrium (Table 4). Bottleneck analysis did not reveal consistent signs of contraction of population size and recent bottleneck, according to the three mutational models tested (

Table 5. Probability (bold indicates significant p-values; p < 0.05) of bottleneck for stone crayfish populations using Wilcoxon sign rank test under three different mutational models: infinite allele model (IAM), stepwise mutation model (SMM), and two-phase model (TPM). Significant differences are marked in bold.

Site	Wilcoxon Sign-Rank Test 1-Tail			Wilcoxon Sign-Rank Test 2t
	IAM	TPM	SMM	IAM	TPM	SMM
1	0.03	0.15	0.41	0.05	0.30	0.81
2	0.34	0.71	0.95	0.69	0.69	0.30
3	0.04	0.34	0.66	0.08	0.69	0.81

).

Genetic differentiation among three sites showed no genetic differentiation, with pairwise F_ST = 0 for all site pairs. Regarding population genetic structure, the Evanno method revealed that the optimal number of clusters was two (ΔK = 2). The Bayesian clustering analysis in STRUCTURE showed a high level of admixture, with all individuals exhibiting the same genetic architecture and ancestry from both genetic clusters (Figure 8).

4. Discussion

The results of this research provide another example of how anthropogenic pressure negatively affects populations of endangered freshwater species. Stone crayfish populations are facing population declines driven by the introduction of invasive alien species, climate change, and anthropogenic disturbance of their habitats [1,3,40]. Recent studies have revealed that anthropogenic impacts on freshwater habitats pose the highest threat for the stone crayfish populations in Croatia, while in the future, currently suitable habitats are predicted to be lost under different climate change scenarios [31]. Since natural habitats of the stone crayfish are usually isolated and often fragmented by habitat destruction [41], and since the stone crayfish has a low dispersal capacity, there is a weak to no possibility for natural repopulation, if a certain population disappears [42]. Although the study was conducted in a short period of time, and we are aware that long-time monitoring would provide a clearer representation of the state of the stone crayfish population in the Kustošak Stream, we nonetheless gained initial insights into the effects of anthropogenic pressures on a stone crayfish population located outside of the protected areas. Habitat alterations, including stream channelization and urban development, have a drastic negative impact on the freshwater ecosystems [43]. There are numerous documented records of population declines due to habitat alteration, including freshwater taxa such as fish and mussels [44,45]. Furthermore, human activities have also affected the quality of the physical habitat for crayfish, subsequently resulting in reduced crayfish abundance [46,47]. For example, from 67 threatened North American crayfish species, 52% were reported to be imperiled by habitat alteration [48]. In addition, the abundance of crayfish has often been positively correlated to the abundance of shelters [49], which are usually provided by heterogeneous (natural) habitats, that include cobbles, macrophytes, and bank vegetation. Previous studies on the white-clawed crayfish, A. pallipes, showed that the presence of bank vegetation provides valuable protection for smaller size classes, resulting in positive effects on the abundance of this species [50]. The results of our study corroborate that habitat alteration (both recent and historical) potentially caused reduction in the stone crayfish population abundance and changes in size structure. Both population size estimation methods (Schnabel–Schumacher and Eschmeyer), as well as relative crayfish abundance measured using CPUE, show the decline in abundance in population at the two anthropogenically impacted sites, the recently modified and long-time modified sites in the Kustošak Stream. These significant differences were observed in a short stretch (ca. 1 km long), suggesting that they could be a result of the habitat alterations and likely subsequent water quality changes, reduced natural shading, and availability of shelters. The observed changes in abundance were accompanied by changes in the population size structure and measured morphological characteristics. According to the size class frequency histograms, all size classes were observed only at the undisturbed site, which could indicate a healthy population [51]. There, crayfish were smaller compared to individuals in the recently modified and long-time modified site, where mature crayfish prevail. Considering that the stone crayfish reach sexual maturity at 50–60 mm [52], we can conclude that sexually mature individuals are present at all three sites, but anthropogenically impacted sites lack smaller individuals. The research regarding the movement rates of stone crayfish has shown that it is a more sedentary species compared to other crayfish [53,54]. Stone crayfish were recorded to have a maximum movement of 133 m in 55 days [54]. Since the three sites in our study are more distant (499.67 m between sites 1 and 2 and 236.73 m between sites 2 and 3) and capturing was done daily, we presume that the changes in the size class frequency were not a result of migration of bigger individuals between sites. Rather, the lack of smaller individuals at anthropogenically impacted sites could be due to stone substrate and bank vegetation removal, as well as streambed alteration, which changed the optimal living conditions for the stone crayfish population present there, making those sites insufficient for supporting all size classes. The results of our study showed no significant impact of habitat alteration on crayfish conditions, since no differences were observed between sites. This could be related to the timing of the research, considering that during the mating season, adult, reproductively active crayfish are at peak condition [52,55]. Previous studies defined the optimal physicochemical parameters of water supporting stone crayfish populations as pH 5.0–8.6, water temperature 11–26 °C, O₂ over 4.0 mg/L, Ca²⁺ 7.0–70.0 mg/L, Mg²⁺ 2.6–21.0 mg/L, Cl⁻ up to 16.7 mg/L, NO₃⁻ up to 44 mg/L, and NO₂⁻ up to 1.7 mg/L [56]. Most of the physicochemical parameters measured in our study met the abovementioned ranges. Furthermore, PLS-R analyses showed the strongest positive correlation of relative crayfish abundance (CPUE) with K⁺, NO₃⁻, O₂, and Ca²⁺ and the strongest negative correlation with Mg²⁺, Na⁺, pH, and water temperature. The oxygen demands of the stone crayfish have been reported to be higher in comparison to other native European crayfish species due to their occurrence in the upper parts of streams distinguished by water with lower temperature [57,58], explaining the positive correlation of CPUE with O₂ and the negative correlation with water temperature. Similarly, the negative correlation of CPUE with pH was expected, as lower pH levels may be toxic in the long term [58]. The positive correlation with calcium ions was an expected outcome, which is typical for crustaceans in general due to their calcified exoskeleton [59].

Contrary to population size and structure results, genetic data revealed no significant effects of habitat alteration on the genetic diversity within the studied area, as genetic variation is still evident based on both mitochondrial haplotypes and nuclear alleles. While levels of genetic diversity were similar according to the microsatellite analysis, the number of recorded COI haplotypes differed among sites, with a higher number in the undisturbed site (1) compared to anthropogenically modified downstream sites (2 and 3).

Phylogenetic reconstruction indicated that three novel haplotypes nested within formerly recognised mitochondrial phylogroup CSE sensu Klobučar et al. [12]. This was expected, since the geographical position of the studied population falls into the region where the CSE mitochondrial phylogroup is distributed, and previous studies have established a strong phylogeographical structure of the stone crayfish, with numerous haplotypes described within each phylogroup [12,14]. This high level of haplotype diversity emerged from the past geo-climatic events, which altered river drainage patterns in the region [12], and low contemporary levels of gene flow among the populations [31]. However, although overall haplotype diversity across studied sites was high and comparable to other CSE populations in the region [14], we observed a decreasing number of haplotypes as habitat alteration increased.

Microsatellites analyses revealed moderate within-population genetic diversity and low differentiation among upstream and downstream sites, suggesting a high level of genetic connectivity, which indicates that the stone crayfish has maintained gene flow despite documented habitat disturbances. Similar results were found for various species, which did not show different levels of genetic diversity in altered and natural habitats [60,61,62]. Despite possible population fragmentation and human disturbance, all sites exhibited similar levels of genetic diversity when compared to other Croatian populations for which microsatellite data have been generated [31]. Generally, lower levels of genetic diversity occurring within the stone crayfish populations suggests that this species may have experienced severe habitat fragmentation or population size reduction, both in the past and recently [31]. The low genetic diversity, along with environmental instability, caused primarily by anthropogenic activities, can decrease species fitness and affect long-term survival [19,63]. Although our study showed that heterozygosity and inbreeding were not strongly affected by habitat alteration, all three sites showed lower observed heterozygosity than expected, suggesting that inbreeding or genetic drift might have an impact in the future. Inbreeding in a small, isolated population of the stone crayfish elevates its extinction risk, giving it substantial conservation significance [19,64]. The genetic clustering approach used to assess population structure showed high admixture and individuals exhibiting assignment to two ancestral genetic groups. Moreover, the geographical proximity of sites probably increased the gene flow and homogenized the present genetic variation, which resulted in all individuals exhibiting a similar genetic background.

Although habitat alteration did not yet result in detectable genetic consequences, we strongly recommend regular genetic monitoring that will establish genetic isolation trends of the studied population and identify early warning signs of further population decreases. Numerous studies have shown that a vast number of species experience genetic erosion due to habitat loss [65,66,67,68]. The negative effects of habitat destruction include increased the isolation of populations or species and the creation of small, isolated subpopulations. In the case of the stone crayfish, additional longer periods of isolation can lead to reduced population size and genetic diversity, consequently causing cascading effects throughout ecosystems due to its keystone role in freshwater habitats [19]. Long-term studies are thus required to fully understand the impacts of the anthropogenic disturbance on the population of the stone crayfish and its genetic diversity in this stream; thus, we advocate that such monitoring should be performed and maintained regularly, since our study indicates that such effects may exist and are visible with short-term monitoring. In conclusion, conservation of this species in general, and specifically in the Kustošak Stream, should focus on maintaining favourable habitat conditions and genetic connectivity. Additionally, we suggest that no further habitat alterations should be planned or performed in this stream, as they will further imperil population connectivity and potentially result in the disappearance of the stone crayfish from these downstream (urbanized) areas outside the park’s borders. Moreover, we strongly suggest that any further stream bed reconstruction activities planned in any of the park’s 45 streams inhabited by the stone crayfish should be avoided or planned only after long-time careful analysis of the current population status.