Asellus aquaticus can be viewed as a species complex and its taxonomy is not fully resolved (Verovnik et al., 2009; Protas et al., 2023). However, it should be noted that it currently contains two described species (Verovnik et al., 2009) as well as numerous subspecies (Sket, 1994). Partially in agreement with the current taxonomy, the latest large-scale phylogenetic study found that there might be as many as thirteen operational taxonomic units in Europe (Sworobowicz et al., 2015). For simplicity, we treat all sampled populations included in this study as A. aquaticus. While taxonomic and phylogenetic status should not affect our predictions regarding within-population variance, it might do so on the among-population level. The implications of this are detailed in the discussion.

A total of 825 individuals were collected from eight cave and nine surface populations of A. aquaticus, spanning Italy, Slovenia, Hungary and Romania (Table 1; Figure 1). The dataset used in this study is identical to that analysed in Balázs et al. (2021), but it is supplemented with two new Romanian populations (Aleea Cetăţii Well, Hagieni Spring) and individuals where missing measures of traits were present. The latter had to be omitted from the previous study as the statistical methods did not allow their inclusion.

In Italy and Slovenia, the surface and cave populations live in cold flowing freshwater. They inhabit sinking rivers, where the surface populations inhabit the surface, while the cave populations inhabit the underground parts of the same rivers. The only exception is the surface population collected from the Ljubljana Moors, which lives in stagnant freshwater. In Hungary, two populations, the Molnár János Cave and the Malom Lake live in thermal water, where the lake is fed by thermal groundwater emerging from the cave. The other two surface populations, the Csömör Stream and Dunakeszi Lake inhabit flowing (former) and stagnant (latter) non-thermal freshwater. In Romania, the Aleea Cetăţii Well, the Dimitru Ana Well and the Kara Oban Lake are parts of the same thermal, sulphidic, ground-surface water system that flows through the chemoautotrophy based Movile Cave (Feurdean et al., 2003; Flot et al., 2014). The Aleea Cetăţii and Dimitru Ana wells represent the cave habitat, while Kara Oban Lake represents the surface habitat. The Hagieni Spring is an access point to another thermal, sulphidic, groundwater system which is hydrologically separated from the former system. The Băile Turceşti Spring represents a non-thermal, non-sulphidic water habitat. Further information on the sampled populations can be found in Konec et al. (2015) and Balázs et al. (2021).

For simplicity we use the term ‘population’ throughout the paper, but we note that the connectedness of the sampling locations in Romania is not fully resolved. Although the populations of the Aleea Cetăţii and the Dimitru Ana Wells are likely not fully separated by physical barriers, their mean trait values, standard deviations, and slopes of trait - body size regressions are markedly different (data not shown). Further, a previous analysis of microsatellite data from the same system revealed that samples from the different sampling localities form separate clusters (Konec et al., 2015). Hence, we treated them as different populations.

Sampling of the populations was carried out using hand-nets and by washing off the rocks. Diving was necessary in Molnár János Cave, where a modified Sket-bottle (Chevaldonné et al., 2008) was used. Sex determination was performed on each individual based on the gonopod morphology with a Zeiss Stemi 2000 stereomicroscope (Carl Zeiss AG, Oberkochen, Germany). The collected individuals were stored at 4 °C in RNALater (Thermo Fisher Scientific Inc., Waltham, Massachusetts, USA) in labelled vials until slide preparation (see below), except for the animals collected from the Aleea Cetăţii Well and the Hagieni Spring, which were stored in 96% ethanol. The difference in the storage methods is highly unlikely to affect estimated variances, especially after we accounted for body size (see below). All samples and slides are deposited at the Department of Systematic Zoology and Ecology at Eötvös Loránd University in Hungary. Only animals larger than 3.5 mm were kept after sorting, as this is the approximate size limit for adult individuals (Bloor, 2010). For sample sizes and sampling site details, see Table 1 and Figure 1. The uneven sample size among populations and sexes is the result of the exclusion of juveniles after sampling, skewed sex ratio during collection, differences in available animals during sampling and the destruction of certain samples during preparation. The sample collections were done under the permits: Slovenian Environmental Agency (ARSO), document nos. 35602-46/2016-4; 35602-41/2021-5 (Slovenia); Regione Autonoma Friuli Venezia Giulia (Italy), Direzione Centrale Infrastrutture e Territorio, PROT. TERINF-GEN-2017 n. 0019455/P; DINPI/5927-0/2020 (Hungary).

We measured 12 functional morphological traits which represent six functional categories (Figure 2). Additionally, we measured body length to account for body size. Body width is expected to positively correlate with female fecundity (Ridley and Thompson, 1979; Vick and Bloom, 2010). Antenna I, where the total length was measured, is responsible for chemical sensing (Heimann, 1984). Antenna II has a mechanosensory function and aids the detection of food and mates (Bertin and Cézilly, 2005). Although flagellum length would be the best approximation of function, this was broken in most of the sampled individuals. Detecting and quantifying the brake is difficult due to the high variability of the flagellar articles (Prevorčnik et al., 2004; Verovnik et al., 2009). Therefore, we measured this trait as the sum of the length of the V^th and the VI^th articles of the peduncle. All the traits measured on pereopod I (propodus length, propodus width, dactylus seta row length, dactylus length, and spine length) are involved in feeding and grooming (except spine length) and probably have a role in both interspecific and intraspecific aggressive interactions (Bertin et al., 2002; Bauer, 2013). The measured traits of pereopod IV (total length, propodus polygon length, and propodus arbitrary length) have different functions in females and males. In females, it serves as a walking leg; in males, it is primarily used to hold the female during precopulatory mate guarding (Ridley and Thompson, 1979). Pereopod VII, measured as total length, has a general walking function and might aid sensing, as elongated appendages increase the available surface for the sensory setae that are aligned along the pereopods (Fišer, 2019).

Body length and width were measured using TpsUtil. v. 1.74 (Rohlf, 2015) and TpsDig2 v. 2.30 (Rohlf, 2017) software on images taken both with a Canon 600D camera (Canon Inc., Tokyo, Japan) and with a Zeiss Stemi 2000C (Carl Zeiss AG, Oberkochen, Germany) microscope with standardised position in the centre, a microscope light and a 3 cm long scale of millimetre precision. Appendage measurements were conducted on microscope slides using a Zeiss Axioscope II microscope and the AnalySIS Program Package (Carl Zeiss AG, Oberkochen, Germany). The microscope slides were prepared by mounting the appendages on regular glass slides in Kaiser’s glycerine gelatine (Merck KGaA, Darmstadt, Germany) with the dorsal side upwards. This was then covered with a cover glass. All morphological measurements were carried out by the same person.

First, we built Bayesian double-hierarchical generalised linear models (DHGLMs) with a Gaussian distribution, separately for each trait and sex. Trait value was the dependent variable, habitat (cave vs. surface) and body length were the independent variables and population nested in habitat was the random effect. Body length was included to account for differences in body size. The model was built using the brms package (Bürkner, 2017, 2018) which is based on the Bayesian software Stan v. 2.26 (Carpenter et al., 2017; Stan Development Team, 2021). The DHGLM approach enables simultaneous estimation of a ‘mean model’ and a ‘residual model’. The mean model estimates mean differences among habitats, and the residual model estimates the extent of variance within each habitat. It is important to note that the intercept of the dispersion model (i.e., habitat mean residual standard deviation; sigma) is estimated on the log scale to ensure that standard deviations are always positive (Cleasby et al., 2015). Each model was run using weakly informative priors. The models were specified as in Equation 1.

Four chains were run, with 80,000 iterations, a warmup of 5,000 iterations and a thinning interval of 40. Although thinning of chains is often unadvised, the high number of iterations could lead to high autocorrelation (Link and Eaton, 2012), therefore it was applied. Model convergence was accessed for each model via inspection of traceplots, Ȓ convergence diagnostics (see Bürkner, 2017; Vehtari et al., 2020), and effective sample sizes (Vehtari et al., 2020). Trace plots indicated model convergence for all DHGLM models, Ȓ was always < 0.01 and the bulk and tail effective sample sizes were generally above 7,000, though always higher than 6,000.

To report standardised measures in cases of both significant and non-significant mean differences between the habitats, we calculated within-population (V_wp) and among-population (V_ip) variances for each habitat type for each sex as specified by Royauté and Dochtermann (2021). We assessed differences in the variances between the habitats by simultaneously inspecting whether the 95% credible interval of variance differences between the habitats exclude 0 and by inspecting the log variance ratios. Log variance ratios > 1 imply high biological relevance while log variance ratio < 1 imply low biological relevance (Royauté and Dochtermann, 2021).

Subsequently, we conducted an analysis on the variances extracted from the DHGLMs using an approach based on the general workflow of meta-analysis, but with data collected from a single study. Here we applied Bayesian generalised linear mixed models (GLMMs) with Markov chain Monte Carlo methods using the MCMCglmm package (Hadfield, 2010). Although the number of individuals sampled per population was uneven, the variance caused by sampling error should be minimal, as the sampling method was the same in each case. The variance caused by measurement error should also be negligible, as all measurements were carried out by the same person. Therefore, we only accounted for the uneven sample sizes in the variance in effect sizes. It was calculated following Holtmann et al. (2017) (Equation 2).

To account for the non-independence of the variance components (Nakagawa and Santos, 2012), we applied multilevel models to estimate effect sizes for all traits, estimating their global effect size. The model was fitted using an uninformative inverse gamma prior following Wilson et al. (2010) with 3,000,000 iterations, a thinning interval of 500, and 1,000,000 burn-in. First, we ran an intercept-only mixed model, with ‘morphological trait’ as a random effect, to determine the mean variance. Next, we constructed meta-regression models to test for the effects of moderators: ‘sex’ (female vs. male) and ‘habitat’ (cave vs. surface) (Nakagawa and Santos, 2012). We conducted fixed-effect mixed models to estimate the mean effect size for each moderator separately (Winternitz et al., 2017), and then we estimated the interaction term in a separate model. Variance differences between the different levels of the moderators were accessed by inspecting whether the 95% credible interval of variance differences excludes 0. To assess which model fits our data best, we compared deviance information criterions (DIC), where smaller DIC indicates a better fit. However, minor differences of up to three indicate similar goodness of model fit (Spiegelhalter et al., 2002). In the within-population variance model, all models performed similarly (DIC: intercept model -186.98; habitat model -186.59; sex model -186.21; interaction model -185.11). In the case of the among-population model, the habitat × sex interaction model performed the best (DIC: intercept model 4.07; habitat model -23.09; sex model 0.58; interaction model -39.66). To estimate the heterogeneity of the effect sizes and to understand the proportion of variance explained by our moderator variables, we partitioned total heterogeneity (I^2total) and residual variance (I^2residual) of all original fixed effect models and the interaction model using I² statistics (Higgins and Thompson, 2002). We evaluated heterogeneity estimates considering those around or below 25% to be low, those around 50% to be medium and those around or above 75% to be high (Higgins et al., 2003).

We report the results for each sex and trait separately here, which results in a combination of a total of 24 cases. Within-population variance differed significantly between the habitats in 12 out of 24 cases; three trait variances differed only in females; one trait differed only in males and four traits differed in both sexes (Table 2). Within-population variance in these cases was always larger in caves. The absolute log variance ratio was > 1 in three traits (all in the case of females) and it was < 1 in nine traits (Table 2). Therefore, the biological significance of the latter nine cases is limited (Royauté and Dochtermann, 2021).

Among-population variance differed between the habitats in eight out of 24 cases; three traits differed only in females; three traits differed only in males and one trait differed in both sexes (Table 2). In all these cases morphological variance among-populations was higher in the cave than in the surface habitat. For trait level patterns, see Table 2. The absolute log variance ratio was always > 1 (Table 2), which indicates that these patterns are biologically relevant (Royauté and Dochtermann, 2021). For the original model estimates of the variances see Supplementary Figure 1 and Supplementary Table 1.

The V_wp did not differ between the sexes (mean estimate for variance differences [95% Credible Interval of the mean estimate]: -0.007 [-0.05 – 0.036]) neither between the habitats (0.018 [-0.025 – 0.06]). There was also no difference between the different levels of the interaction term of habitat × sex. This means that variances of neither females nor males differed between the two habitats (cave female vs. surface female: -0.024 [-0.036 – 0.082]); cave male vs. surface male: 0.011 [-0.05 – 0.07]) and there was also no difference between the sexes within the cave or within the surface habitat (cave female vs. cave male: -0.001 [-0.062 − 0.06]; surface female vs. surface male: -0.14 [-0.077 – 0.46]) (Figure 3). Total heterogeneity among traits in V_wp was low (2.68%) and remained close to constant at approximately 1% in all models.

The V_ip differed between the sexes (-0.155 [-0.29 − -0.017]) with males being more variable than females. It also differed between the habitats (0.314 [0.214 − 0.425]) with caves being more variable than the surface. There were also differences between the different levels of the interaction term of habitat × sex. The interaction model showed that V_ip differs between the habitats in both males and females (cave female vs. surface female: 0.189 [0.058 – 0.313]; cave male vs. surface male: 0.434 [0.318 – 0.562]), being larger in the cave environment (Figure 4). The model also revealed differences in the variance between the sexes within the cave habitat (cave female vs. cave male: -0.276 [-0.4 − -0.147]) with males being more variable, though no such difference was observed in the surface habitat (surface female vs. surface male: -0.031 [-0.152 – 0.091]). Total heterogeneity in V_ip among traits was medium (45.97%). Residual variance remained moderate in the sex model (43.31%), but after including habitat (28.72%) or habitat × sex (20.23%) as moderators, residual variance decreased to relatively low and low, respectively.

Our original prediction was based on the assumption that environmental fluctuation in a population’s habitat is expected to positively correlate with both genetic and morphological variation (Sultan and Stearns, 2005; Bell, 2010). This is because fluctuating environments (such as surface habitats) are expected to result in changes in the strength and direction of natural selection (Dingemanse et al., 2004; Siepielski et al., 2009), which should lead to an increased variance. Conversely, the relative stability of caves should lead to decreased variance. Kralj-Fišer et al. (2025) found decreased variance in cave A. aquaticus populations when comparing the variation in movement activity between cave and surface ecomorphs, supporting the fluctuating selection hypothesis. While the overall morphological variance was similar between the habitats, in contrast with our predictions, we found higher variance in cave populations in certain traits, namely the traits connected to pereopod IV. Quantification of morphological variance is generally scarce (Chevin et al., 2015), therefore, there is no predetermined limit to what constitutes small and large variance. Differing overall variances would have been straightforward to interpret, but in the case of similar variances, defining morphological variance on an absolute scale (small, medium or large) as well as on a relative scale (by comparing the habitats) would be crucial to disentangle the observed pattern.

In addition to the decrease in morphological variance a decrease in genetic diversity would also be expected in caves (Poulson and White, 1969). This is not only due to the absence of fluctuating selection but also because of the small number of colonisers from the surface (founder effect), small population size, and high degree of isolation, which results in limited or absent gene flow with other populations (Garant et al., 2007; Trontelj, 2018). Such a phenomenon has been observed in cave populations when compared to surface populations in several taxa (Carlini et al., 2009). However, genetic diversity is not always low in cave populations. Substantial genetic diversity was found in some cave animals at times with even higher genetic diversity than surface conspecifics or closely related taxa (Stepien et al., 2001; Lejeusne and Chevaldonné, 2006). Although this was not tested here, genetic results regarding A. aquaticus are mixed, but the here studied cave populations generally show no decrease in diversity. Verovnik et al. (2003) found that most cave populations are similarly diverse to surface populations, such as the Trebiciano Cave and Planina Cave populations included in this study. Konec et al. (2015) also showed that in Romanian cave populations genetic diversity is similar to that on the surface. Additionally, they found that it is only slightly reduced in the Slovenian Planina Cave (the latter result is incongruent with Verovnik et al. (2003)). Higher genetic variation than expected is probably attributable to relaxed selection (Cui et al., 2019). Such processes have mostly been studied in relation to eye and pigment loss (Protas et al., 2011). When comparing the gene expression patterns of optimal eye function in blind cave and sighted surface fish, Stern and Crandall (2018) found that relaxed selection was responsible for the differences in gene expression. The role of this mechanism was previously proposed for two morphological traits of A. aquaticus by Konec et al. (2015).

Taken together we propose that differing selection intensity/efficiency on the different traits paired with similar genetic variance in caves as on the surface and no gene-flow between caves leads to similar variance as inherent of fluctuating selection in the connected surface populations. At the same time understanding the reported V_wp patterns warrants further targeted studies.

Regarding among-population variance, this study was exploratory, as two contrasting mechanisms with contrasting outcomes could be at play in caves. The similar environmental characteristics across caves should result in lower variance, whereas the long time since isolation coupled with mutation and drift, should lead to higher variance. Along with the predictions of the latter, we found that V_ip is higher in caves in both sexes. This paired with the previously found mean differences between populations within the habitats (Balázs et al., 2021) clearly implies that morphological diversification within the cave habitat is more pronounced than on the surface. Similar pattern was found by Zhao et al. (2024) who detected larger interspecific variation in body shape and appendage length in cave-dwelling than in surface-dwelling diving beetles (Dytiscidae). They attributed these differences to genetic drift and neutral evolution.

In caves - where the population size is expected to be lower due to sparse resources (Poulson and White, 1969) - drift or inbreeding, coupled with a lack of gene flow between the populations can also result in large differences in the gene pool and genetic variation (Garant et al., 2007; Willi et al., 2007). However, the observed pattern also suggests that diversifying selection based on small environmental differences can be highly effective when the environment is stable and the populations are genetically isolated. On the contrary, dispersal is expected to decrease genetic and morphological variation among surface populations (Wright, 1978).

Additionally, we must note that, while in this study we did not measure environmental parameters, at least in terms of water temperature and chemistry (thermal vs. cold; sulphidic vs. non-sulphidic), the studied cave habitats are not uniform. This could result in divergence among populations. In the case of the fish, Poecilia mexicana, which occupy cave and surface habitats with both sulphidic and non-sulphidic water, Tobler et al. (2008) found that gene expression and morphology differed systematically between both the habitats and the water types. Given the low number of known cave-adapted A. aquaticus populations resulting in few replicates, such a comparison is not possible in our system, thus generalisation of our results require caution. At the same time environmental differences are also present among the surface habitats, as some of the surface populations inhabit the water emerging from the studied caves. This suggests that explanations for the differences solely based on water temperature and chemistry are highly unlikely, but it is possible that other unmeasured parameters play a certain role in the observed pattern. Taken together, we suggest that the diversifying selection scenario is more likely in our system than strong relaxed selection and neutral evolution because of (i) the marked morphological adaptations that have been reported following cave colonisation in our system (Balázs et al., 2021), which is against the assumption that solely relaxed selection shapes the traits and (ii) the marked sexual differences in V_ip within the cave habitat (see below).

We should also mention that, at first glance, the presence of a single described species in the surface and two described species in caves could seem like a straightforward explanation for the results obtained. The latest detailed phylogenetic study of the species complex (Sworobowicz et al., 2015) implies that operational taxonomic unit diversity is higher in surface than in cave populations. Therefore, we argue that the observed pattern is not inherent of phylogenetic diversity.

Trait difference between the sexes can stem not only from differences in the means, but also from differences in variance. Sex-biased differences in phenotypic variance have lately resulted in scientific debate (Harrison et al., 2022; Del Giudice and Gangestad, 2023), with empirical studies yielding controversial results. Certain studies reject the greater male variability hypothesis (Harrison et al., 2022), others only detect it in certain traits. A detailed study on mice revealed no general variance difference between the sexes when looking at a large set of phenotypic traits, although the hypothesis was supported when solely analysing morphological traits (Zajitschek et al., 2020). In our study, we could not detect overall variance differences between the sexes in V_wp in either habitat. This means that the greater male variability hypothesis is not supported in our system. It should also be added that although discussing the results trait-by-trait would be too speculative, here an interesting pattern arose that cannot be disregarded. While no overall difference was detected in V_wp between the sexes, biologically relevant habitat differences, with cave adapted populations having higher variance, were found in females in traits connected to pereopod IV. This trait has a walking function in females, but in males this trait has a prime function in mate acquisition (Ridley and Thompson, 1979). This pattern suggests that while the general elongation of this trait happens after colonising caves (Balázs et al., 2021) and the trait becomes more variable in females, both elongation and variation is restricted in males as it could distort the function.

The V_ip was the highest in cave males, but cave females also showed substantial variance when compared to surface conspecifics. This suggests a key role of sexual selection (besides survival and fecundity selection) in adaptation to the cave environment and can be seen as an additional argument against neutral evolution being a dominant source of among-cave morphological variation in A. aquaticus. Studies on subterranean species generally found that the sex ratio is female-biased, but with substantial differences across populations or closely related species (Culver and Holsinger, 1969; Vonk and Nijman, 2006; Premate et al., 2021). Although no study has yet explored the sex-ratio in cave-adapted A. aquaticus, sex-ratio variation could result in differences in the strength of sexual selection in males between the different populations and explain their higher variance.

As an addition to the results, it should be noted that this study does not allow us to clearly distinguish between the genetic and the environmental components of the observed variance. However, presumably the observed patterns reflect their interplay (Vogt, 2020). The sampled cave populations are genetically isolated from surface populations and from each other (Protas et al., 2023). The split between the cave and surface populations dates back at least 60,000 years in Hungary (Pérez-Moreno et al., 2017) and at least 200,000 years in Slovenia, Romania and Italy (Verovnik et al., 2009). Such long isolation coupled with the marked environmental differences between cave and surface habitats has led to a considerable morphological divergence (Prevorčnik et al., 2004; Konec et al., 2015; Balázs et al., 2021). While the genetic basis of eye and pigment loss is relatively well studied in A. aquaticus (Lomheim et al., 2023; Fišer et al., 2024), there was only a single study on another morphological trait namely the antenna II (Mojaddidi et al., 2018). Their results imply that the differences have a genetic basis. Taken together, the long genetic isolation and the proven complex adaptations make cave-surface morphological differences highly likely to be the result of evolutionary (genetic) mechanisms rather than phenotypic plasticity. At the same time, as caves are characterised by relatively stable environments (Lauritzen, 2018), even under the unlikely assumption that all observed patterns are solely due to phenotypic plasticity the predictions of this study would have remained the same.