Life History Traits of the Stygophilous Amphipod Synurella ambulans in the Hyporheic Zone of the Lower Reaches of the Upper Sava River (Croatia)

Abstract

Estimating the growth and population dynamics of crustaceans is important for understanding the role of species in ecosystems and for conservation and monitoring purposes. This study investigated temporal and spatial variation in population dynamics (growth rate, instantaneous mortality rates, and longevity) of the stygophilous freshwater amphipod Synurella ambulans (F. Müller, 1846), and influencing environmental factors. Seasonal sampling was conducted from December 2018 to October 2019 at two sampling sites (Jarun and Medsave) along the Sava River in northwestern Croatia. A Bhattacharya cohort analysis was applied to the length-frequency data using the FISAT software package. At least four cohorts were distinguished during the year. The possible influence of environmental factors on growth parameters was investigated using Canonical Correspondence Analysis (CCA), which revealed that the important environmental factors influencing S. ambulans abundance were water temperature and dissolved oxygen concentration. The average total body length (TBL) of females was 4.25 mm, while the TBL of males was 3.34 mm. The asymptotic length (L_∞) and growth coefficient (K) were estimated to be 6.30 mm and 0.92 year⁻¹ for females and 5.40 mm and 0.59 year⁻¹ for males, respectively, at the Jarun site. L_∞ and K were estimated to be 7.20 mm and 0.34 year⁻¹, respectively, for females and 4.00 mm and 1.20 year⁻¹, respectively, for males at the Medsave site. Estimated total mortality (Z) for each sex had higher values for males than females at both sampling sites. The lack of ecological knowledge on S. ambulans populations and relationships with invertebrate communities in the hyporheic zone of the large river ecosystems emphasises the importance of a detailed study for protecting this species and its vulnerable ecotonal groundwater-connected environment.

1. Introduction

Amphipods (class: Malacostraca; order: Amphipoda) belong to the most abundant, diverse, and taxonomically challenging crustacean order, including over 10,500 species occurring in various aquatic and even terrestrial ecosystems, but most live in marine environments [1,2,3]. However, amphipods are a significant component of various freshwater ecosystems, playing a crucial ecological role in nutrient cycling, bioindication, ecotoxicity testing, hosting parasites, and supporting numerous ecosystem services [4,5]. Due to their critical functional role as prey for many macroinvertebrates, fish, birds and even mammals, they are essential for transporting energy and nutrients to higher trophic levels [6], highlighting their ecological importance in aquatic food webs and connecting aquatic and terrestrial environments.

1.1. Distribution Patterns and Ecology of Synurella ambulans

Synurella ambulans (F. Müller, 1846) belongs to the family Crangonyctidae, a diverse, almost exclusively freshwater Holarctic family comprising 249 species and ten subspecies inhabiting mainly subterranean and surface biotopes with groundwater connections [3,7,8]. It may represent a species complex [9], based on its morphological variability [10,11,12]. Even though the species belongs to the cold-stenothermal, photonegative and thigmotactic organisms [13], it is a euryoecious amphipod species widely distributed throughout central, eastern and south-eastern Europe that can live in a wide range of surface ecosystems (lakes, springs, rivers, swamps, ponds) and subterranean realms (hyporheic zone (HZ), hypotelminorheic zone, phreatic interstitial zone, and caves), in both lowland and mountainous regions [12,14,15,16,17] up to 1600 m a.s.l. [18,19]. As a highly adapted subterranean water dweller, the species belongs to eustygophilous inhabitants [20], although the eyeless stygobiotic populations were also documented in European countries [8].

S. ambulans has been recorded in Europe in many countries: Albania [15], Austria [21,22], Belarus [23], Belgium [24], Bulgaria [11], Croatia [15,25], Czech Republic [26,27], Estonia [28], Germany [16], Hungary [15,29,30], Italy [31,32], Latvia [33,34], Lithuania [17], Montenegro [15,35], Poland [11,36], Romania [37,38], Russia [18], Serbia [15], Slovakia [26,39], Slovenia [20,40,41], Switzerland [4,42], Turkey [43], and Ukraine [44]. In Croatia, it has been found in the Dinaric karst region (intermittent Krčić spring in central Dalmatia) [45] and in the Pannonian region [25]. In addition, this amphipod inhabits the hypotelminorheic zone with underlying flysch at higher elevations in Croatia (Medvednica Mountain, 927 m a.s.l.) [19]. Synurella ambulans rarely inhabits deep lakes, but it was found in the karstic oligotrophic lake Vrana on the island of Cres at depths of 50 to 70 m [46]. It has also been recorded in the HZ of the Sava River in eastern Slovenia [40]. Recently published data confirm that S. ambulans is a suitable bioindicator of chronic metal contamination in the HZ of the Sava River, as it accumulates metals from the interstitial water from the HZ [47]. As a small burrower-like inhabitant of deeper sediment, the species is often overlooked during classical macrozoobenthos field surveys [48], and therefore directed and targeted research is crucial for understanding the ecology and life history traits of this species.

1.2. Amphipod Life History Strategy in the Subterranean Environment

The life history traits of subterranean species have been poorly studied due to many obstacles and challenges to accessing their habitats, difficulties in animal collecting, unsuccessful breeding in laboratory conditions, and the rarity of most subterranean species [49,50,51]. There are only a few articles addressing some aspects of the life history traits of S. ambulans [7,11,12], their mouthparts morphology [52], distribution patterns, e.g., [9,17,24,43], and taxonomy [12]. Despite their importance in hyporheic communities [40], there are still significant gaps in our knowledge of groundwater-connected crustaceans’ biology and life history strategies. Moreover, in the Slovenian and Croatian parts of the Dinaric karst, stygobiotic populations (races) of S. ambulans were recorded in assemblages with stygobiotic Niphargus spp. and the cave salamander Proteus anguinus [20,53]. Therefore, detailed data on the biological characteristics of the species are crucial to understanding its position in the functional structure of the groundwater-connected ecosystems (GWCE).

Knowledge of the population structure and dynamics of crustaceans is fundamental to understanding population ecology, the ecological stability of a species, and the protection of aquatic ecosystems [54]. In addition, the effects on growth and reproduction in amphipods can be measured as one of the endpoints in assessing the effects of contaminants [55]. Life history traits such as the number, age and size at maturity, growth pattern, longevity, timing of reproduction, size and sex ratio of offspring, and fecundity of species have been shown to be related and interconnected to many factors and subjects, from habitat geography and behavioural ecology to natural selection [51,56,57,58]. For example, species from cold climates have greater longevity and greater body size, while species from warm climates have shorter longevity and smaller body size [56]. Life histories are also related to lifestyle and adaptive abilities of crustaceans to specific environmental realms [50,51,56,58].

1.3. Methodology for Life History Analyses

Crustaceans lack typical age determination structures such as scales, otoliths, spines, or shells due to the complex process of ecdysis and restricted growth associated with moulting [59,60]. Because their hard exoskeleton does not grow, an estimation of growth and age structure is severely limited. Therefore, indirect approaches have been developed to calculate the age of crustaceans. In recent decades, several new methods have been used to analyse ageing in crustaceans, such as the widely used content of the brain pigment lipofuscin as an accurate indicator of age [61,62], the counting of endocuticle growth bands in gastric mills and eyestalks [63], and the reflex action mortality predictor (RAMP) method used to predict mortality [64]. The amount of lipofuscin pigment is related to physiological age rather than chronological age [65]. However, analysis of this fluorescent age pigment in amphipods (in slow-growing amphipods) has been shown to be a more accurate method for chronological age determination than body size analysis [62]. In some decapods and slow-growing amphipods, lipofuscin has been shown to be a more adequate age predictor than size [66,67,68,69].

One of the most commonly used tools for assessing population parameters of shrimps is FISAT (FAO-ICLARM Stock Assessment Tools) [70,71,72]. FISAT software (version 1.2.2, FAO, Rome, Italy) was initially developed for fish and is still widely used to assess fish stock exploitation, e.g., [73,74,75], or to study changes in the growth rate of fish [76]. The Von Bertalanffy growth function (VBGF) [77] is an appropriate mathematical model used to describe the growth patterns of crustaceans during their life span. Several studies have used this model to measure the growth of crustaceans, such as economically important shrimp [78,79], amphipods [80,81,82,83], mysids [84], isopods [85,86], and decapods [87,88,89,90,91].

Knowledge of the life history of European gammarids in freshwater environments is sparse and fragmentary, with only a few studies in recent years [92]. Nowadays, analyses of population dynamics in crustaceans include not only life history traits (e.g., body size, population density, breeding period, and sex ratio), but a much broader dataset that includes analyses of longevity, mortality, growth rate, age information, and number of generations (cohorts).

The temporal and spatial changes in environmental factors such as temperature and salinity can strongly influence crustaceans’ growth rates and ageing analysis, as specimens of the same age may belong to different length classes [62,93]. Unfortunately, age analysis based on size-frequency histograms in crustaceans can be a density-dependent method in addition to age [94]. Nevertheless, length-frequency analysis is fast, inexpensive, and efficient, especially for small, non-commercial crustaceans. Moreover, neurolipofuscin as an age marker, mostly used in commercial decapods or larger amphipods, requires the dissection of the brain [66,94], which is very difficult to perform in small species such as S. ambulans (adults 3–6 mm, according to Sidorov and Palatov [12]).

1.4. The Aims of the Research

The general aim of the present research was to assess the relationship between the life history traits of the analysed species S. ambulans and their hyporheic habitat requirements. The specific aims of the study were to (i) estimate key population parameters of S. ambulans, including growth and mortality rates, using length-frequency distribution data; (ii) determine temporal and spatial differences in population dynamics and reproductive strategies (sex ratio and the number of cohorts) at two sampling sites in the HZ of the Sava River, Croatia; and (iii) analyse the possible influence of environmental variables on the measured population parameters. Two sampling sites close to each other with different environmental parameters in the HZ of the Sava River were selected to analyse the differences in life history traits and population dynamics according to previously published data on the metal bioaccumulation and energy status of S. ambulans [47,95]. This study provides a first insight into the population dynamics of S. ambulans inhabiting the HZ of the Sava River and a detailed complex overview of their life history traits in general.

2. Materials and Methods

2.1. Study Area

With a total length of 945 km, the Sava is the longest river in the Balkans and one of the most ecologically valuable rivers [96]. Due to the high biodiversity in the wetlands, six sites in the Sava River Basin have been declared Ramsar sites [97]. The Sava River Basin covers a total area of 97,713 km² [96], spanning the territories of six countries—Slovenia, Croatia, Bosnia and Herzegovina, Serbia, Montenegro, and Albania. It represents one of the most significant sub-basins of the Danube River [97]. Within the Sava River Basin, there are four ecoregions: Alps (ER 4), Dinaric Western Balkan (ER 5), Hungarian Lowlands (ER 11), and Hellenic Western Balkan (ER 6) [98]. High-quality groundwater is the primary drinking water supply in the Sava River basin, which includes two main units: the Pannonian Basin, with intergranular aquifers and the Dinarides, with predominantly limestone aquifers [97]. The Zagreb aquifer system is the only source of drinking water in the city of Zagreb. However, it is threatened by declining groundwater levels and increasing groundwater pollution [99]. In addition, the integrity of the upper and middle sections of the Sava River Basin is threatened by proposed new hydroelectric power plant projects [100,101]. Two study sites in northwestern Croatia, the Medsave and Jarun gravel bars, were selected to determine the biology and habitat ecology of S. ambulans in the HZ of the Sava River. The first sampling site is located in the village of Medsave (45°50′04″ N, 15°46′32″ E), in the rural zone, near the town of Zaprešić (Figure 1). The second sampling site, near the city district of Jarun (45°46′24″ N, 15°55′56″ E), is located in an urban area of the city of Zagreb (Figure 1). Since the Sava River is regulated, it was challenging to find gravel bars that were fully accessible during all seasonal sampling campaigns and where we would not have daily flooding of the gravel bars. The enormous impact of HEPP in Slovenia was evident in the extreme daily oscillations of the water level in the Sava River, which made the plans for collecting animals in the HZ unpredictable. The sediment structure at two sampling sites in the HZ of the Sava River was investigated using particle size distribution (PSD) analysis. Our results showed that sand was the prevailing fraction (66.2%), followed by silt (29.8%) and clay (4.0%). Previous biological studies of the HZ of the Sava River in Croatia were conducted only in the area around Zagreb, e.g., [15,102,103].

2.2. Sample Collection and Field Measurements

Sampling was conducted once per season in December (winter) 2018 and April (spring), July (summer), and October (autumn) 2019 at each of the study sites located in the HZ of gravel bars (GB) on the Sava River. Sampling sites were selected along the 50–100 m GB. Because hyporheic environmental conditions and the invertebrate community in the river sediments are highly variable and depend on the hydrologic patterns, at least three replicates were collected at each site. The position of each replicate was determined by the distance from the watercourse in meters. For each sampling replicate, 50 L of interstitial water mixed with sediment and extracted aquatic invertebrates was collected from various depths of at least 40 cm using the Bou-Rouch method, using a piston pump attached to a movable steel tube (Ø 50 mm, 110 cm long) with a perforated distal end (35 holes 5 mm in diameter) [104]. The sampling of interstitial fauna from the extracted mixture was performed using a set of fine-mesh hand nets (mesh sizes of 100, 200, and 500 µm). The sampled organisms were stored in plastic bottles and fixed with 80% ethanol. In the laboratory, specimens of S. ambulans were identified using a stereomicroscope (Zeiss Stemi 2000 C, Jena, Germany) based on the morphological characteristics of the species described in Karaman [105], sorted by gender and ontogeny (adult, adolescent, juvenile) and counted. In addition to S. ambulans, other taxa were identified in small numbers, for example, other amphipod species from the genus Niphargus spp., Copepoda, Ostracoda, Turbellaria, larvae of Ephemeroptera, larval and adult stages of Coleoptera, dipteran larvae Chironomidae, and bivalves of the family Sphaeriidae and Oligochaeta. After sampling, physico-chemical water properties (water temperature (°C WT), pH, conductivity (μS cm⁻¹ Cond), dissolved oxygen concentration (mg L⁻¹ DO), and oxygen saturation (% Sat)) were measured in pumped water at each study site at each subsampling using the respective portable field meters, which were calibrated before samplings (oximeter OXI 96, WTW GmbH, Weilheim, Germany, for WT and DO; pH-meter 330i, WTW GmbH, Weilheim, Germany, for pH and conductometer Sension 5; Hach, Loveland, CO, USA, for Cond). Alkalinity (mg L⁻¹ CaCO₃ Alk) was also measured in situ via the acid titration (0.1 mol L⁻¹ HCl) of water samples using methyl orange as a colour indicator. For the subsequent water analysis in the laboratory, an additional 1 L sample of the interstitial water was collected at each site.

2.3. Laboratory Measurements

2.3.1. Determination of Environmental Parameters

One litre of hyporheic water was pumped from the HZ for laboratory physicochemical analyses of eight parameters. Total water hardness (mg CaCO₃ L⁻¹ TWH), concentrations of nitrite (mg L⁻¹ N-NO₂⁻), nitrate (mg L⁻¹ N-NO₃⁻), and orthophosphate (μg L⁻¹ P-PO₄³⁻) were determined using the standard APHA methods [106], while total chemical oxygen demand (mg O₂ L⁻¹ COD_KMnO4) was measured using the potassium permanganate method [107]. Dissolved organic carbon (mg L⁻¹ DOC) and particulate organic carbon (mg L⁻¹ POC) were measured in the interstitial water, while total organic carbon (% TOC) was measured in the sediment samples. DOC concentrations were determined in triplicate using the sensitive High-Temperature Catalytic Oxidation (HTCO) method at 680 °C [108,109]. A total organic carbon analyser TOC-V_CPH (Shimadzu, Kyoto, Japan) with a Pt/Silica catalyst (Elemental microanalysis, Okehampton, UK) and Non-Dispersive Infrared (NDIR) detector for CO₂ were used for DOC measurements. The average instrument blank and Milli-Q blank were equivalent to 0.07 mg C L⁻¹ (n = 30); the reproducibility of measurements was high (1.9%). POC was analysed with an SSM-5000A solid sample module associated with a TOC-V_CPH carbon analyser, calibrated with glucose. POC concentrations were corrected on the basis of blank filter measurements. The average filter blank, including the instrument blank, was equivalent to 0.005 mg C L⁻¹. The reproducibility obtained for the glucose standard was high (3.0%). TOC content was determined from lyophilised, well-homogenised sediment samples using an SSM-5000A module connected to a TOC-V_CPH analyser through a high temperature (900 °C) catalytic (mixture of Pt/silica and CoO) oxidation method with IR detection of CO₂ and calibrated with D (+) glucose (Merck, Darmstadt, Germany) [110]. The methods for the determination of organic carbon (DOC, POC, TOC) are accredited according to HRN EN ISO/IEC 17025:2017.

2.3.2. Amphipod Measurements

Each amphipod was measured with a stereomicroscope (Zeiss Stemi 2000 C, Jena, Germany) using an ocular micrometre with the accuracy of 0.1 mm. Total body length (TBL) was measured from the anterior margin of the head to the anterior margin of the telson. Gender was identified by the presence of penes on the 7th sternite of the male pereon and the presence of oostegites in females. Adults were defined as individuals > 3 mm with wholly developed penes or oostegites, while adolescents (Ad. Male and Ad. Female) were 2–3 mm in TBL without wholly developed penes or oostegites. The ovigerous females had fully developed, setose oostegites and eggs in the brood pouch, while postovigerous females (Post. Female) had fully developed marginal setae but an empty brood pouch. Individuals without penes or oostegites less than 2 mm in TBL were classified as juveniles (Figure 2). The sex ratio (number of females/number of males) was also calculated seasonally at each study site.

2.4. Data Analysis

Data were tested for normality (Shapiro–Wilk test) and variance homogeneity. Because normality and variance homogeneity were not met, non-parametric statistics were chosen. Differences between the mean TBL of ontogeny classes (female, female adolescent, male, male adolescent, juvenile) among sampling sites were tested using a Mann–Whitney U test. Ovigerous females were excluded from analyses due to their small number (in total, 3 individuals). The significance of variability in TBL between sampling seasons was tested using Kruskal–Wallis ANOVA by ranks, followed by Fisher LSD post hoc tests, separately for each sampling site. The chi-square test (χ²-test) with Yeats correction was used for testing whether the sex ratio differed significantly from the expected equilibrium 1:1. All statistical analyses were performed using the Statistica software package (v.14.0.0.15, TIBCO Software Inc., Palo Alto, CA, USA; 2020) and values of p < 0.05 were considered significant.

2.4.1. Cohort and Growth Analyses

Specimens were pooled into 0.25 mm length classes (starting at 0.90 mm and ending at 6.38 mm) for each sampling season and study site. The VBGF (Von Bertalanffy Growth Function) was used to estimate population growth parameters in the FAO-ICLARM Stock Assessment Tools II (FISAT II software, ver. 1.2.2, provided by the Food and Agriculture Organization, http://www.fao.org/fishery/topic/16072/en, accessed on 27 January 2022). Growth parameters of each population were estimated using the Electronic Length Frequency Analysis I (ELEFAN I) routine integrated into the FISAT II [111]. ELEFAN I was used to determine the growth coefficient K with units of reciprocal time (year⁻¹) and asymptotic length (L_∞). The K scan routine was performed to obtain a reliable estimate of the K value. Length-frequency distributions for each population, sex, and study site (Jarun and Medsave) were used to detect and separate cohorts using the Bhattacharya method [112] in FISAT II. The mean total body length (TBL), standard deviation, and number of individuals were estimated using the Bhattacharya method and NORMSEP analyses of Modal Progression Analysis (MPA) in FISAT II. Modal separation for the study site Jarun and Medsave with a Separation Index (SI) greater than 2 indicated a single cohort. Each cohort was identified seasonally and recorded with length-frequency histograms. Seasonal size-frequency distributions at both study sites were compared. The hypothetical age at zero size (t₀) was calculated according to the standard equation by Pauly [113].

The growth performance index (φ′) was calculated as follows: φ′ = 2 log₁₀(L_∞) + log₁₀(K) [114].

2.4.2. Mortality and Longevity

Total mortality (Z, years⁻¹) was calculated based on the seasonal length converted catch curve (FISAT routine) using linear regression [115,116]. Natural mortality (M, year⁻¹) was calculated using the following equation [115]:

log(M) = −0.0066 − 0.279 log(L_∞) + 0.6543 log(K) + 0.4634 log(T)

where L_∞ is the asymptotic length, K is the growth coefficient of the VBGF, and T is the mean annual water temperature of the habitat. We recorded 14.99 °C and 15.10 °C as the mean annual water temperature in the hyporheic interstitial for the Jarun and Medsave study sites, respectively, where the shrimps were collected during the research.

The estimated potential longevity (t_max) of the species was calculated using the relationship t_max = 3/K [115].

Growth, mortality, and longevity were estimated separately for males and females and the combined sex.

2.4.3. Principal Component Analysis (PCA) and Canonical Correspondence Analysis (CCA)

Prior to the analyses of environmental parameters, all variables were normalized. Posteriorly, a resemblance matrix based on the Euclidean distance was calculated. The PCA was performed in order to detect seasonal sample differences among study sites based on environmental data. The PCA analyses were performed using PRIMER v6 [117].

CCA was performed to determine the relationship between the abundance of each gender (male and female) and ontogeny classes (female, female adolescent, male, male adolescent, juvenile) of S. ambulans and the measured environmental variables using CANOCO 5 software [118]. In the ordination figures, environmental variables are represented as arrows, and the arrow’s length indicates each variable’s relative importance [119]. The angle between the arrows indicates the degree of correlation among environmental variables, and the arrangement of gender and ontogeny classes concerning the arrow represents their responses to the specificity of the environmental preferences. The environmental variables used in CCA were WT, pH, Cond, DO, Sat, TWH, and concentrations of N-NO₂⁻, N-NO₃⁻, P-PO₄³⁻, COD, DOC, POC, TOC, and Alk. The sampling depth and position of the Bou-Rouch pump were also added to the analysis. The Monte Carlo randomisation test (499 permutations under the reduced model) was performed to determine the statistical significance of the correlations [119].

3. Results

3.1. Population Structure

A total of 616 specimens (adults and adolescents) of S. ambulans (185 female, 431 male) from the Jarun study site and 892 (196 female, 696 male) from the Medsave study site were analysed (

Table 1. Number of S. ambulans specimens per season, total body length (TBL) range (min–max), mean TBL ± SD for females and males, and sex ratio at the Jarun site (* indicate values significantly different from 1:1, χ²-test, p < 0.05). Adolescent females and males are included.

	Total Body Length (mm)								Sex Ratio (F:M)
	Females				Males
	n	Min–Max	Mean	±SD	n	Min–Max	Mean	±SD
Winter	26	3.30–6.00	4.39	0.50	65	2.65–4.15	3.38	0.32	0.40:1 *
Spring	40	2.80–6.00	5.16	1.04	118	2.40–4.50	3.70	0.38	0.34:1 *
Summer	15	2.80–4.10	3.26	0.49	131	1.75–4.50	2.78	0.55	0.11:1 *
Autumn	104	2.80–5.10	3.97	0.54	117	2.30–4.00	3.28	0.40	0.89:1
Total	185	2.80–6.00	4.20	0.79	431	1.75–4.50	3.29	0.38	0.43:1 *

). A total of 123 juveniles were sampled at the Jarun site and 38 juveniles were sampled at the Medsave site. At the Jarun sampling site, the most distinct population peaks occurred in the spring (April 2019) and autumn (October 2019). At the Medsave sampling site, the highest population abundances were observed in summer (July 2019) and autumn (October 2019) (Table 1). There was no statistically significant difference in the S. ambulans abundances between the two sampling sites (p = 0.221).

The frequency distribution of the different life stages is shown in Figure 3. At the Jarun site, the peak of juveniles represented 26% and 29% of the population in spring and summer, respectively. The highest peak of adolescents was in the summer and accounted for 71% of the population. The two peaks of adults were in winter and spring and accounted for 87% and 67% of the population, respectively. At the Medsave site, an increase in juveniles was observed in the summer. The highest peak of adolescents was in the summer and accounted for 91% of the population. Adults accounted for more than 90% of the population in winter and summer. Intersexuality was not detected in any of the amphipod populations in this study.

The statistical comparison of TBL between study sites showed that the TBL of females and males in winter and spring and that of males in autumn were significantly higher at the Medsave site compared to the Jarun site (Table 1 and

Table 2. Number of S. ambulans specimens per season, total body length (TBL) range (min–max), mean TBL ± SD for females and males, and sex ratio at the Medsave site (* indicate values significantly different from 1:1, χ²-test, p < 0.05). Adolescent females and males are included.

	Total Body Length (mm)								Sex Ratio (F:M)
	Females				Males
	n	Min–Max	Mean	±SD	n	Min–Max	Mean	±SD
Winter	20	4.13–5.78	4.92	0.50	83	3.00–5.00	3.77	0.34	0.24:1 *
Spring	29	4.20–6.38	5.43	0.41	195	2.75–4.80	3.84	0.35	0.15:1 *
Summer	85	2.10–4.50	3.03	0.59	179	1.80–4.62	2.94	0.63	0.47:1 *
Autumn	62	2.75–4.90	3.83	0.54	239	1.80–4.40	2.99	0.48	0.26:1 *
Total	196	2.10–6.38	4.30	1.08	696	1.80–5.00	3.39	0.49	0.28:1 *

). A significantly higher TBL of adolescent males was found only in autumn at the Jarun site (Mann–Whitney U test, level of significance p < 0.05). A comparison between TBLs of males at both study sites revealed no statistically significant difference, whereas females at the Jarun site had significantly higher TBLs compared with females at the Medsave site (p < 0.001). During all four seasons, females had significantly higher TBL than males at both study sites (p < 0.001). Significant seasonal differences were found at the Medsave site for females, female adolescents, males, and male adolescents, whereas at the Jarun site differences were found for females, female adolescents, males, and juveniles (Kruskal–Wallis One Way Analysis of Variance on Ranks, p < 0.05).

3.2. Length-Frequency Analysis

For the length-frequency analysis, female adolescents were pooled with adult females and male adolescents with adult males. The frequency distribution of female TBL ranged from 2.80 to 6.00 mm (mean female TBL = 4.23 mm; SD = 0.87) in the Jarun study site, whereas it ranged from 2.10 to 6.38 mm (mean female TBL = 3.84 mm; SD = 1.04) in the Medsave study site in all seasons (Table 1 and Table 2). Male TBL ranged from 1.75 to 4.50 mm (mean male TBL = 3.26 mm; SD = 0.56) in the Jarun study site, and from 1.80 to 5.00 mm (mean male TBL = 3.31 mm; SD = 0.64) in Medsave (Table 1 and Table 2). At both study sites and in all four seasons, females had a greater size range (2.10–6.38) than males (1.75–5.00). At the Jarun site, females had a wider range of size classes in winter, spring, and autumn. Males had a wider range of size classes during the summer (Figure 4a). At the Medsave site, females had a wider range of size classes in spring, whereas males had a wider range of size classes in winter, summer, and autumn (Figure 4b).

The highest frequencies of female TBL were between 3.50 and 3.80 mm in Jarun in autumn and between 2.60 and 2.80 mm, and between 3.10 and 3.30 mm in Medsave in summer, while the highest TBL frequencies of males were recorded between 3.50 and 3.80 mm in Jarun in spring and between 2.50 and 2.80 mm in Medsave in autumn. Overall, females were more abundant in the smaller TBL size classes, while males were more abundant in the larger TBL size classes. The values of sex ratio on the Jarun study site were significantly different from 1:1 in most seasons (χ²-test, p < 0.05), and they were male based, except in autumn, when the population was in equilibrium (Table 1). The sex ratio was significantly different from the expected equilibrium in all seasons in the Medsave study site, always favouring males (Table 1 and Table 2). The sex ratio was exceptionally biased towards males in spring at Medsave and in summer at the Jarun site. The breeding period could not be estimated because ovigerous females were recorded in extremely low numbers (only three specimens in December).

The estimated parameters of VBGFs for Jarun and Medsave sites differed between sexes (

Table 3. Growth parameters of the von Bertalanffy growth function of S. ambulans during the period from winter (December) 2018 to autumn (October) 2019.

	Jarun			Medsave
Growth Parameters	Females (n = 185)	Males (n = 431)	Combined Sex (n = 616)	Females (n = 196)	Males (n = 696)	Combined Sex (n = 892)
Growth coefficient (K) (year−1)	0.92	0.59	1.50	0.34	1.20	0.50
Hypothetical age at zero size (t0) (year)	−0.57	−0.36	−0.75	−0.14	−0.64	−0.31
Asymptotic length (L∞) (mm)	6.30	5.40	4.40	7.20	4.00	6.60
Longevity (tmax) (year)	3.26	5.08	2.00	8.82	2.50	6.00
Growth index φ′	1.56	1.23	1.46	1.25	1.28	1.34

). For females, K and φ′ were higher at the Jarun site than for females at the Medsave site, whereas L_∞ and t_max were higher at the Medsave site. For males, L_∞ and t_max were higher at the Jarun site than for males at the Medsave site, while K and φ′ were higher at the Medsave site. At the Jarun study site, females had higher values of K, L_∞, and φ′ than males, while males had higher t_max compared to females. At the Medsave study site, females had higher L_∞ and t_max, while males had higher K and φ′ compared to females. The estimated growth parameters of different species of Amphipoda obtained from other studies are compared in

Table 4. Comparison of growth parameters L_∞ (asymptotic length), K (growth coefficient), natural mortality coefficient (M) and/or total mortality (Z) of S. ambulans with different species of Amphipoda. The data are presented alphabetically, by the species name.

Species	Sex	Mean Size (mm)	L∞ (mm)	K (year−1)	t0 (year)	M or Z (year−1)	Sampling Period	Location	Reference
Gammarus wilkitzkii **	Females	19.5–49.7	64.6	0.48			Summer 1994	Franz Josef Land, Arctic	[120]
	Males	62.3 (max) (min–max)	69.2	0.47
Nototropis minikoi **	Females	5.12 ± 0.59					June 2015	Campechen Lagoon, Quintana Roo, Caribbean Sea, Mexico	[82]
	Males	4.04 ± 0.42
	Sex combined	0.94–6.34 (min–max)	6.66	0.67	0.18	/
Pseudoichestoidea brasiliensis **	Females	7.3–12.2 (min–max)					May 1996–December 1997	Arachani beach, Uruguay	[80]
	Males
	Sex combined	8.74 3–12.2 (min–max)	11.64	1.61	−0.175	2.35 (M)
	Females	7.8–12.1 (min–max)						Barra del Chuy beach, Uruguay
	Males
	Sex combined	8.92 3–12.2 (min–max)	11.27	1.88	−0.175	2.47 (M)
Synurella ambulans *	Females	4.23 ± 0.87	6.75 +	0.63	−0.36	0.86 (M)	Winter 2018–autumn 2019	Sava River, Croatia	Present study
	Males	3.26 ± 0.56	4.70 +	0.90	−0.50	2.70 (M)
	Sex combined	3.28 ± 0.96	5.50 +	1.00	−0.53	1.25 (M)
Waldeckia obesa **	Females	5–31 (min–max) sex combined	7.47	0.50		(0.27/year) (Z)	January–March 1998	eastern Weddell Sea, Antarctic	[66]
	Males		6.92	0.60		(0.43/year) (Z)

Note: * freshwater, ** marine, ⁺ all values are means of two sampling sites.

.

3.3. Cohorts and Growth

The four cohorts of S. ambulans were identified during one year at each study site (Jarun and Medsave) based on length-frequency distribution and confirmed with the Bhattacharya method (Figure 5, Figure 6 and Figure 7). At the Jarun study site, relatively large overwintering adults (Cohort I and Cohort II) were detected from December 2018 to late April 2019 (Figure 5 and Figure 7a). Adults of Cohort III were detected from July to October 2019, and adults of Cohort IV from April to October 2019. The largest amphipods disappeared from the population after the spring season (Figure 7a). The largest adult males and females were found in the population in April 2019, with a mean total body length of 4.40 mm and 6.00 mm, respectively.

At the Medsave study site, amphipods at the upper level consisted of two parallel cohorts (Cohorts I and II) from December 2018 to April 2019 (Figure 6 and Figure 7b). During the spring and early summer (April to July 2019), a new cohort, Cohort III, appeared in the population (Figure 7b). Cohort IV increased in size and number from April to October 2019. Juveniles from Cohort III and Cohort IV probably appeared in the population in the early summer of 2018 and autumn 2018, respectively. The largest adult males were detected in the population in December 2018, with a mean total body length of 4.75 mm, and disappeared from the population by summer 2019. The largest adult females were recorded in April 2019, with a mean total body length of 5.70 mm, and disappeared from the population during the summer. The longest cohort life span of both populations was eight months (Figure 7a,b).

3.4. Mortality

In the present study, amphipod mortality was calculated for the entire population of S. ambulans (

Table 5. Mortality parameters of S. ambulans in the hyporheic zone of two study sites (Jarun and Medsave) at the Sava River.

	Jarun			Medsave
Parameters	Females	Males	Combined Sex	Females	Males	Combined Sex
Total mortality (Z, year−1)	3.13	4.60	3.48	1.77	1.92	2.10
Natural mortality (M, year−1)	1.96	1.53	2.98	0.99	2.65	1.30
M/Z (%)	62.62	33.26	85.63	55.93	138.02	61.90
M/K	2.13	2.59	1.99	2.91	2.21	2.60
Z/K	3.40	7.80	2.32	5.21	1.60	4.20
Annual mean temp. (°C)	14.99			15.10

). The total mortality coefficient (Z) was estimated from the slope value of the length-converted catch curve, separately for the combined sex, females, and males of the Jarun and Medsave populations. The estimated total mortality (Z) for separate sexes had higher values for males than females at both sampling sites. The estimated natural mortality (M) was higher for males compared to females at the Medsave site and higher for females at the Jarun site. The highest value of natural mortality (2.65 year⁻¹) and total mortality (4.60 year⁻¹) was calculated for males at the Medsave and the Jarun study sites, respectively. A comparison of length-converted catch curves showed that M was higher for females at the Jarun site, and for males at the Medsave site. Total mortality was higher for both females and males at the Jarun site than at the Medsave site. The natural mortality of the Jarun population (M = 2.98 year⁻¹, combined sex) corresponded to about 86% of Z, revealing the low propagation of this species in the HZ of the Sava River at the Jarun study site. The natural mortality of the Medsave population (M = 1.30 year⁻¹, combined sex) corresponded to about 62% of Z, revealing the extremely low repopulation of the species in the HZ in this part of the Sava River. The natural mortality to growth factor (M/K) ratio was highest at the Medsave site for females and lowest at the Jarun site for the combined sex (males and females combined).

Both the M/K ratio and the Z/K ratio were higher for females than for males at the Medsave site, and for males at the Jarun site.

3.5. Relationship between Environmental Parameters and Gender/Ontogeny Classes’ Abundance

A previous paper already presented the detailed environmental parameters [47]. In this research, new parameters (sampling depth and position of the Bou-Rouch pump) were added that were the key factors in the analysis of S. ambulans populations in the HZ, and they were all combined using the new types of analyses (PCA and CCA). A seasonal temperature gradient was observed in the HZ of the Sava River, ranging from 10.4 °C in winter (December) to 23.3 °C in summer (July) at the Medsave site and from 12.7 °C in winter (December) to 16.5 °C in summer (July) at the Jarun site (

Table 6. Mean, standard deviation (SD), minimum (min), and maximum (max) of the environmental (water) parameters in the HZ of two study sites, Jarun and Medsave, on the Sava River.

Study Site	Jarun				Medsave
Season	Winter	Spring	Summer	Autumn	Winter	Spring	Summer	Autumn
Position (m)	0.22 ± 0.73 −0.30–3.00	0.83 ± 1.35 0.00–3.00	0.14 ± 0.07 0.00–0.20	1.86 ± 2.05 0.10–5.50	1.34 ± 2.29 0.30–8.00	0.16 ± 0.69 −0.30–2.00	0.00 ± 0.00 0.00–0.00	0.10 ± 0.00 0.10–0.10
Depth (cm)	65.28 ± 4.23 57–70	56.88 ± 4.32 50–65	51.50 ± 2.07 49–55	56.38 ± 8.75 45–65	65.43 ± 1.34 62–67	54.00 ± 1.15 52–56	51.25 ± 0.96 50–52	48.67 ± 3.21 45–51
WT (°C)	12.97 ± 0.77 12.30–13.80	12.66 ± 0.06 12.55–12.75	16.35 ± 0.16 16.00–16.50	16.56 ± 0.76 16.10–18.40	10.84 ± 0.94 9.50–11.50	12.14 ± 0.04 12.10–12.20	23.30 ± 0.00 23.30.−23.30	14.50 ± 0.00 14.50–14.50
DO (mg O2 L−1)	6.09 ± 0.94 5.07–6.91	5.72 ± 0.16 5.45–5.88	4.37 ± 0.10 4.20–4.50	6.45 ± 1.46 4.40–7.59	8.41 ± 1.17 7.64–10.18	5.73 ± 0.15 5.50–5.90	6.04 ± 0.14 5.90–6.20	8.94 ± 0.02 8.92–8.95
Sat (%)	44.63 ± 24.24 18.30–65.70	56.35 ± 1.38 54.60 ± 58.60	44.50 ± 1.41 43–47	65.90–14.43 45–77	77.24 ± 8.06 72.20–89.50	56.59 ± 1.21 54.90–58.60	6.04 ± 0.14 5.90–6.20	8.94 ± 0.02 8.92–8.95
pH	8.05 ± 0.07 7.97–8.11	7.45 ± 0.04 7.39–7.52	7.42 ± 0.05 7.30–7.46	7.80 ± 0.24 7.51–8.00	8.30 ± 0.09 8.24–8.44	8.33 ± 0.04 8.25–8.37	8.51 ± 0.08 8.40–8.58	7.96 ± 0.06 7.92–8.03
Cond (µS cm−1)	541 ± 27.61 517–571	578.75 ± 2.31 575–580	584.88 ± 1.55 583–587	576.25 ± 42.25 522–634	447.86 ± 39.57 388–488	420 ± 2.38 418–425	431 ± 1.29 430–433	439 ± 1.73 437–440
Alk (mg CaCO3 L−1)	261.67 ± 15.34 245–275	278.13 ± 3.72 270–280	212.50 ± 8.02 200–225	279.38 ± 10.92 265–290	242.86 ± 11.72 225–250	157.25 ± 0.26 157–157.50	201.75 ± 2.36 200–205	215 ± 25.98 200–245
TWH (mg CaCO3 L−1)	225.27 ± 7.28 218.94–233.18	277.68 ± 0.00	270.56 ± 0.00	235 ± 33.38 202.92–267.00	242.08 ± 0.00	243.86 ± 0.00	137.06 ± 0.00	267.00 ± 0.00
P-PO43− (mg L−1)	0.08 ± 0.02 0.06–0.10	0.24 ± 0.00	0.17 ± 0.00	0.16 ± 0.05 0.08–0.20	0.06 ± 0.00	0.10 ± 0.00	0.04 ± 0.00	0.10 ± 0.00
N-NO2− (mg L−1)	0.12 ± 0.02 0.10–0.13	0.18 ± 0.00	0.07 ± 0.00	0.14 ± 0.01 0.12–0.15	0.14 ± 0.00	0.06 ± 0.00	0.09 ± 0.00	0.16 ± 0.00
N-NO3− (mg L−1)	0.69 ± 0.10 0.60–0.80	1.96 ± 0.00	0.24 ± 0.00	1.53 ± 0.86 0.65–2.34	0.46 ± 0.00	2.31 ± 0.00	1.48 ± 0.00	2.72 ± 0.00
COD (mg O2 L−1)	1.76 ± 0.08 1.69–1.85	5.34 ± 0.00	2.55 ± 0.00	2.63 ± 0.97 1.69–3.54	1.89 ± 0.00	3.77 ± 0.00	5.07 ± 0.00	4.48 ± 0.00
DOC (mg L−1)	0.63 ± 0.02 0.60–0.65	0.82 ± 0.00	0.70 ± 0.00	0.90 ± 0.00	1.36 ± 0.00	1.37 ± 0.00	1.37 ± 0.00	0.94 ± 0.00
POC (mg L−1)	0.71 ± 0.15 0.55–0.84	0.59 ± 0.00	0.91 ± 0.00	2.63 ± 0.00	0.13 ± 0.00	0.37 ± 0.20 0.13–0.52	0.56 ± 0.05 0.52–0.62	0.69 ± 0.00
TOC (%)	2.32 ± 0.00	1.77 ± 0.00	1.92 ± 0.00	1.73 ± 0.00	1.07 ± 0.00	2.68 ± 0.00	2.97 ± 0.00	2.83 ± 0.00

).

The PCA analysis indicated that the first two axes of the PCA explained 48.9% of the variance among samples (PC1 = 27.2% and PC2 = 21.7%) and revealed slight differences between environmental parameters on the study sites, mainly due to Alk, Cond, depth, DO, DOC, pH, and P-PO₄³⁻. The pH (eigenvalue = −0.432) and Cond (eigenvalue = 0.406) presented the strongest correlations in PC1, while Alk (eigenvalue = 0.404) and N-NO₃⁻ (eigenvalue = −0.323) were the most important variables in explaining the PC2 (

Table 7. Statistical summary of axes generated by PCA for environmental parameters.

	PC1	PC2
Eigen values	4.63	3.69
% of variance	27.2	21.7
Position	−0.067	0.206
Depth	−0.163	0.366
WT	0.164	−0.26
DO	−0.31	0.193
Sat	−0.226	0.073
pH	−0.432	−0.097
Cond	0.406	0.191
Alk	0.156	0.404
TWH	0.209	0.114
P-PO43−	0.381	0.008
N-NO2−	0.077	0.313
N-NO3−	0.076	−0.323
COD	0.136	−0.31
DOC	−0.332	−0.159
POC	0.27	0
TOC	0.044	−0.402

, Figure 8).

The influence of environmental parameters on the amphipod analysed classes was investigated using CCA, which showed that overall, different classes of S. ambulans (gender (male and female) and ontogeny (adult, adolescent and juvenile)) were associated differentially with environmental variables. All classified categories were distributed along the four axes of the CCA. The first axis had an eigenvalue of 0.3 and the second axis had an eigenvalue of 0.1, with both axes exhibiting gender/ontogeny–environment correlations of 0.9. The Monte Carlo test was significant (test of significance of all canonical axes: trace = 0.56; F-ratio = 9.9; p < 0.005) and cumulatively, axes 1 and 2 accounted for 77.6% of the total variance, with correlations between the density of S. ambulans classes (gender and ontogeny) and environmental variables of 0.973 (axis 1) and 0.894 (axis 2). According to the correlation coefficients, DO, DOC, position, and Sat had a high positive correlation to the first CCA axis, while WT and POC were negatively correlated. The depth of samples, TWH, and concentration of P-PO₄³⁻ were positively correlated with the second CCA axis, while Cond and COD were negatively correlated. Adult males were primarily associated with the concentration of N-NO₂⁻, while adult females were associated preferably with the depth of samples and TWH. The juveniles were associated with Cond, while adolescent females were primarily associated with POC (Figure 9).

The ordination of the samples showed that they tended to cluster following the seasons, revealing a temporal gradient according to the CCA axes (Figure 10). Winter samples (W) with higher DO and Sat and pH value clustered on the positive side of the first CCA axis. Spring samples (SP) with high Cond and P-PO₄³⁻ clustered on the negative side of the first CCA axis. Summer samples (SU) with high Cond and the lowest concentration of DO clustered mainly on the negative side of the second CCA axis. Autumn samples (A) were characterised with the highest TOC concentration and clustered on the negative side of the first CCA axis and the positive side of the second CCA axis (Figure 10).

4. Discussion

4.1. Life Cycle

The life cycle of S. ambulans has previously been studied in populations from Romania, Poland, and Germany, respectively [7,11,33]. These are the only data on the life history of S. ambulans in the literature. The population structure of S. ambulans from Poland and Romania differed from the population in the Sava River. In our study, the recruitment of the new generation began in April, while in the population in Romania it began in May and in Poland it began later, in June [7,11]. Thus, the population in Romania became bimodal (consisting of adults and juveniles) in May, the population in Poland became bimodal in July and August, and the population in the Sava River became bimodal in April. In all populations (Romania, Poland and Croatia), the parent generation had disappeared in summer (June, July, and August). The juveniles that appeared in Poland in summer grew to sexual maturity in the following spring. It was concluded that S. ambulans in Romania have an annual univoltine (one generation per year) semelparous (single-breeding) life cycle [7]. In the study from Poland, it was suggested that the life cycle of S. ambulans is univoltine and iteroparous (females produce only a few broods in their lifetime) [11]. It has been shown that gammarids living in stable water temperatures with a groundwater supply are iteroparous and breed continuously [121]. Our study showed that the life cycle of S. ambulans may be univoltine, with the extremely low number of ovigerous females found only in December (three specimens). Such differences in life cycle between S. ambulans populations at different latitudes or altitudes (Poland, Romania, and Croatia) may be due to seasonal variations (temperature regime, different longevity of unfavourable conditions) [7]. In Poland, the water temperature ranged from 0.8 °C in January to 18.8 °C in July. In Romania, the temperature range was the highest, with water temperature ranging from 2.5 °C in December to 24 °C in July and August. In our study, the temperature ranged from 10.4 °C in December to 23.3 °C in July (Table 6).

Thus, different temperature minima and maxima may result in different recruitment periods among these S. ambulans populations. This was particularly evident in the later reproduction and recruitment in the population from Poland, where temperatures were lowest and the latitude highest. Depending on whether they live at high or low latitudes, gammarids are semelparous or iteroparous, respectively [56]. In the amphipods Monocorophium insidiosum (Crawford, 1937) and Corophium multisetosum Stock, 1952, the lifespan and number of generations were also affected by latitude [122,123]. The number of cohorts of the same crustacean species, Neomysis americana (S.I. Smith, 1873), was shown to be more strongly related to temperature regime (the higher number of cohorts was recorded at higher temperatures) and to a lesser extent to latitude [84]. In addition, the S. ambulans populations from Poland and Romania were epigean—found in a shallow pond and flowing rivers in lowlands, respectively. On the other hand, the S. ambulans population in Croatia was hypogean—inhabiting the river’s HZ (GWCE). Thus, the different habitats could also be the reason for the differences in the life cycle of S. ambulans.

Data on the sex ratio of S. ambulans indicate that males were more abundant than females in all seasons, while the opposite was found in populations from Poland and Romania [7,11]. The overall sex ratio in Poland and Romania was 1:0.4 (634 females and 247 males in Poland), while in our study it was 0.4:1 (185 females and 431 males) at the Jarun site and 0.3:1 (196 females and 696 males) at the Medsave site. Compared to other amphipod species, such as Crangonyx pseudogracilis Bousfield, 1958 and Waldeckia obesa (Chevreux, 1906), the sex ratio was in favour of females [66,83], while the proportion of males was higher than that of females in Nototropis minikoi (A.O. Walker, 1905) [82]. Several factors may influence the changes in sex ratio, such as different mortality rates between males and females [88,124,125,126], different longevity or differences in the minimum size at which secondary sexual characteristics appear [125,127]. One of the reasons for the differences in sex ratios in amphipods could be chemicals from agriculture and from untreated sewage [83,128]. In our study, the male-biased ratio could be due to the females (which are larger than males) being more vulnerable to the lower oxygen concentration in the HZ and requiring larger amounts of food, and therefore inhabiting the shallower part of the HZ with more oxygen and food. Moreover, it is possible that the male-biased ratio is due to the fact that males require a greater number of moults to become sexually mature. Synurella ambulans has a longer life cycle compared to other epigean amphipods [83,129], resulting in more extended maturation and growth. It has been previously shown that amphipods disappear from the population after mating (males) or after the release of the brood (females) [130]. Thus, males of S. ambulans remain in the population in high numbers until females are sexually mature and until they have completed breeding.

The observed TBL values of S. ambulans were similar to the TBL values reported in previous studies [11]. The maximum size of males (4.7 mm) was comparable to the maximum size of 4.8 mm recorded in our study, whereas the maximum size of females (7.0 mm) was higher than in our study (6.4 mm). The size difference between females and males might be related to mating behaviour, since there is no pre-copula pairing (where males are larger than females) in S. ambulans [7]. In our study, the maximum TBL of females and males were higher at the Medsave site than at the Jarun site. It was suggested that the significant differences in the length of individuals of Dikerogammarus haemobaphes (Eichwald, 1841) (Amphipoda) inhabiting two sites with different hydrological characteristics could be due to different hydrological conditions and different predation pressures [81]. In our study, the water levels and discharges of the Sava River were similar at the two sampling sites. Another important factor influencing the differences in maximum TBL values between sampling sites could be nutrient level. DOC and TOC concentrations were higher at the Medsave site than at the Jarun site, which could result in higher maximum TBL values. Similarly, it was shown that the largest size of Emerita analoga (Stimpson, 1857) females was associated with food availability [131].

Regarding seasonality, the mean TBL of females and males at both sites was lowest in summer because the population consisted only of juveniles and adolescents. When we observed only sexually mature individuals, the lowest mean TBL was in autumn for most age groups (females at the Medsave site and males at the Medsave and Jarun sites). The opposite pattern was observed in the study on adult individuals of the amphipod Dikerogammarus haemobaphes, who had the lowest size in summer and the highest size in autumn [81].

4.2. Growth and Mortality

The analysis of growth parameters revealed that both males and females had different growth rates depending on location. Recent studies have shown that population density can affect crustacean growth, as they were negatively correlated [86,132,133]. However, in our study, no significant correlation was found between amphipod abundance (n of individuals in 50 L) and mean animal size (Pearson coefficient r_female = −0.425, r_male = −0.307; p > 0.05). The estimated growth parameters (K and φ′) for the two populations (Jarun and Medsave) showed that females grew faster than males at the Jarun site, while the opposite was observed at the Medsave site. The females of Crangonyx richmondensis Ellis, 1941, C. pseudogracilis (both belonging to the same family of Crangonyctidae as S. ambulans) and Dikerogammarus villosus (Sowinsky, 1894) had a higher growth rate than males [83,134,135]. In contrast, a higher growth rate was recorded for males than females in the Antarctic amphipod W. obesa [66].

In the populations of S. ambulans from the Jarun and Medsave site, four cohorts were identified during one year, which is an extremely low number of generations in comparison to small surface river dwelling amphipods such as Echinogammarus cari (S. Karaman, 1931), with as many as 12 cohorts per year [136]. At the Jarun study site, the large overwintering adults belonging to Cohort I and Cohort II would produce the new offspring, which would appear in the population at the beginning of the next year (2020). Since adults of Cohort III and of Cohort IV were detected in the population in spring, summer, and autumn 2019, juveniles from Cohort III and Cohort IV likely appeared in the population in early summer 2018 and autumn 2018, respectively. At the Medsave study site, Cohorts I and II would produce new offspring, which would appear in the population at the beginning of the next year (2020). Based on the observations of Cohort IV, which gradually increased in size and number from April to October 2019, it likely reached the sexually mature phase in winter 2020. The growth of the third cohort at the Jarun study site was faster at lower temperatures than at the Medsave study site, with higher average water temperatures in the HZ during the spring and summer (Table 6). The present research showed that the individual cohort life span in S. ambulans in the HZ of the Sava River could be more than 12 months long, which was more similar to the surface freshwater amphipod E. cari adapted to cold water streams, with a cohort life span of 6 to 12 months [136].

The influence of environmental factors (e.g., hydrological conditions, sediment type, temperature, food availability, water mass movements in the sea) on crustacean growth parameters has been reported by several authors [79,85,123,127,137,138,139,140]. A comparison of the estimated growth parameters of VBGF with other amphipod species revealed similar K values of S. ambulans with the polar species W. obesa from Antarctica, Gammarus wilkitzkii Birula, 1897 from the Arctic, and with N. minikoi from the Caribbean Sea (Table 4). This generally lower growth rate of S. ambulans compared to another amphipod, Pseudorchestoidea brasiliensis (Dana, 1853), may be related to the lower metabolic rate and oxygen consumption regularly found in hypogean amphipods [141,142]. Meagre oxygen consumption rates have also been measured in polar amphipods [143]. This supports the suggestion that a lower growth rate is related to lower oxygen consumption in species living in Arctic regions under sea ice [120].

Considering the population types defined by Schwamborn [144] (“type A” and “type B”) and using the “L_max approach” (maximum length approach, i.e., the largest organism in the sample) [145], we analysed the two populations in the Sava River. At the Medsave site, L_∞ (7.20) for females was greater than L_max (6.38), and L_∞ (6.30) for females was greater than L_max (6.00) at the Jarun site. This points to a “Type B” population characterised by slow growth and a very high Z/K ratio [144], but not a typical one, since L_max was not far below L_∞.

The total mortality rate of S. ambulans differed between males and females at both sites, with males exhibiting a higher mortality rate. High daily water level amplitudes are the first hypothesis that may explain the differences between mortality rates. Adult males (individuals with lower TBL) may remain in the upper layer of the sediment, while adult females (individuals with higher TBL) are more profound. In this way, males may remain without water and nutrients for a period of time, which could eventually lead to higher mortality. It was suggested that small individuals (juveniles and males) of the anomuran crab Emerita brasiliensis Schmitt, 1935 that remain near the substrate surface are more exposed to tide changes compared to larger females [87]. Consequently, this stratification by sex and size may lead to size-dependent mortality. Another hypothesis is that males are more vulnerable to predation because they need to find females to mate with [146]. The higher total mortality rates of females and males at the Jarun site compared to the Medsave site could be explained by lower food availability, as DOC and TOC concentrations were lower at the Jarun site (Table 6). The importance of food has been documented as one of the factors explaining the mortality rate of E. brasiliensis [87].

The natural mortality of amphipods can be caused by predation, disease, age, and environmental factors [93,147,148]. The previous study confirmed the high mortality rate of amphipods in the HZ under extreme hypoxic conditions [148]. This could be for the same reason as the high natural mortality of males (2.65 year⁻¹) in the HZ of the Medsave study site, due to strong water level fluctuations, and very often a low flow velocity. However, those conditions affected the Jarun population less, since males had lower natural mortality compared to males at Medsave. In addition, the Medsave and Jarun populations were under the extreme influence of water level fluctuations caused by hydroelectric power plants (HEPPs) in Slovenia, resulting in daily changing rewetting and drying fluctuations.

The growth–mortality continuum was introduced for benthic invertebrates, where ΔZ/k is a measure of relative predation pressure [149]. Calculated ΔZ/k values for S. ambulans were positive for males (0.56) and females (0.19) at the Jarun site and for females (0.09) at the Medsave site, suggesting higher than average predation pressure. It was suggested that S. ambulans may be preyed upon by the larger amphipods [7]. However, in the Sava River, only turbellarians were found in the samples from the HZ as potential predators of S. ambulans.

Differences in populations and species longevity in crustaceans could be influenced by variations in growth rate, mortality, and local environmental conditions such as temperature, concentration of organic matter, and oxygen [83,150,151,152]. The lifespan of amphipods can vary from a few months to over 10 years [56]. Species with long life spans usually live under harsh environmental conditions (e.g., low temperatures, high seasonality of food availability) [120,153]. S. ambulans showed high longevity, similar to the amphipods W. obesa (>5 years) [66] and G. wilkitzkii (5 years) [120]. S. ambulans had a higher longevity compared to other amphipod species that live in temperate waters, such as Echinogammarus marinus (Leach, 1816) and C. pseudogracilis, whose longevity can be as long as 1 and 2 years, respectively [83,129]. The higher longevity of females at the Medsave site compared to females at the Jarun site may be related to the lower total and natural mortality of females at Medsave. However, this pattern could not be explained in males, which had higher longevity and total mortality at the Jarun site but lower natural mortality than at Medsave. The influence of higher mortality rate on shorter life span was shown in E. brasiliensis [127].

4.3. Environmental Factors and Life History Traits of S. ambulans

Several environmental factors in the Sava River were recorded to have a substantial impact on population dynamic and general life history traits of the species S. ambulans. Temperature, as a key environmental factor, has a great influence on the growth and number of moults of crustaceans [60,154], with the growing rate related to high temperature values. Species from cold climates have a longer life span and larger body size, while species from warm climates have a shorter life span and reach smaller body sizes [56]. In addition, temperature has been shown to affect the macroinvertebrate community structure [155,156]. Synurella ambulans is a euryoecious organism (adapted to different habitats) [9] that can tolerate a wider range of physico-chemical parameters, including temperature, although it belongs to the cold stenothermal group of amphipods [13], with some races highly adapted to stable subterranean environments [20]. The higher maximum water temperature at the Medsave study site (23.3 °C) could lead to a faster growth rate of males at this site, compared to the Jarun population with a narrower temperature range and lower maximum water temperature values (16.5 °C) (Table 6). Optimal water temperature is important for crustacean longevity [83]. In GWCE, such as the HZ, the optimal water temperature is equal to the mean annual air temperature. Accordingly, the optimal physiological temperature of S. ambulans in the Sava River plain is 11.6 °C (for the period from 1862 to 2016) [157]. According to the Climate Atlas of SFRY (for 1931–1960) [158], the mean annual temperature was 10–11 °C and, according to Zaninović et al. [159], 11–12 °C. Environmental temperature may affect the duration of the reproductive period; for example, in the isopod Excirolana braziliensis Richardson, 1912, reproduction was seasonal and may have been interrupted due to temperature fluctuations in the sea [160]. In our study, reproduction was also seasonal, as ovigerous females occurred only in winter. This can be justified by the claim that S. ambulans, as a representative of the family Crangonyctidae, is in fact a cold-water species [13]. The number of recruitment events that occurred in spring and summer did not differ between Jarun and Medsave populations. The temporal scale (seasonal sampling) used in this study could affect the interpretation of growth results, as S. ambulans populations could change over shorter periods of time. Therefore, a more detailed (monthly) study with longer data series is needed to provide a more accurate estimate of cohort numbers and growth parameters. However, due to the decreasing number of gravel bars along the Sava River and the potential vulnerability of the species to extinction, we believe that such an intensive sampling campaign would not be acceptable. It was shown that populations of E. braziliensis are influenced by daily changes [161]. Therefore, short-term hydrological variations, such as abrupt daily changes in water level amplitudes of the Sava River, may also influence changes in the population structure.

Our results highlighted that S. ambulans are distributed unevenly in the HZ throughout the year based on the influence of various environmental factors on gender and ontogeny classes. Among crucial environmental factors in the HZ were food amount belonging to the community and abundance balancing variables where higher nitrate promoted the decomposition of leaf litter via stimulated detritivore feeding activities [162], which also includes S. ambulans. Determined relations in distribution patterns of different genders and groups of organisms according to ontogeny and environmental parameters in our analysis were confirmed through the oscillation of organic matter with the depth of the samples. Moreover, part of the organic matter may remain below the streambed during flooding or sediment movements at higher water discharge [162], which also occurs in the investigated study sites due to the large influence of the HEPP in the upstream parts of the Sava River. The investigated species S. ambulans, like most hyporheic invertebrate communities, responds to the seasonal flow recession, which was documented even for macroinvertebrates entering the HZ in a groundwater-dominated stream [163]. CCA analysis showed strong relations with temperature and oxygen as the most important environmental parameters in the HZ, influencing growth and development in both gender and ontogenetic stages in S. ambulans. Both environmental factors have been shown as regulation parameters for growth and development in crustaceans [60,83]. General suboptimal environmental conditions in the HZ of the Sava River during the summer (higher temperature values and lower oxygen, Table 6) due to lower hyporheic flow could be the reason for the negative relation with the abundance of lower ontogenetic classes (adolescent males and females) as more sensitive to warming and hypoxic conditions.

4.4. Advantages and Limitations of Length-Based Methods for Amphipod Studies

Size-frequency growth models are an important source of information on age and longevity in crustaceans because they are powerful, cheap, easy to apply, and require very few resources [139,164]. Although other, more accurate methods for age determination exist (e.g., tagging and recapture, lipofuscin-based ageing, radiometric ageing), these size-frequency-related growth models are still extremely relevant, especially for many small shrimps [165]. Size-frequency growth models have been shown to be more reliable for short-lived species (e.g., mysids such as the opossum shrimp N. americana (10 months)) than for long-lived species (e.g., decapods such as the crayfish Pacifastacus leniusculus (Dana, 1852) (16.7 years)) [91,139,166]. However, such growth analyses have limitations and must be used and interpreted with caution [164]. The need for additional data or assumptions is a common limitation of length-frequency based models, and alternative conclusions are possible depending on them [166]. In addition, it has been shown that methods used to fix or constrain L_∞ a priori can introduce bias in their L_∞ estimates [144]. Methods used for the unconstrained search for optimal combinations of growth parameters, such as Response Surface Analysis, do not provide 95% confidence intervals for estimates (e.g., K and L_∞) [165]. Current length-based methods were critically evaluated and a new pseudo-R² index for the goodness of fit of von Bertalanffy growth models (a low index implies a “poor” dataset and small sample size) was proposed [165]. This new bootstrap-based method offers a promising new direction for length-based analyses of growth and mortality in natural populations, especially due to its high reproducibility and accuracy.

5. Conclusions

The environmental adaptation of S. ambulans to the hyporheic interstitial zone in the present study could originate from eco-physiological adaptations of the species to better cope with the habitat intermittency (huge daily rise and fall of the water level in the Sava River coming from the upstream parts in Slovenia due to the influence of HEPPs), huge oscillations of the water current, and consequently fluctuating environmental parameters and steep abiotic and biotic (life history) gradients. Another factor influencing the adaptation capacity could be connected to the directional selection of life in subterranean environments where energy demand is lower relative to surface habitats.

Thanks to the present study, we noticed that S. ambulans is a very active burrower in the HZ and, as such, belongs to the key bioengineers unclogging hyporheic interstitial space because the spatial and functional suitability of the HZ is impaired when the interstitial pore space becomes clogged. Thanks to burrowing, the activity of the species plays a vital role in the vertical connection and exchange between free-flowing surface water, pore water of the HZ, and groundwater (upwelling and downwelling), and thus we can declare the species to be a guard of the HZ and a sentinel of climate change.

In our study, intrinsic (result of aging) and extrinsic (environmental caused) factors influenced the natural mortality rates in S. ambulans populations, increasing in the environmental conditions with expected higher metabolic rates (higher temperature, lower oxygen), and lowering the longevity of populations. The size of females reaching sexual maturity changes seasonally and mainly depends on water temperature. The shallow hyporheic interstitial zone may only represent a refuge for S. ambulans, which means the species does not complete its entire life cycle in this environment because it uses the hyporheic interstitial as a spawning habitat, which was confirmed by the presence of juveniles, but not as a nursery place due to the lack of ovigerous females. A fine-scale life history trait analysis covering at least three years of observations is needed for a complete overview of cohorts.

The species S. ambulans acts as an essential link in the transfer of matter and energy in the hyporheic interstitial zone, and knowledge of its population dynamic and life history traits is crucial for the conservation strategy of the species and its environment. Due to the ecological relevance of the species to the hyporheic ecosystem, it is important to protect the gravel bars it inhabits along the Sava River and European rivers in general. Our results may help future researchers to conserve the GWCE as a whole, especially regarding the conservation purposes of S. ambulans as a species inhabiting the vulnerable HZ of the Sava River.