Sample collection

A total of 195 individuals of C. gibelio comprising 52 diploid sexual males (2nM), 52 diploid sexual females (2nF) and 91 triploid gynogenetic females (3nF) (total length and body weight, mean ± SD: 28.80 ± 2.39 cm and 399.77 ± 97.24 g for triploid females, 28.21 ± 2.91 cm and 374.29 ± 103.47 g for diploid females, 27.80 ± 2.66 cm and 328.00 ± 83.56 g for diploid males) were sampled by electrofishing from the natural habitat. Fish were sampled in the same season of the year (the first half of August, corresponding to a hot summer in central Europe) across four consecutive years. All individuals originated from the River Dyje near Břeclav (48.8027356 N, 16.8385806 E; the River Morava basin, Czech Republic). A blood sample was collected from each individual by puncture of the caudal blood vessel using a heparinized syringe and stored in 1.5 ml tubes. Each fresh blood sample was used for haematological profile analysis and oxidative burst analysis. For the following immunological analyses (i.e. complement activity and IgM level), a portion of each blood sample was centrifuged, and the plasma was collected and preserved in a freezer. A fin clip (1 cm²) was collected from each individual and placed in 75% ethanol for ploidy determination and DNA genotyping (Pakosta et al. 2018).

After blood and fin collection, fish were immediately placed in a tank containing the original water and transported to the laboratory under aeration. Individual fish were sacrificed humanely by a sharp blow to the cranium and severing the cervical spine. The sex, standard length (in cm), and body weight (in g) were determined for each individual. In addition, the hepatopancreas, spleen, and gonads were weighed. The following indices were calculated (see Šimková et al. 2015): Fulton's condition factor (FK), reflecting the overall somatic condition and potential difference in growth performance; the spleen-somatic index (SSI), reflecting the immunocompetence; the gonado-somatic index (GSI), reflecting the investment into reproduction; and the hepato-somatic index (HSI), reflecting fish vitality and energy resources.

The sampling, transportation, and maintenance of experimental fish, as well as the killing method, complied with legal requirements in the Czech Republic (§ 7 law No. 114/1992 about the protection of nature and landscape and § 6, 7, 9 and 10 regulation No. 419/2012 about the care, breeding and using experimental animals). The experiment was approved by the Animal Care and Use Committee at the Faculty of Science, Masaryk University in Brno (Czech Republic).

Parasites

Experimental fish were examined for the presence of parasites using the methodology described by Ergens & Lom (1970) within 36 h after transporting fish to the laboratory. Fin and gills were checked for the presence of metazoan ectoparasites, and all internal organs were examined for the presence of metazoan endoparasites using Olympus SZX7 stereomicroscope. In addition, the numbers of Ichthyophthirius multifiliis on fish gills were counted. Using fine needles, the individual parasites (except Myxozoa and I. multifiliis) were removed from external and internal organs immediately after fish killing and dissecting and were fixed.

Different preservation methods were used for the fixation of different parasite groups. Monogeneans were removed from gill arches and fins and fixed using a mixture of glycerine and ammonium picrate (GAP) (Malmberg 1957). Representatives of Digenea, Cestoda, Hirudinea and Crustacea were fixed using 4% formaldehyde (Ergens & Lom 1970). Finally, representatives of Nematoda were fixed by a mixture of glycerine and 70% ethanol (Moravec 1994). Individual parasites were identified using an Olympus BX50 light microscope (with an Olympus DP71digital microscopic camera) equipped with phase contrast and DIC (differential interference contrast, according to Nomarski) and digital imaging software (MicroImage 4.0 for Windows). For details on parasite data, see Pakosta et al. (2018).

Haematological analyses and glucose

Erythrocyte count, leukocyte count, and haematocrit value were assessed according to Svobodová et al. (1991). Immediately after blood sampling, the erythrocyte count (in T.l^–1) and leukocyte count (in G.l^–1) were calculated using a Bürker haemocytometer after staining with Natt-Herrick solution at a ratio of 1:200 (Svobodová et al. 1991). Haematocrit was measured using heparinized microcapillaries (75 mm long). Samples were centrifuged in microcapillaries using a haematocrit centrifuge (12,000 g for 3 min).

The glucose level in plasma was analysed using a commercial enzyme kit following the manufacturer's instructions (Glu L 1000, PLIVA-Lachema, Czech Republic). Samples in duplicates were analysed using a plate reader (Tecan Sunrise, Switzerland), and the concentration of glucose (mmol.l^–1) was calculated according to the standard concentration of glucose.

Immunological analyses

Oxidative burst refers to the process in which phagocytes of the immune system produce reactive oxygen metabolites, which they use to attack pathogens that overcome the natural protective barriers of an organism. Following Buchtíková et al. (2011), the reaction mixture included 50 times-diluted blood in Hank's balanced salt solution, luminol (Molecular Probes, Eugene, Oregon, USA; dissolved in borate buffer; pH = 9; final concentration, 10^–3 mol.l^–1) and Zymosan A as an activator (Zymosan A from Saccharomyces cerevisiae; Sigma-Aldrich, USA; final concentration, 0.25 mg.ml^–1). Oxidative burst was measured as the luminol-enhanced chemiluminescence recorded during one hour at room temperature (20 °C) using a LM01-T luminometer (Immunotech, Czech Republic). The integral of oxidative burst (measured in relative light units, i.e. RLU) is represented by the area under the kinetics curve (RLU*s).

Complement activity measures a nonspecific humoral immune response and plays a significant role in either classic (antibody-dependent) or alternative defence against pathogens. Due to its complexity, it participates in lytic, pro-inflammatory, chemotactic, and opsonization processes and connects nonspecific humoral and cellular mechanisms induced by phagocytes (Ellis 1999). The total complement (bacteriolytic) activity (i.e. comprising all activated pathways) was measured luminometrically following Kilpi et al. (2009), Nikoskelainen et al. (2002) and Virta et al. (1997). The recombinant strain of Escherichia coli K12 containing the luxABCDE gene, which originally comes from the terrestrial bioluminescent bacteria Photorhabdus luminescens (Atosuo et al. 2013), was exposed to fish plasma. The emitted lightning radiation was detected by an LM01-T luminometer (Immunotech, Czech Republic) and positively correlated with the viability of transformed E. coli. Complement activity is expressed as the period in which the viability of E. coli decreases in 50% of the total amount of bacteria; a shorter time represents higher complement activity in a fish plasma sample of the same concentration. Thus, complement activity is finally expressed in terms of inverted values (in h^–1).

The lysozyme concentration in skin mucus was determined in vitro using the radial immunodiffusion method in agar gel mixed with the chemoorganotrophic bacteria Micrococcus luteus (CCM 169). Mucus samples and lysome standards (Sigma-Aldrich), each of a total volume of 15 µl, were incubated in a wet chamber under laboratory temperature (20 °C) for 24 h. Afterwards, following Poisot et al. (2009), diffusion zone averages were calculated to determine the concentration of lysozyme in skin mucus (the resulting value is given in mg.ml^–1). As immunoglobulins represent essential components of specific immunity, IgM antibodies are primarily measured in the case of poikilothermic organisms. Total IgM in plasma was determined using specific precipitation with zinc sulphate (0.7 mM ZnSO4·7H2O, pH 5.8) (McEwan et al. 1970), following Šimková et al. (2015). Ig originating from the solution was removed by centrifugation, and IgM quantification (in g/L) was determined based on the total level of proteins in the sample using a commercially available kit (DC Protein Assay, Bio-Rad, USA) and a plate reader (Tecan, Sunrise, USA), both before and after precipitation.

Ploidy determination and DNA genotyping

A fin clip from each individual was used to determine the ploidy status of the fish. Before analysis, this tissue was homogenised on a Petri dish in a 2 ml solution of CyStain DNA 1 step PARTEC, and the relative DNA content was estimated using a Partec CCA I continuous flow cytometer (Partec GmbH; www.sysmex-partec.com). The same procedure was applied previously for ploidy measurements using blood samples (Flajšhans 1997, Halačka & Lusková 2000).

The control region of mitochondrial DNA (D-loop) was amplified for each specimen using the primers published by Papoušek et al. (2008) and protocols by Pakosta et al. (2018). Dneasy^TM Blood and Tissue Kit (Qiagen) was used for DNA extraction. The PCR products were purified using a High Pure PCR product purification kit (Roche) and directly sequenced using a BigDye Terminator Cycle sequencing kit on an ABI 3130 Genetic Analyser (Applied Biosystems, USA). The sequences were aligned using MEGA X (Tamura et al. 2013) and analysed using TCS 1.21 (Clement et al. 2000).

Statistical analyses

The measured variables were divided into two groups: immune (leukocyte count, complement activity, lysozyme, oxidative burst amplitude, SSI, IgM) and physiological (FK, GSI, HSI, glucose, haematocrit, erythrocyte count). Following Rohlenová et al. (2011), who showed PCA as a useful method when interpreting the overall physiological and immunological status of fish hosts, PCA was conducted on each of the variable groups to obtain indicators of the physiological and immune status of fish (hence decreasing the number of predictors). These indicators were represented by row coordinates of axes with eigenvalues > 1 (first two axes in both cases).

Non-metric multidimensional scaling (NMDS) was used to describe parasite assemblage composition (PCA was not used in this case, being prone to a problem of double-zeros). NMDS was performed using Bray Curtis distance on relative abundance (%) data, as the aim was to obtain ordination reflecting different types of parasite infection. The coordinates of the first two NMDS axes were considered indicators of parasite assemblage composition. Only parasite species with a prevalence greater than 10% were considered (see Pakosta et al. 2018 for prevalence).

Four main parasite load parameters were considered response variables: parasite abundance (log-transformed), parasite species richness, and the first two NMDS axes (representing proxies for parasite assemblage composition). We searched for the effects of four predictors for each response variable: two immunity indicators (the first two axes of PCA based on immune variables) and two physiology indicators (the first two axes of PCA based on physiology variables).

Mixed linear models (LMM) were constructed to reveal inter-group differences in all immunological and physiological variables represented by the PCA scores and also in individual physiological and immune variables, parasite abundance (log-transformed), parasite species richness, and parasite assemblage composition represented by NMDS scores with year as a random effect and with total length as covariate (in case of abundance and parasite species richness). Tukey HSD approach was conducted to control for type II error in multiple post-hoc pairwise comparisons (using functions glht and mcp from multcomp package, Hothorn et al. 2008).

Two different analytical approaches were conducted to determine whether parasite load was related to physiological and immunological predictors in fish groups. In the first analysis, we investigated whether there are inter-group differences in the effects of the predictors. We constructed a linear mixed model for each response variable that included fish group, four physiological and immunological variables (PH1, PH2, IM1, IM2) and their interactions with fish group as predictors and year as a random factor. We used a backward stepwise model selection procedure to reach the final model for each response variable. Term removal was based on Akaike information criterion (AIC) differences between the two nested models being < 2 and confirmed by conducting likelihood ratio tests between the two models. The presence of an interaction term in a final model demonstrated that the response of the parasite load parameter to the physiological and immune parameters present in the interaction term differed significantly between the fish groups. Fish total length (TL) was used as a covariate in parasite abundance and richness models. Furthermore, we have tested the effects of the same set of predictors on the overall assemblage structure in a multidimensional space (distance matrix based on Bray-Curtis distances), using multivariate permutational analysis of variance (PERMANOVA, Anderson 2001), following the same procedure (except for conducting likelihood ratio tests).

Fig. 1.

PCA of physiological (upper panel) and immune (lower panel) variables measured in three fish groups of C. gibelio, showing projections of physiological and immune variables into the first two PCA axes (a, d) and fish group coordinates along the first (b, e) and second (c, f) principal components (using scaling type 4, showing case scores and factor loadings). Percent of explained variability for each axis is shown in axis labels. Horizontal bar – median, box – interquartile range, whiskers – non-outlier range (1.5*interquartile range), points – outliers. FK – Fulton's condition factor, GSI – gonado-somatic index, HSI – hepatosomatic index, SSI – spleen somatic index, complement – complement activity, RLUa – oxidative burst amplitude, 2nM – diploid sexual males, 2nF – diploid sexual females, 3nF – triploid gynogenetic females.

Fig. 2.

NMDS of parasite assemblage composition: a) species scores projected onto the first two NMDS axes – only species with an occurrence > 10% are shown; b) and c) fish group scores along the first two axes: horizontal bar – median, box – interquartile range, whiskers – non-outlier range (1.5*interquartile range), points – outliers, 2nM – diploid sexual males, 2nF – diploid sexual females, 3nF – triploid gynogenetic females.

Table 1.

Factor loadings of the two principle components representing variables of fish host physiology (PH1, PH2) and immunity (IM1, IM2). Significant correlations with the PCA axes are in bold.

In the second analysis, we investigated which predictors best explained the variability of a response variable in each fish group. For each combination of fish group, response variable, and predictor, we constructed an LMM with sampling year as a random effect (for parasite abundance and parasite species richness, total length (TL) was used as a covariate). Thus, twelve tests (four for each fish group) were conducted for each response variable. For each model, we calculated AIC and marginal R², and we tested the significance of the predictors. We used the Benjamini-Hochberg procedure to control the false discovery rate in multiple tests (Benjamini & Hochberg 1995), with the procedure applied to all tests conducted with the same response variable and FDR set to 15%.

Allstatisticalanalyseswereperformedusingsoftware R, version 4.1.2 (R Core Team 2021), using libraries ade4 (Thioulouse et al. 2018), multcomp (Hothorn et al. 2008), MuMIn (Bartoń 2020), nlme (Pinheiro et al. 2021), lmtest (Zeileis & Hothorn 2002) and vegan (Oksanen et al. 2020).

Physiology and immunity in asexual and sexual gibel carp

The first two axes of physiology PCA explained 60.9% of overall variability. While the first axis (PH1) was correlated mainly with FK, HSI, and GSI (all positively) and erythrocyte count (negatively), the second axis (PH2) was positively correlated with glucose and haematocrit (Fig. 1a, Table 1). There were significant differences among fish groups along PH1 and PH2 (LMM, both df = 2,189 and P < 0.001). The three fish groups were differentiated along PH1 (post-hoc tests, all P < 0.001), with diploid males situated in the negative (left) part of the gradient, triploid females in the positive (right) part and diploid females in the middle (Fig. 1b). The differentiation was less clear along PH2, with diploid males reaching significantly higher values than both female groups (post-hoc test, P = 0.002 and P < 0.001), which were not distinguishable from each other (post-hoc test, P = 0.238, Fig. 1c). Individual data for physiological variables are shown in Fig. S1. Statistically significant differences were shown between diploid males and each of diploid and triploid females for haematocrit, GSI and FK (LMM, post-hoc test, all P < 0.001), whilst differences between all groups compared were found only for erythrocytes and HSI (LMM, post-hoc tests, all P < 0.005).

Fig. 3.

Relationship between the first two axes of immunity (IM1, IM2) and physiology PCA (PH1, PH2) and parasite abundance in the three fish groups (sexual diploid females – 2nF, sexual diploid males – 2nM, and gynogenetic triploid females – 3nF). Significant relationships are marked by predicted lines. The dashed line depicts trends for predictors that became non-significant after correction for multiple testing.

Fig. 4.

Relationship between the first two axes of immunity (IM1, IM2) and physiology PCA (PH1, PH2) and parasite species richness in the three fish groups (sexual diploid females – 2nF, sexual diploid males – 2nM, and gynogenetic triploid females – 3nF). Significant relationships are marked by predicted lines. The dashed line depicts trends for predictors that became non-significant after correction for multiple testing.

Table 2.

Parameters of final models predicting parasite abundance a) and parasite assemblage composition (b-d) of three fish groups of gibel carp. Parasite assemblage composition is predicted via linear mixed models using coordinates of first and second NMDS axes (b,c) as a response variable or using permutational multivariate analysis of variance (PERMANOVA) to determine predictor effects on parasite assemblage composition in multidimensional space d). Linear mixed models were used to predict response variables (a-c). For null models see Table S1.

The first two axes of immunological PCA explained 46.7% of overall variability. The first axis (IM1) was positively correlated with leukocyte count and complement activity and negatively with oxidative burst and lysozyme; the second axis (IM2) was negatively correlated mostly with SSI, IgM, and lysozyme and positively with oxidative burst (Fig. 1d, Table 1). There were no differences between fish groups along IM1 or IM2 (LMM, df = 2,189, P = 0.301 and df = 2,189, P = 0.071, Fig. 1e, f). Individual data for immune variables are shown in Fig. S2. The statistically significant differences between diploid males and each of diploid females and triploid females were found for lysozyme and SSI (LMM, post-hoc tests, lysozyme: P = 0.047 and P < 0.001, SSI: P = 0.009 and 0.007 for comparisons with diploid and triploid females, respectively), whilst the difference in IgM between triploid females and each of diploid males and diploid females were found (LMM, post-hoc tests, P = 0.019 and 0.045).

Parasite load in asexual and sexual gibel carp

While there was no significant difference in parasite species richness between the fish groups (LMM, df = 2,188, P = 0.104, Fig. S3), there were significant differences in parasite abundance (LMM, df = 2,188, P = 0.011, Fig. S3), with triploid females having significantly fewer parasites than diploid females (post-hoc test, P = 0.015) and marginally non-significantly fewer parasites than diploid males (post-hoc test, P = 0.083).

The first NMDS axis (NMDS1) was primarily associated with Gyrodactylus spp. and S. petruschewskii (positive values), whilst the second NMDS axis (NMDS2) was mainly positively associated with I. multifiliis and negatively with S. petruschewskii. There were significant differences among fish groups along NMDS1 (LMM, df = 2,189, P = 0.038) but not along NMDS2 (LMM, df = 2,189, P = 0.208). The coordinates of diploid females were shifted significantly more to the right (positive) side along NMDS1 than triploid females (Fig. 2, post-hoc test, P = 0.031).

Table 3.

Relationship between coordinates of the first two axes of immunology and physiology PCA (IM1, IM2, PH1, PH2) and four parasite assemblage characteristics: parasite abundance (log-transformed), parasite species richness, and coordinates of the first two NMDS axes determining the percentage composition of the parasite assemblage. For each combination of response and predictor, a P-value determining the significance of the predictor in a linear mixed model is presented. Regression coefficients are presented only for those predictors that were significant in the model. P-values and coefficients that were significant are in bold, those that became non-significant after correction for multiple testing are in italics. All significant terms represent the best models according to AIC and marginal R² (Table S2).

Effects of physiology and immunity on parasite load

The final model for predicting parasite abundance included fish total length, IM1, fish group and fish group:IM1 interaction (LMM, Table 2, Table S1), revealing different effects of IM1 on parasite abundance. The same IM1 effect on NMDS1 was observed with a final model, including IM1, fish group and fish group:IM1 interaction (LMM, Table 2, Table S1). None of the predictors significantly affected parasite richness (LMM, Table S1). NMDS2 coordinates were affected by IM1 only (LMM, Table 2, Table S1). Overall differences in parasite assemblage structure were affected by all four predictors (IM1, IM2, PH1, PH2), fish group and fish group:IM1 interaction (PERMANOVA, Table 2, Table S1).

Parasite abundance in both diploid males and diploid females increased with increasing IM1 (Table 3, Table S2, Fig. 3). Parasite abundance in triploid females was related to IM2 instead (Table 3, Table S2, Fig. 3). The parasite richness of both diploid males and diploid females were unrelated to any of the four immune and physiological predictors (Table 3, Table S2, Fig. 4). The parasite richness increased with decreasing PH2 in triploid females, this relationship was not significant after correction for multiple testing (Table 3, Table S2, Fig. 4).

Parasite assemblage structure, as expressed by NMDS1 was positively related to IM1 in diploid females and negatively to IM2 in diploid males, with no significant relationships in triploid females (Table 3, Table S2, Fig. 5). The coordinates of NMDS2 were negatively associated with PH1 and IM1 and positively with PH2 in diploid males, with the relationship with IM1 and PH2 becoming non-significant after correction for multiple testing (Table 3, Fig. 6). NMDS2 coordinates were also negatively associated with PH2 in diploid females, whilst no significant relationship was found for triploid females (Table 3, Table S2, Fig. 6).

Physiological components in asexual and sexual gibel carp

Considering the physiological components of fish vigour, we showed the clear differentiation between three fish groups, with sexual females positioned among sexual males and gynogenetic females along the first physiological axis (Fig. 1 and Fig. S1). Specifically, the differentiation along this axis was explained mostly by erythrocyte count and hepatosomatic index, reflecting the energy resources of organism, however, was also well explained by gonado-somatic index, reflecting reproductive potential, and Fulton's condition index reflecting body condition, growth performance, feeding activity, and nutrient availability (Fig. S1). All these variables except for erythrocyte count tended to reach their lowest values in sexual males when compared to both females, and in addition, higher values in hepato-somatic index in gynogenetic females when compared to sexual females. These findings may suggest some advantages for gynogenetic females when competing with sexuals in the same habitats and/or even when competing with sexual females for the sperm of sexual males by means of higher energetic resources stored in the hepatopancreas (revealing also feeding performance in gynogenetic females). The better condition status of gynogenetic females of gibel carp when compared to the sexual form was previously revealed according to the total protein concentration in blood (Vetešník et al. 2013). However, a study of gibel carp focussed on seasonal changes in vigour-related traits in gynogenetic and sexual forms (Šimková et al. 2015) revealed differences only between males and females for somatic condition and gonado-somatic index; no obvious differences between sexual and gynogenetic females for these indexes or even in levels of oestradiol were shown, indicating that similar amounts of energy are invested in condition and reproduction by both sexual and gynogenetic females. Furthermore, Vetešník et al. (2013) highlighted the role of calcium important for egg development in both gynogenetic and sexual females of gibel carp. In contrast, Przybyl et al. (2022) showed differences in sex androgen concentrations between sexual and gynogenetic females of gibel carp, indicating that androgen physiology may contribute to the coexistence of gynogenetic-sexual complexes. Concerning reproductive output, Mee et al. (2013) hypothesized the lower fitness of the gynogenetic form when compared to the sexual one; however, their study revealed a weak fitness advantage for gynogenetic Phoxinus eos-neogaeus. Similar fecundities revealed in gynogenetic and sexual females in gibel carp (Šimková et al. 2015, Przybyl et al. 2022) may suggest that males do not discriminate against asexual females, as was revealed in a study on sexual Chrosomus eos (Cope, 1861) and asexual C. eos-neogaeus (Barron et al. 2016); however, the role of potential mate discrimination has never been investigated in gynogenetic-sexual gibel carp.

In agreement with Šimková et al. (2015), our study revealed higher HSI in gynogenetic females when compared to sexual females and sexual males (and we even showed higher HSI in sexual females when compared to sexual males, Fig. S1), indicating higher energy resources in the liver, higher liver functions, and/or a higher potential to mobilize liver resources in asexually reproducing fish, likely reflecting their advantages in habitats shared with sexuals.

The physiology of asexual organisms (typically polyploids) differs from that of sexual diploids, e.g. by the production of hormones or in the number and size of cells, including erythrocytes (Benfey 1999, van de Pol et al. 2020). The lower number of cells of triploids is compensated by the cell size, as was shown, for example, for erythrocytes in tench, this lower number having, in fact, no effect on the haematocrit value (Svobodová et al. 1998). The lower number of erythrocytes in gynogenetic females of gibel carp was previously documented (Šimková et al. 2015) and confirmed in the present study (Fig. S1). In the present study, we revealed higher number of erythrocytes and haematocrit values in males when compared to both sexual and gynogenetic females, which was previously evidenced by Vetešník et al. (2006) and which may be interpreted as indicating the better functional status of aerobic performance (Nussey et al. 1995, Wood et al. 2019). However, Šimková et al. (2015) documented a higher oxygen-carrying capacity per erythrocyte in gynogenetic females, supporting their ecological tolerance, and revealed a trade-off between the number of erythrocytes and the oxygen-carrying capacity per erythrocyte in the sexual-gynogenetic complex of gibel carp as one of the potential mechanisms contributing to the coexistence of gynogenetic and sexual forms. Šimková et al. (2015) documented low aerobic performance expressed by the heart index as a physiological disadvantage in gynogens, and proposed that a limited heart functional capacity in gynogens, reflecting their aerobic performance, organ aerobic activity, and ability to engage in physical activities, may represent one of the mechanisms balancing the evolutionary costs of sexual reproduction in gibel carp.

The adaptation of glucose metabolism is important for a hypotoxic environment. Species of the genus Carassius have been considered as some of the most hypoxia-tolerant species, exhibiting physiological adaptation that even facilitates anoxic survival (Johnston & Bernard 1983, Bickler & Buck 2007, Vornanen et al. 2009). Gong et al. (2020) analysed differential gene expression in artificially-induced gynogenes and natural sexuals of the cyprinid species blunt snout bream (Megalobrama amblycephala Yih, 1955) and showed differences in the expressions of multiple genes associated with hypoxia response, indicating enhanced hypoxia tolerance in gynogens. Our data on blood glucose do not seem to support higher hypoxia tolerance in gynogens of gibel carp; therefore, further studies focused on hypoxic or anoxic natural environments are necessary to clarify the hypoxia mechanisms in gibel carp. Glucose level is also considered as a measure reflecting the level of environmental stress (e.g. Hoar et al. 1992, Martínez-Porchas et al. 2009). Our data revealed slightly higher values of glucose in sexual males (glucose contributed to the data variability expressed by the second physiological axis); however, no differences between gynogenetic and sexual females were found, which is in line with Vetešník et al. (2013).

Immune components in asexual and sexual gibel carp

Considering immune components related to fish vigour, no obvious differences were reported between the three groups of fish in overall immunity analysis using PCA (Fig. 1) with respect to the index of immunocompetence (spleen-somatic index), nonspecific immunity (lysozyme, leukocyte count, oxidative burst, complement activity), and specific immunity(IgMlevel).Thefirstimmunityaxisresulting from PCA was mostly explained by nonparametric immune variables (with leukocytes and complement activity acting in the opposite direction to oxidative burst and lysozyme), and the second immunity axis was mostly explained by specific immunity measured by IgM, and by the index of immunocompetence, and was also explained by oxidative burst and lysozyme with lysozyme acting in the opposite direction than other three variables. However, using individual data (Fig. S1) we showed high IgM level in gynogenetic females when compared to sexuals, and higher lysozyme and spleen-somatic index in sexual males when compared to both females (Fig. S1). Some weak differences were revealed also for parasite load.

Hakoyama et al. (2001) showed that the high prevalence of Metagonnimus (Trematoda) in Japanese crucian carp (C. auratus) was due to the lower immune activity of the phagocytes (representing nonspecific immunity) in the gynogenetic form. However, Šimková et al. (2015) did not reveal differences in nonspecific immunity between gynogenetic and sexual gibel carp, though they documented a difference in specific immunity between the two reproduction forms consistent across different seasons, and proposed that different specific immune responses may contribute to the coexistence of gynogenetic and sexual gibel carp. However, their study did not analyse the level of parasite load for a given sample of gynogenetic and sexual gibel carp, which should potentially induce a higher production of IgM in gynogenetic females. Šimková et al. (2013) focussed on the genetic diversity of the MHC complex and the level of parasite infection in gynogenetic and sexual gibel carp. They revealed that the most common MHC genotypes of gynogenetic females harboured the highest parasite species load (expressed by species richness or abundance), which again supports the role of specific immunity in the coexistence of gynogenetic and sexual forms in the same habitats, and is also in line with the prediction of the Red Queen hypothesis, i.e. that a high level of parasite infection in common asexual clones favours genetically diverse sexuals and promotes the short-term coexistence of sexual and asexual populations. In fact, some studies suggest that asexual fish suffer from a higher parasite load than sexual fish (e.g. Lively et al. 1990, Hakoyama et al. 2001, Mee & Rowe 2006), while others showed no difference in parasite load between gynogens and sexuals (Weeks 1996, Tobler & Schlupp 2005).

Links between parasite load and immune or physiological components in asexual and sexual gibel carp

Finally, we investigated the links between parasite load (expressed by abundance, species richness, and parasite assemblage composition) and immune or physiological components using two approaches: 1) to reveal inter-group differences in the effect of the predictors, and then 2) to investigate which of the predictors explained best the variability of a response variable in each fish group (here, we inferred the inter-group differences based on comparison of the trends observed). Using both approaches, we revealed that abundance and parasite assemblage composition representing two different measures of parasite load were affected by gibel carp immunity and physiology, whilst parasite species richness was almost unaffected by the immunity and physiology of gibel carp. Basically, we revealed different effects of IM1 on parasite abundance and parasite assemblage composition (Table 2 and Table S1), with IM1 explained by nonspecific immunity, i.e. IM1 was positively correlated with leukocyte count and complement activity and negatively with oxidative burst and lysozyme. Using individual immune data analysis, lysozyme seems to be the main immune variable contributing to the difference in nonspecific immunity among fish group (Fig. S2). Analysing each fish group separately, sexuals and gynogens exhibited different associations between immunity components and parasite abundance. Whilst an increase in parasite abundance induced the activation of nonspecific immunity (associated with IM1, Table 1) in sexual fish, parasite abundance in gynogens was associated with specific immunity and the index of immunocompetence (having higher contribution to IM2, Table 1, Fig. 3). Specifically, gynogenetic females, with a high induced specific response and smaller investment in immunocompetence (Fig. S2), were more susceptible to parasites. Šimková et al. (2013) revealed a high parasite species richness and a high total parasite abundance in the gynogenetic form of gibel carp with low MHC variability (as a measure of adaptive immune response). Parasite species richness in gynogenetic females was negatively associated with physiology components representing glucose level and haematocrit (measures of organismal stress induced by environmental factors and oxygen carrying capacity, see above) mostly contributed to the PH2 (Table 1, Fig. 1), indicating that, from the point of view of total parasite abundance, parasites may represent a more important stressor for genetically homogenous gynogenetic fish when compared to highly genetically variable sexuals.

Parasite assemblage composition reflecting parasite interactions was affected by both immunity and physiology. This was clearly evidenced by PERMANOVA (Table 2, Table S2) performed to test the effects of the set of predictors on the overall assemblage structure in a multidimensional space. Considering parasite assemblage composition (reflecting the presence of parasite species or parasite groups in parasite assemblages), our results focused on the analyses of individual groups (Table 3, Table S2) suggest that immune and physiological components in sexuals are more affected by some of the parasite species or parasite groups than in gynogenetic females. In this respect, obvious differences in the associations between parasites and immunity or physiology were documented between genders of sexual gibel carp, while fewer differences were found between gynogenetic and sexual females, i.e. IM1 representing the variables of nonspecific immunity and PH2 representing haematocrit and glucose level respectively were affected by some parasite groups (mostly gill and fin monogeneans of Gyrodactylus) only in sexual females. The associations between parasite assemblage composition and immunity and physiology (expressed by IM2, PH1 and PH2) in males (Fig. 6) may specifically indicate the role of male androgens, mainly testosterone. Specifically, negative associations between PH1 (reflecting mainly body condition, reproductive performance, energy performance, and oxygen-carrying capacity performance expressed by erythrocyte count) and parasitism in males provides evidence that multiple parasitism, estimated by parasite assemblage composition, is costly for male condition and performance. In contrast, the fact that no associations between parasite assemblage composition and immunity or physiology (expressed by IM and PH) were revealed in gynogenetic females (Figs. 5, 6) may suggest their advantage when being infected by multiple dominant parasite species (in contrast to the situation of overall parasite abundance).

The investigation of parasite selection and immune and physiological constraints as potential mechanisms contributing to the coexistence of gynogenetic and sexual forms of gibel carp (and potentially also other asexual-sexual complexes) indicates some variability in parasite load and variability in host immune/physiological responses in gynogens and sexuals. This pattern may partially be driven by the genetic structure of fish populations in regard to adaptive genetic markers which are under different environmental selection. However, some of the observed patterns in immunity and physiology (mainly reflecting the inter-group differences in gonado-somatic index, the hepato-somatic index, erythrocyte count, glucose level and haematocrit level) seem to contribute to the coexistence of asexual and sexual forms. Finally, we suggest that the coexistence of asexual and sexual forms in fish is locality-dependent (habitat-dependent), and that various evolutionary, ecological, and behavioural processes and mechanisms and their multiple interacting effects can facilitate/stabilize the real coexistence of asexuals and sexuals on different temporal and spatial scales.