Vigour-related traits of gibel carp (Carassius gibelio): do they represent reproduction-associated costs facilitating the coexistence of asexual and sexual forms?

Abstract. The coexistence of sexual and asexual lineages is rarely documented in vertebrates. Various ecological and evolutionary mechanisms have been proposed to explain their stable coexistence. Gibel carp (Carassius gibelio), a highly invasive freshwater fish in Europe, is one such example of a sexual and asexual vertebrate, combining gynogenesis (sperm-dependent parthenogenesis) and sexual reproduction. In this study, we focused on vigour-related traits in gynogenetic and sexual forms of gibel carp coexisting in the same habitat, to reveal whether there is a link between parasite load and vigour-related traits reflecting the potential advantage of one reproductive form over another, which may eventually facilitate the obligatory coexistence of sperm-dependent gynogenetic females with sexual males. Using physiological parameters (indexes of somatic condition, energetic condition, and reproductive performance; glucose levels; and erythrocytes-related variables), diploid sexual males, diploid sexual females, and triploid gynogenetic females were found to be differentiated from one another. However, based on immune variables representing innate immunity, specific immunity, and index of immunocompetence, mostly sexual males were found to be weakly differentiated from both groups of females. We revealed different patterns of interactions between parasite abundance and immune variables between sexual and gynogenetic forms. Using parasite assemblage composition, different relationships between parasite assemblage and immunity or physiology in sexual males and sexual females were evidenced, potentially related to male reproduction biology. In contrast, gynogenetic females exhibited the advantage that parasite assemblage composition did not affect their immunity or physiology. Our study suggests that reproduction mode-associated costs of physiology and immunity may facilitate the coexistence of the asexual-sexual complex. We highlight that multiple ecological processes and evolutionary mechanisms contribute to the coexistence of asexual and sexual gibel carp.


Introduction
Asexual reproduction is a rare phenomenon reported in a small number of vertebrates, mostly in fishes, amphibians and reptiles (Dawley 1989, Schmidt et al. 2011).However, if it occurs, the asexual vertebrate lineage usually coexists with the sexual one in natural habitats, which, in fish, has been documented for Poeciliidae, Atherinidae, Cyprinidae and Cobitidae (Lamatsch & Stöck 2009, Lampert 2009).
The coexistence of asexual and sexual forms/species of fish is related to the fact that asexual reproduction is achieved by gynogenesis or hybridogenesis.In the case of gynogenesis, unisexual females use the sperm of conspecific or sympatric closely-related sexual males; however, sperm only triggers embryogenesis without incorporating the paternal genome, which results in the production of clonal all-female offspring (Dawley 1989, Vrijenhoek 1994, Avise 2008).In the case of hybridogenesis, unisexual females mate with sexual males; however, the eggs are produced without recombination and consist exclusively of the maternal genome because premeiotic cell division excludes the paternal genome (i.e. the products of hybridogenesis are hemiclones) (Dawley 1989, Avise 2008).
Gibel carp (Carassius gibelio Bloch, 1782), a member of the Carassius auratus complex, is one of the few fish representatives and the only cyprinid combining gynogenesis and sexual reproduction.Both reproductive forms of gibel carp have been documented in natural conditions, often coexisting in the same habitats (Šimková et al. 2013, Fuad et al. 2021, Przybyl et al. 2022).Asexual reproduction, i.e. gynogenesis, an extended period of reproduction, and a high degree of ecological tolerance are considered to be responsible for gibel carp's successful and rapid invasion into novel habitats (García-Berthou et al. 2005, Przybyl et al. 2022).The impact of gibel carp on freshwater ecosystems has been widely evidenced, especially in the decline of native crucian carp populations and the degradation of the habitats of this native species (Paulovits et al. 1998, Navodaru et al. 2002, Gaygusuz et al. 2007).In the Czech Republic, four mitochondrial DNA lines of the C. auratus complex have been documented -C.gibelio, C. auratus, C. langsdorfii and the so-called M line, with C. gibelio representing the dominant line of the C. auratus complex (Papoušek et al. 2008(Papoušek et al. , Šimková et al. 2015)).A recent study of gibel carp from the River Dyje (in the southeastern part of the Czech Republic) revealed that 96% of sampled C. auratus complex specimens belonged to the C. gibelio line (Pakosta et al. 2018).The highly invasive gibel carp is currently documented across most countries of Eurasia (Kottelat & Freyhof 2007) and represents a severe risk to native fish diversity.Concerning its invasion history in the Czech Republic, the first specimens entered the hydrologic system by migration from the River Danube in 1975 (Lusk et al. 1977).All-female populations of C. gibelio were recorded up to 1992 when the first males were found in natural habitats.A few years later, C. gibelio formed stable mixed sexual-asexual populations composed of sexually reproducing diploid specimens and gynogenetically reproducing typical triploid females and rare triploid and tetraploid males (Lusková et al. 2004).
Concerning gibel carp, the different costs relating to physiological and condition-related traits for gynogenetic and sexual forms of gibel carp have been hypothesized to promote their coexistence (Vetešník et al. 2013, Šimková et al. 2015, Przybyl et al. 2022).For example, lower aerobic performance was found in gynogenetic gibel carp, suggesting that this physiological disadvantage potentially balances the costs of sexual reproduction (Šimková et al. 2015).In addition, differences in immune performance between gynogenetic females and the sexual form of gibel carp were shown (Šimková et al. 2015), which may be potentially explained by the different parasite load (Šimková et al. 2013).According to the Red Queen hypothesis (Van Valen 1973, Jaenike 1978, Bell 1982), asexuals with reduced genetic diversity are an evolutionary target of parasite attack, and, therefore, the higher parasite load in asexuals when compared to sexuals may represent one of the mechanisms stabilizing the coexistence of sexuals and asexuals.The role of parasites in stabilizing the coexistence of asexual and sexual forms of Japanese crucian carp (C.auratus Linnaeus, 1758) was suggested by Hakoyama & Iwasa (2004).The Red Queen hypothesis was tested using a field study on a population of gibel carp, where the most common MHC (major histocompatibility complex) genotype of gynogenetic females exhibited a higher parasite load compared to rare gynogenetic genotypes or sexual genotypes (Šimková et al. 2013).However, Pakosta et al. (2018) highlighted that the parasite load in gibel carp exhibits temporal and spatial variability and is probably driven by multiple abiotic and/or biotic factors rather than Red Queen coevolutionary host-parasite dynamics.
In line with previous studies investigating the ecological factors and evolutionary mechanisms stabilizing the coexistence of gynogenetic and sexual forms of gibel carp, the aims of the present study were 1) to analyse vigour-related traits, i.e. physiological and immune traits, in gynogenetic and sexual forms of gibel carp coexisting in the same habitat, and 2) to reveal whether there is a link between vigourrelated traits and parasite load (measured by parasite abundance, parasite species richness and parasite assemblage composition), reflecting a potential disadvantage in using asexual reproduction relative to the sexual reproduction, which may eventually facilitate/stabilize the coexistence of gynogenetic and sexual forms of gibel carp in the same habitats.

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 2 ) 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 - ) 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 timesdiluted 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 luminolenhanced 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 ).

Coexistence of gynogenetic and sexual gibel carp
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 (logtransformed), 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 (logtransformed), 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 posthoc 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).
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 2 , 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%.

Results
Sequencing the control region of mitochondrial DNA (D-loop) revealed that all specimens sampled during the four-year study belonged to the C. auratus complex.For the present study, only specimens determined as the C. gibelio mitochondrial line were selected (corresponding to 96% of the specimens examined).Both the sexual form (63% of specimens corresponding to diploid males and females) and the gynogenetic form (37% of specimens corresponding to triploid females and 1% to triploid males) were found.
The following parasite species within the following taxonomic unit were recognized in gibel carp and used in the calculation of parasite species richness and abundance: Protozoa -I.

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 (posthoc 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  ) 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 nonsignificant after correction for multiple testing.
Coexistence of gynogenetic and sexual gibel carp 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 nonsignificant after correction for multiple testing.
Coexistence of gynogenetic and sexual gibel carp 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).
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, posthoc tests, P = 0.019 and 0.045).
The first NMDS axis (NMDS1) was primarily associated with Gyrodactylus spp.and S. petruschewskii (positive values), whilst the second NMDS axis 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.Coexistence of gynogenetic and sexual gibel carp (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).

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 nonsignificant 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).

Discussion
The persistence and coexistence of asexual and sexual forms of fish within a single species or species 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 2 (Table S2).  ) and the first NMDS axis (a proxy for parasite assemblage composition) 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.

Response
Fig. 6.Relationship between the first two axes of immunity (IM1, IM2) and physiology PCA (PH1, PH2) and the second NMDS axis (a proxy for parasite assemblage composition) 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.
Coexistence of gynogenetic and sexual gibel carp complex occupying the same aquatic habitats have been explained by various ecological processes and evolutionary mechanisms.The necessity of the coexistence of asexual and sexual forms is even most pronounced in the case when asexuals reproduce by gynogenesis, i.e. by sperm-dependent parthenogenesis, as the persistence of these asexuals is directly interlinked with the presence of sexual sperm donors (i.e.males of conspecifics or closelyrelated species).Therefore, the present study was focused explicitly on selected vigour-related traits expressed by immune variables and physiological variables in gynogenetic and sexual forms of gibel carp and on parasite load measured by parasite abundance, parasite species richness and parasite assemblage composition.
From the evolutionary point of view, because of the two-fold costs of sexual reproduction (Maynard Smith 1978, Bell 1982), the sexual form risks being outcompeted by asexual forms that use the same resources and produce eggs at twice the rate.Such sexassociated costs should be compensated by certain long-term evolutionary or short-term ecological and behavioural disadvantages faced by asexuals, mainly mating discrimination against asexuals (Mee & Otto 2010, Riesch et al. 2012), differential susceptibility to food stress (Tobler & Schlupp 2010), a fitness disadvantage among asexuals (Mee et al. 2013), parasite selection of the most common asexual genotype (Hamilton et al. 1990, Mee & Rowe 2006, Šimková et al. 2013), and the limited physiological performance of asexuals (Šimková et al. 2015(Šimková et al. , Przybyl et al. 2022)).Concerning gibel carp forming populations composed of coexisting gynogenetic females and sexual male and female specimens, we can hypothesise that the stable coexistence of asexuals and sexuals should be facilitated by mechanisms compensating the two-fold costs faced by sexuals and mechanisms regulating the potential competition between gynogenetic and sexual females for the sperm of sexual males.Therefore, we analysed the immune and physiological components of fish vigour and parasite load in three fish groups -gynogenetic triploid females, sexual diploid females, and sexual diploid males.To analyse the overall differences in immunity and physiology among three fish groups, we applied a PCA approach to collapse the number of predictors (following Rohlenová et al. 2011), and we supplemented this approach with the analysis of individual physiological and immune components by demonstrating the consistent results using both approaches.In the present study, we hypothesized that some vigour-related traits may represent potential advantages or constraints for asexuals which contribute differentially to the inevitable coexistence of sperm-dependent gynogens and their sexual partners.

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.
Coexistence of gynogenetic and sexual gibel carp 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(Šimková et al. , 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 oxygencarrying 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 (IgM level).The first immunity axis resulting 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 Coexistence of gynogenetic and sexual gibel carp 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 shortterm 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 Coexistence of gynogenetic and sexual gibel carp 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.

Fig. 1 .
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. 3 .
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 nonsignificant after correction for multiple testing.

Fig. 5 .
Fig.5.Relationship between the first two axes of immunity (IM1, IM2) and physiology PCA (PH1, PH2) and the first NMDS axis (a proxy for parasite assemblage composition) 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 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.
Coexistence of gynogenetic and sexual gibel carp