Host Genetics and the Skin Microbiome Independently Predict Parasite Resistance

TL;DR

Host genetics and skin microbiome each independently influence resistance to parasites in guppies, with both factors together providing better predictions of infection severity.

Contribution

Demonstrates that host genetics and microbiome independently predict parasite resistance, with combined effects explaining infection severity better than either alone.

Findings

Resistant guppies had lower infection severity in females, while susceptible guppies had higher severity.

Three distinct bacterial community types were identified, independent of selection lines, that explained infection severity.

Microbiome and host genetics together explain infection severity better than either alone.

Abstract

Host responses to parasite infection involve several interacting systems. Host genetics determine much of the response, but it is increasingly clear that the host‐associated microbiome also plays a role. Host genetically determined systems and the microbiome can also interact; for example, the microbiome can modulate the immune response, and vice versa. However, it remains unclear how such interactions between the host immune system and the microbiome may influence the host's overall response to parasites. To investigate how host genetics and the microbiome interact to shape responses to parasites, we imposed truncation selection on Trinidadian guppies (
Poecilia reticulata
) for low and high resistance to the specialist ectoparasite Gyrodactylus turnbulli. After 3–6 generations of breeding without parasites, we sampled the skin‐associated microbiome and infected fish from each line. We applied Dirichlet Multinomial Modeling (DMM) machine‐learning to identify bacterial community types across lines and evaluated how selection line and community type explained variations in infection severity. Our findings showed that among females, the resistant line had significantly lower infection severity, while the susceptible line had higher infection severity. Among males, only the susceptible line experienced higher infection severity compared to the other lines. Line did not explain skin microbial diversity, structure or composition. Our DMM analysis revealed three distinct bacterial community types, independent of artificial selection lines, which explained just as much variation in infection load as selection line. Overall, we found that the microbiome and host genetics independently predict infection severity, highlighting the microbiome's active role in host–parasite interactions.

Host genetic selection and skin microbiome independently predict infection severity of ectoparasites. Furthermore, both factors together explain infection significantly better than either alone.

Figures

Experimental design. In Step 1, we infected 200 F1 lab‐bred offspring of wild caught Trinidadian guppies and imposed one round of truncation selection. Sixty fish founded each line (30 females and 30 males), which we then bred in the absence of parasites for 3–6 generations. In Step 2, we took skin microbiome samples from a subset of fish before experimentally infecting 87 individuals. Fish were weighed and measured before and after infection. Parasites were counted under the dissecting scope every other day for 18–19 days.

Artificial selection lines significantly predicted G. turnbulli infection severity in Trinidadian guppies. (A) Infection severity (worm days) reflects the area under the infection load curve. Female guppies (n = 44) in susceptible and resistant lines differed from the control line. Male guppies (n = 43) only showed differences in the susceptible line. Box plots illustrate median, interquartile range, and values within 1.5× the interquartile range of infection severity. (B) Parasite counts over time by line and fish sex. The number of parasites is modeled as a non‐linear function of day, with separate smooths for each combination of line and sex, using a Generalized Additive Mixed Model. Points are raw data from infected individuals, and the solid lines represent the model's predicted values.

Artificial selection lines also predicted relative tissue‐specific tolerance in Trinidadian guppies. Tolerance measured as per‐parasite mass change across infection minus max tolerance. Female guppies (n = 44) in resistant lines significantly differed from control and susceptible lines. Male guppies (n = 43) did not show significant differences. Box plots illustrate the same metrics as in Figure 2.

The distribution of skin bacterial community types across guppy selection lines showed no significant association. (A) A mosaic plot presents bacterial community types (A–C) across control (c), resistant (r), and susceptible (s) lines, with no significant distribution differences. (B) PCoA analysis (n = 50) indicated bacterial community type was significant in variation, but no significant clustering was detected for selection lines.

Selection lines and bacterial community types predicted G. turnbulli infection severity. (A) Lines significantly predicted severity in swabbed fish (n = 31). (B) Community type A had lower infection severity than types B and C. No interactions between community type and selection line were detected. A combined model of line and bacterial types explained infection severity better than either factor alone.

Full text

Introduction

1

Hosts use multiple strategies to defend themselves against parasites and pathogens (parasites hereafter). Host genetic background is commonly credited with much of the defense against parasites through the immune system, behavior, or other genetic/protein alterations, for example, thalassemia (Labro 2022). Recent research has highlighted that the host‐associated microbial communities, or “microbiome,” may also influence parasite defense indirectly through the host immune system or directly via niche occupation or secreting anti‐parasitic compounds (Bernardo‐Cravo et al. 2020), and may also exhibit heritable components (A. H. Morris and Bohannan 2024). Indeed, the host‐associated microbiome can dramatically impact host survival and defense against parasites. However, the question of whether the microbiome plays an active role in host–parasite infection or merely reflects the host's health/immune status is commonly posed (Trevelline et al. 2020). While there is abundant evidence that the microbial community can have a direct or indirect role in host–parasite interactions (Becker et al. 2009; Nelson and May 2020), there is also evidence that the microbiome is merely a bystander in some diseases and a reflection of the host's genetics, diet, and general health (Ni et al. 2017).

Host genetics and the microbiome interact throughout host development, with complex feedback complicating the separation of their roles in traits like parasite defense. On one hand, hosts exert control over their symbionts, shaping their associated microbiomes, via immune responses, barrier function, and physiological regulation (Wilde et al. 2024). On the other hand, the host immune system relies on early microbial colonization for proper development and ultimately shapes immune cell composition (Amenyogbe et al. 2017). The host's genetic makeup therefore influences microbiome composition and may lead to the selection of specific microbial communities that can be inherited across generations (Bordenstein and Theis 2015; A. H. Morris and Bohannan 2024). However, distinguishing host genetic effects on the microbiome from factors such as priority effects, drift, and other stochastic elements is an ongoing challenge (Debray et al. 2022), and the influence of external factors such as environment and temporal changes in shaping host‐associated microbiomes must also be carefully considered. Direct manipulation of the host genetics in traits like parasite resistance, where the microbiome likely plays a role, thus provides a useful tool in disentangling host genetic effects from microbial ones.

Both the host‐associated microbiome and the immune system commonly differ between host sexes (Markle and Fish 2014). These differences are likely due to other correlated differences, for example, high testosterone levels in the blood both suppress the immune system and significantly alter gut microbial communities (Santos‐Marcos et al. 2023), and both the immune system and microbiome change during pregnancy and menstruation (Chen et al. 2023; Gorczyca et al. 2022). Given that host sexes across taxa also differ in their response to parasites (Morrow 2015), it is essential to test for sex‐specific interactions between genetic resistance factors, the microbiome, and parasites.

Trinidadian guppies ( Poecilia reticulata ) and their specialist ectoparasite Gyrodactylus turnbulli provide a rich context to study host–parasite dynamics and interactions with the microbiome. Gyrodactylus spp. feed and reproduce (both sexually and asexually, generation time 24 h) on host skin. Guppies show a large variation in parasite resistance across populations, sexes, and ecological contexts (Dargent et al. 2016; K. P. Phillips et al. 2021; Stephenson et al. 2015; Van Oosterhout et al. 2003). Guppy resistance is at least partly genetically determined and heritable (Dargent et al. 2016; Karl P. Phillips et al. 2018; K. P. Phillips et al. 2021; Van Oosterhout et al. 2003; Weiler et al. 2025). Innate and acquired immunity are important for G. turnbulli resistance; however, the precise mechanisms involved in guppy parasite resistance remain unclear (Cable and van Oosterhout 2007; Mohammed et al. 2016). Research into genetically based resistance to G. turnbulli has largely focused on the major histocompatibility complex (MHC) (Fraser et al. 2010; Fraser and Neff 2009; Karl P. Phillips et al. 2018). Notably, skin bacterial alpha diversity correlates with MHC heterozygosity in frogs (Cortazar‐Chinarro et al. 2024), and probiotic supplementation results in significant upregulation of MHC genes in fish (Tanpichai et al. 2025). However, other immune pathways have been indicated in gyrodactylid infection that are also linked to microbial communities in host mucus (Konczal et al. 2020; Ivanov et al. 2008). Similar associations and covariation between host genetics and the skin microbiome of fish would set up a clear avenue for the skin microbiome and host genetics to interactively shape parasite resistance phenotypes in our system.

Our primary objective was to understand how host genetic background and the skin microbiome interact or shape each other as predictors of parasite resistance and susceptibility. To investigate this, we imposed truncation selection for host resistance and susceptibility to the parasite, quantified as the mean number of parasites over infection. After 3–6 generations of breeding, we infected offspring of the truncation selection with the addition of sampling the skin microbiome before infection. To explore alternative explanations of selection patterns, we also examined the levels of tolerance (ability to minimize the fitness losses per parasite) within our lines. Through these investigations, we found that host genetics and the microbiome significantly and independently predicted infection severity, with parasite resistance and tolerance showing sex‐specific responses to selection.

Materials and Methods

2

Data and Sample Collection

2.1

Truncation Selection to Establish Guppy Lines That Differ in Genetically Based Resistance to Gyrodactylus Turnbulli

2.1.1

We used laboratory‐bred descendants of wild guppies from a high‐predation population on the Caura River, Trinidad (UTM: 20 P 679527.7 m E, 1180376.4 m N; elevation 112 m). Wild guppies (n ~600) were transported to Cardiff University in June 2012 (Cefas APB authorization number CW054‐D‐187A) and treated prophylactically with Binox (Nitrofurazone). The site was previously confirmed Gyrodactylus spp.‐free by previous surveys, and no Gyrodactylus spp. were found in a subset of 80 fish examined. We collected fin clips from 66 parental‐generation fish for microsatellite genotyping. Fish were housed in 70 L mixed‐sex tanks under a 12 h light:12 h dark cycle at 24°C ± 1°C and fed daily with Aquarian tropical fish flakes, supplemented with Artemia and bloodworm. Fry were transferred to juvenile tanks soon after birth, checked for ectoparasitic infection, and reared until sex determination (6–8 weeks). They were then kept in sex‐specific tanks for at least another month before use, ensuring all fish were sexually mature and unmated. Gyrodactylus turnbulli were from an established isogenic culture in the lab (Gt3), maintained on commercially sourced guppies (“culture fish”).

We experimentally infected 100 male and 100 female F1 guppies to evaluate their resistance and select founders for our “resistant” and “susceptible” lines. To initiate infections, a heavily infected culture fish was euthanized with MS222 (0.02% tricaine methanesulfonate) and placed next to an anesthetized experimental fish until two G. turnbulli parasites were transmitted, as observed under a dissecting scope. Experimental fish were housed individually in 1 L tanks, and parasites were counted every other day on anesthetized fish. After 9 days, all fish were treated with Levamisole and screened under anesthetic several times to confirm parasite clearance. Fin clips were collected for the microsatellite genotyping described below. The 30 males and 30 females with the lowest mean number of worms (mean ± SEM: 3.38 ± 0.28; range: 0–8) founded the “resistant” line, and those with the highest mean number of worms (mean ± SEM: 21.3 ± 0.75; range: 8–36) founded the “susceptible” line, and randomly selected uninfected F1 fish (30 of each sex) founded the control line. Each line was allowed to breed freely in separate 120 L tanks. Tanks were monitored for fry, which were moved to F2 rearing tanks. The F2 generation was split into two subpopulations to mitigate potential effects of the tank and other stochastic processes post‐selection (named “A” and “B” within each of the susceptible, resistant, and control lines), maintaining 60 fish per tank. This population size was kept consistent across the six subpopulations, with new tanks set up for each subsequent generation as above. During preliminary analysis, we confirmed that the A and B subpopulations did not differ in either infection severity or parasite tolerance within each line (see Supporting Information: Section 3.1) and therefore combined them for the subsequent analyses described below.

Lines were bred in the absence of the parasite under the laboratory conditions described above in Cardiff for ~3 generations, then shipped to Zürich in 2017 and bred for a further ~3 generations. In Zurich, guppies were housed in mixed‐sex groups at densities of 1–2 fish per L in 4.5 L tanks in a recirculating system at 24°C ± 1°C, on a 12 h light:12 h dark lighting schedule (overhead fluorescent lighting) and fed daily on commercial flake food supplemented with Artemia. Each subpopulation was therefore housed across at least 10 separate tanks, again with multiple new tanks set up for each subsequent generation 3–6. The fish were held and bred under these conditions from August 2015 until their experimental infections in August–November 2017.

A Note on Replication in Our Study Design

2.1.2

As described above, we imposed one round of truncation selection for resistance and susceptibility on one group of F1 fish (Experimental design: Figure 1). Typically, selection experiments using guppies to characterize the genes underlying phenotypic variation (e.g., van Wijk et al. 2013), or evolutionary change in correlated traits (e.g., Bartuseviciute et al. 2022; Kotrschal et al. 2013) replicate this process three times because the scale of inference is at the level of the selection event: to draw conclusions about how underlying genes or correlated traits respond to selection, the selection event should be replicated (Kawecki et al. 2012).

Experimental design. In Step 1, we infected 200 F1 lab‐bred offspring of wild caught Trinidadian guppies and imposed one round of truncation selection. Sixty fish founded each line (30 females and 30 males), which we then bred in the absence of parasites for 3–6 generations. In Step 2, we took skin microbiome samples from a subset of fish before experimentally infecting 87 individuals. Fish were weighed and measured before and after infection. Parasites were counted under the dissecting scope every other day for 18–19 days.

By contrast, we do not identify genes, or focus on correlated evolutionary trait changes: instead, we use truncation selection to maximize the variation in response to infection in our focal population, and to compare the role of host genetics and associated microbiomes in infection phenotypes. That resistance to Gyrodactylus spp. is heritable is not our focus, and has already been established several times (Gotanda et al. 2013; Madhavi and Anderson 1985; Karl P. Phillips et al. 2018; K. P. Phillips et al. 2021; Weiler et al. 2025). To ensure that the difference between our lines was due to selection, and thus underlying genetic differences, we kept consistent population sizes across generations, split each line into two subpopulations, housed fish across several replicate tanks per line, bred the fish for at least 3 generations in the absence of parasites, and used genotyping to eliminate the role of inbreeding differences between them. Because of these features of our design, we can exclude several potential explanations of our results: the effects of the tank or other environmental conditions; stochastic differences post selection; inbreeding; or maternal effects.

Our selection process thus produced individual guppies that differed in genetically based resistance to our focal parasite; here we show that the microbiomes associated with these individuals contributed independently to their infection phenotype. Our main conclusion, that the host‐associated microbiome predicts infection phenotype independently of and just as strongly as host genetics, is therefore robustly supported by our study design.

Microbiome Sampling Protocol

2.1.3

To evaluate how host genetics and the microbiome combine to affect host response to parasites, we collected guppies 3–6 generations post truncation selection from the breeding tanks. We collected skin microbiome samples from 50 fish split between the infection and sham‐infection treatment groups described below: control (infected n = 10, sham n = 5), resistant (infected n = 12, sham n = 6), and susceptible (infected n = 9, sham n = 8). We used the methods previously published to collect and inventory the microbial communities on the fish skin (Kramp et al. 2022). Briefly, fish were washed twice with two separate 20 mL aliquots of sterile water to ensure that the sample primarily consisted of skin‐associated microbes rather than transient microbes (Kueneman et al. 2014; Kohl 2017). We then swabbed the body surface with a sterile cotton swab for 30 s. Between individual fish, all glassware was sterilized with 70% ethanol. This methodology may not have removed all the DNA of dead or transient bacteria (doing so is rarely possible). However, any contamination and resulting homogenization of the microbiome samples across fish would have reduced our ability to detect bacterial type clusters.

Experimental Infection

2.1.4

We infected 87 guppies 3–6 generations post truncation selection with G. turnbulli in five batches, following the same procedure as for the F1 generation (described above). Gyrodactylus turnbulli for this infection were collected from a local pet shop in June 2017. An isogenic strain was initiated by infecting a parasite‐naïve ornamental guppy with a single parasite individual. This strain was subsequently maintained on parasite‐naïve culture fish and identified to species level using molecular methods. We aimed to infect each with 2 individual parasites, but in practice the infecting dose varied (mean ± SEM: 2.94 ± 0.2). Fourteen guppies from each line were “sham” infected (following the same infection protocol mentioned above without exposure to parasites) and were only used in microbiome and genetic analysis. To minimize over‐handling, fish were only weighed and measured on day 0 of infection and their last day of infection.

What Is the Magnitude of Host Genetic Contribution to Guppy Resistance to Gyrodactylus Turnbulli?

2.2

Statistical Analysis

2.2.1

We used R v.4.3.2 (R Core Team 2024) for these analyses and provide our code, output, and validation in the Supporting Information. We used ggplot2 (Wickham 2016) to plot all the figures. For all models, we assessed collinearity among predictor variables using the pairs function (Zuur et al. 2009), validated model fits using the DHARMa package (Hartig 2022), and used Type II Wald chi‐square tests to test for significance using the Anova() function from the “car” package (Fox and Weisberg 2019).

For our test metric of resistance and susceptibility as “infection severity”—the area under the curve of parasite load across the 18–19 days of infection. This metric provides a comprehensive summary of each individual's infection intensity over time, and because it incorporates multiple observations, it is less sensitive to observation errors compared to more commonly used metrics like day of maximum parasite load or maximum parasite count. Accordingly, it has been widely and effectively used in this system (Stephenson et al. 2017, 2018; Karl P. Phillips et al. 2018) as well as in others (Adelman et al. 2013).

To test response to selection, we used infection severity as the response variable in a generalized linear mixed model (GLMM) in the glmmTMB package, v. 1.1.2.3 (Brooks et al. 2017) (Gamma error family and logit link function). As fixed effects, we included selection line (resistant n = 30, control n = 28, susceptible n = 29), fish sex, length (pre‐infection, using residuals from a regression on sex to control for sex differences), body condition (scaled mass index; Peig and Green 2009, 2010), generation, number of parasites establishing the infection (dose), and trial date (Julian date). Length and body condition residuals were calculated with all fish including sham infected fish to obtain a more accurate distribution. We also included the interaction between sex and line.

We quantified tolerance as the host's relative per‐parasite percentage mass change during infection, calculated by dividing the percentage mass change from day 0 to day 18–19 by the maximum number of worms on the fish, relative to the highest observed percentage change. To assess the impact of selection on tolerance, we used tolerance as a response variable in a GLMM, including the fixed effects of line, generation, sex, length, dose, and Julian date. We excluded scaled mass index due to its overlap with mass change. Furthermore, we did not include infection severity due to high collinearity with other fixed effects within the model, model fit issues, and the non‐linear relationship between infection severity and ppmass.

In a different model, we tested for a tradeoff between resistance and tolerance with a quadratic polynomial regression model to account for non‐linear effects. We used tolerance as the response variable and log‐transformed infection severity (our resistance metric), sex, their interaction, length, and generation as fixed effects.

Calculating Heritability

2.2.2

We calculated narrow‐sense heritability (assuming that the response to selection is primarily due to additive genetic variance) of resistance (using our mean worm metric—the resistance metric we used during truncation selection) and tolerance following the breeder's equation (Lush 1937). For mean worm, we used the average number of parasites between 1 and 9 days of infection, centered to a mean of zero to account for the differences in worm growth between the infections of the F1 and subsequent generations. Tolerance was centered in the same way (Supporting Information: Section 3.8).

P b: Average starting phenotype (mean phenotype of the entire F1 population before selection).

P s: Average selected phenotype (mean phenotype of the selected parents, resistant or susceptible).

P r: Average response phenotype (mean phenotype of the offspring infected population, resistant or susceptible).

h ^2^: Narrow‐sense heritability.

$$[eqn]$$

Do Differences in Inbreeding Contribute to Line Differences?

2.2.3

We evaluated the genetic diversity and variance in inbreeding within the wild‐caught populations, F1, and F3‐F6 generations, to determine if the observed difference between the resistant and susceptible lines could be attributed to differences in inbreeding (control: n = 42; resistant: n = 37; susceptible: n = 36; wild caught: n = 66). We extracted DNA from fin clips with the HotSHOT method: using 50 μL of alkaline lysis reagent and neutralizing solution followed by incubation at 95°C for 30 min (Truett et al. 2000). We used two multiplex assays to amplify eight previously designed microsatellite markers Pret‐69, Pret‐77 (Watanabe et al. 2005), Pre‐9, Pre‐15, Pre‐18 (Paterson et al. 2005), PP‐GATA‐5 (Nater et al. 2008), Pr‐92 (Becher et al. 2002), Hull 9‐1 (Watanabe et al. 2005) at the Genetic Diversity Centre (GDC), ETH Zurich. Samples were genotyped on an ABI Genetic Analyzer 3730 (Applied Biosystems, USA), scored using GeneMarker (Version 2.6.4; Softgenetics), and binned using the package MSatAllele (Alberto 2009) in R (R Core Team 2024). See the Supporting Information for more details (Supporting Information: Section 3.7).

We tested allelic richness (A R), gene diversity (H S), F IS, F ST, and identity disequilibrium (g 2) among our artificial selection lines to evaluate whether the differences in phenotype were attributable to variation in inbreeding and/or genetic drift (Smallbone et al. 2016). We estimated identity disequilibrium (g 2), the correlation of heterozygosity across loci resulting from inbreeding, which is robust in the presence of null alleles (David et al. 2007). Concurrence among population genetic metrics consistent with patterns inbreeding is needed to determine if phenotypic differences among populations are attributable to inbreeding. The fixation index (F ST) was assessed to evaluate the degree of population differentiation potentially caused by genetic drift within lines. Allelic richness (A R) (rarefied to the lowest population sample size), gene diversity (H S) and F IS was calculated in FSTAT 2.39 (Goudet 2001). Significance levels for F IS were assessed in GENEPOP (web version 4.2; Raymond and Rousset 1995) with 5000 iterations followed by a Bonferroni correction for multiple comparisons. Tests for statistically significant differences in allelic richness (A _ R _), gene diversity (H _ S _) and F _ IS _ among groups were conducted in FSTAT with 2000 permutations. Identity disequilibrium (g 2) was calculated in inbreedR (Stoffel et al. 2016). Pairwise population subdivision (F ST) was assessed in adegenet (Jombart et al. 2011).

Does the Host‐Associated Microbiome Contribute to Guppy Resistance to Gyrodactylus Infection?

2.3

Microbial Inventory

2.3.1

We used QIAamp PowerFecal DNA Isolation Kit (Qiagen, Hilden, Germany) to extract DNA, following the manufacturer's protocol. To control for contaminants, we included four control blank samples in extractions and sequencing (Salter et al. 2014). Extracted DNA was shipped on dry ice to the DNA Services Facility at the University of Illinois at Chicago (Chicago, IL, USA). All library preparation, amplification, and sequencing were performed on the Illumina Miseq platform (V4 region of the bacterial 16S rRNA gene), using primers 515F (5′‐GTGCCAGCMGCCGCGGTAA‐3′) and 806R (5′‐GGACTACHVGGGTWTCTAAT‐3′). We filtered the raw sequences for quality, removed chimeric reads, merged forward and reverse reads, and assigned the remaining reads to amplicon sequence variant (ASVs) using the DADA2 pipeline in QIIME2 2019.7 (Bolyen et al. 2019; Callahan et al. 2016). We trimmed each read to 220 base pairs, removed all singleton reads, and used FastTree to build a phylogenetic tree (Price et al. 2010), using the SILVA database classifier release 138 to assign taxonomy (Quast et al. 2013). Untargeted sequences were removed. We removed all non‐bacterial reads (archaea, chloroplast, or mitochondria) and any ASVs detected in all control blank samples prior to analysis with 33 features removed, for complete list (see Dryad Repository) (Salter et al. 2014; Hornung et al. 2019). We rarefied the resulting ASV table to 3654 sequences per sample, which resulted in excluding one sample.

Microbial Alpha and Beta Diversity

2.3.2

We used the DADA2 pipeline to calculate the unweighted and weighted UniFrac distances (Lozupone and Knight 2005) pairwise between all samples to create distance matrices. We used a PERMANOVA from the adonis package in R (McArdle and Anderson 2001) to test for significant differences in microbial community across fish line, sex, length, and body condition. The significance of each term was evaluated sequentially using marginal tests (Type II sums of squares), which account for the variance explained by each predictor after all other predictors are included in the model. We used QIIME2 to calculate the observed number of ASVs, Shannon's diversity index, and Faith's phylogenetic diversity, and compare the mean values of these metrics across line and sex using Pairwise Kruskal–Wallis tests.

Bacterial Community Typing

2.3.3

We used Dirichlet multinomial mixture models (DMM) (Holmes et al. 2012) in the Dirichlet Multinomial v1.16.0 package (Morgan 2023) to assign all samples (n = 50) to three bacterial community types based on the relative abundance of taxa. The DMM model was trained on all microbial features present after decontamination and a relative abundance threshold of 0.01%. To determine the optimal number of clusters (K), we compared the Laplace approximation of the negative log models evidence for values of K, as this fully Bayesian approach to model comparison is superior for this model type and can, under certain circumstances, provide more robust results than AIC or BIC (Holmes et al. 2012; Bishop 2006). The Laplace approximation indicated that the model with K = 3 components had the strongest statistical support (i.e., the highest model evidence) for our dataset (Supporting Information: Section 4.4). All samples were assigned to one of three distinct bacterial community types with a posterior probability of 98% or higher (see Dryad Repository).

Does Bacterial Community Type Predict the Severity of Gyrodactylus Infection?

2.3.4

To test if the bacterial community type significantly predicted infection severity, we used bacterial community type as a fixed effect in the GLMM model used for the selection line resistance analysis described in Section 2.1.2, with the data subset to just swabbed fish.

What Is the Relative Importance of Host Genetics and the Microbiome to Guppy Resistance to Gyrodactylus Infection?

2.4

To compare how well our GLMM models with (a) bacterial community type only, (b) artificial selection line only, and (c) both factors explained variation in infection severity, we used the anova() and summary() functions. All three models shared the same fixed effects, differing only in whether line or bacterial community type or both were included as a fixed effect or the line interaction with sex. Models involving skin microbiome and infection were run on the reduced dataset using only the infected and swabbed fish (n = 31).

Results

3

Host Genetics Contribute to Guppy Resistance to, and Tolerance of, Gyrodactylus Turnbulli Infection

3.1

Resistance Responded to Selection

3.1.1

Selection line was a significant predictor of infection severity (Figure 2; χ ^2^(2) = 37.74, p = 6.4e‐09) Additionally, an interaction between selection line and sex was observed (χ ^2^(2) = 13.51, p = 0.0012). Females had significantly higher infections than males (χ ^2^(1) = 24.29, p = 8.3e‐07) and high starting worm “dose” resulted in higher infection severity (χ ^2^(1) = 5.11, p = 0.024). Post hoc tests to evaluate whether selection line significantly explained infection severity among each sex found that line was still significant in both sexes (female; χ ^2^(2) = 24.33, p = 5.2e‐06 and male; χ ^2^(2) = 17.57, p = 0.00015). Post hoc Tukey tests also revealed that the female guppies from the susceptible line experienced significantly higher infection severity (Cohen's d = 3.49, p < 0.001) and from the resistant line experienced significantly lower infection severity (Cohen's d = 2.74, p = 0.016) than the control. Male guppies from the susceptible line experienced significantly higher infection severity than control line males (Cohen's d = 4.19, p < 0.001), but resistant line males did not differ significantly from control line males (Cohen's d = 1.39, p = 0.15). The heritability of the average number of parasites throughout infection was h ^2^ = 0.79 in the susceptible line, and h ^2^ = 0.41 in the resistant line.

Artificial selection lines significantly predicted G. turnbulli infection severity in Trinidadian guppies. (A) Infection severity (worm days) reflects the area under the infection load curve. Female guppies (n = 44) in susceptible and resistant lines differed from the control line. Male guppies (n = 43) only showed differences in the susceptible line. Box plots illustrate median, interquartile range, and values within 1.5× the interquartile range of infection severity. (B) Parasite counts over time by line and fish sex. The number of parasites is modeled as a non‐linear function of day, with separate smooths for each combination of line and sex, using a Generalized Additive Mixed Model. Points are raw data from infected individuals, and the solid lines represent the model's predicted values.

Tolerance Responded to Selection and Trades Off With Resistance

3.1.2

The selection lines significantly differed in relative tolerance (Figure 3; χ ^2^(2) = 14.56, p = 0.00069) and significantly interacted with sex (χ ^2^(2) = 10.28, p = 0.00586). Additionally, larger fish were more tolerant (χ ^2^(1) = 5.94, p = 0.014). A post hoc test found that selection line was only significant in female guppies (female; χ ^2^(2) = 23.81, p = 6.7e‐06), and smaller females were less tolerant (χ ^2^(2) = 6.96, p = 0.0083). Post hoc Tukey tests revealed that resistant line females were significantly less tolerant than control and susceptible line females (control; Cohen's d = 7.04, p = 0.0003; susceptible; Cohen's d = 8.54, p < 0.00016). However, no fixed effect explained variation in tolerance among males. We found that the bacterial community types did not significantly predict tolerance. The heritability of tolerance was h ^2^ = 0.52 in the susceptible line and h ^2^ = 0.24 in the resistant line.

Artificial selection lines also predicted relative tissue‐specific tolerance in Trinidadian guppies. Tolerance measured as per‐parasite mass change across infection minus max tolerance. Female guppies (n = 44) in resistant lines significantly differed from control and susceptible lines. Male guppies (n = 43) did not show significant differences. Box plots illustrate the same metrics as in Figure 2.

The polynomial regression model revealed a tradeoff between tolerance and resistance at the individual level: higher tolerance was correlated with lower resistance (Figure S1; F = 55.49, p < 7e‐15). No other fixed effects explained significant variation in tolerance.

Selection Lines Did Not Differ in the Extent of Inbreeding or Genetic Drift

3.1.3

We saw a significant reduction in allelic richness (p = 0.005) among all experimental infection groups relative to the wild‐caught founders, likely reflecting a post‐capture reduction in population size (Table S1). However, we found no significant difference in H s (p = 0.09), F IS (p = 0.83), or g 2 (Table S2) among the three lines. Pairwise FST values between wild‐caught and experimental groups were not consistent with population subdivision due to genetic drift (Table S3).

The Host‐Associated Microbiome Is Associated With Guppy Resistance to Gyrodactylus Turnbulli Infection, Independently and to a Similar Degree as Host Genetics

3.2

DMM modeling found that the lowest Laplace approximation model fit was three bacterial clusters (Figure S2). All samples had a posterior probability of 98% or higher. Each bacterial community type contained 10.0%–66.8% of all samples (n = 5–33), and no samples were undefined. The three bacterial community types (A, B, and C) in our fish skin samples were distinct (PERMANOVA: R ^2^ = 0.074, F = 3.33, p = 0.002, Figure 4a). Males and females were balanced among the bacterial community types A and B, and bacterial community types clustered evenly across the artificial selection lines (Figure 4a; Fisher's exact test, p = 0.8). Additionally, we found no interaction between line and bacterial community types in our model (χ ^2^(4) = 3.07, p = 0.546). All skin swab samples were assigned to a community type with the greatest posterior probability. For top contributing bacterial genera, see Figure S3 and Dryad Digital Repository.

The distribution of skin bacterial community types across guppy selection lines showed no significant association. (A) A mosaic plot presents bacterial community types (A–C) across control (c), resistant (r), and susceptible (s) lines, with no significant distribution differences. (B) PCoA analysis (n = 50) indicated bacterial community type was significant in variation, but no significant clustering was detected for selection lines.

The selection lines did not significantly differ in any skin microbial diversity metrics tested. The microbiomes associated with females had significantly higher alpha diversity than those associated with males (Kruskal–Wallis pairwise test, Shannon's index: H = 4.06, q = 0.04; observed ASVs: H = 4.88, q = 0.03).

We found that the bacterial community type A of the skin microbiome on Trinidadian guppies had significantly reduced infection severity compared to bacterial community type B and C (χ ^2^(2) = 15.15, p = 0.0005, Figure 4b). A post hoc test showed that a model including both selection line and bacterial community type is significantly better at explaining variation in infection severity than models with either variable alone (Figure 5; χ ^2^(14) = 13.40, p = 0.0012). Effect sizes for each fixed effect are listed with Supporting Information (Section 5).

Selection lines and bacterial community types predicted G. turnbulli infection severity. (A) Lines significantly predicted severity in swabbed fish (n = 31). (B) Community type A had lower infection severity than types B and C. No interactions between community type and selection line were detected. A combined model of line and bacterial types explained infection severity better than either factor alone.

Discussion

4

Our selection lines demonstrated striking heritability of parasite infection resistance and tolerance. Despite no significant differences in microbiome between the three lines, we found that bacterial community type, found by DMM machine‐learning, predicted infection severity. Host genetics and the host‐associated microbiome thus appear to independently predict infection severity, with similar effect sizes. We also observed notable differences between males and females in skin microbiome and in response to selection. We discuss each of these results in turn.

We found that both guppy resistance and susceptibility to G. turnbulli responded readily to a single generation of truncation selection. We found no notable differences between the three selection lines in our genetic metrics, indicating that the differences we found were not due to differences in inbreeding among them. Our narrow‐sense heritability estimates for resistance and tolerance were high, particularly for our susceptible lines. However, our values are comparable to heritability estimates in previous studies (Weiler al. 2025) and to the heritability of male guppy ornaments (J. Morris et al. 2020), which are correlated with G. turnbulli resistance (Stephenson et al. 2020).

Previous literature has demonstrated that tolerance is an essential factor to consider in guppy response to Gyrodactylus parasites (Jog et al. 2022; Stephenson et al. 2015; Tadiri et al. 2021). However, whether tolerance of infection is heritable has rarely been explored in this system. More broadly, there is now good evidence that tolerance has a heritable genetic basis across systems (Kutzer et al. 2023). Parasite resistance and susceptibility likely responded to selection through different mechanisms. The response to selection for susceptibility could have involved either an increase in tolerance, a loss of resistance, or some combination. Indeed, our experimental design may have favored tolerance in the susceptible line, as only individuals surviving severe infections would have contributed to future generations. Tolerance and resistance likely have different genetic bases, and hosts appear to invest in one strategy at the expense of the other (Råberg et al. 2007). This observation is supported by the fact that our resistant line fish had lower tolerance than control and susceptible line fish—at least among females—and that the two defense strategies were negatively correlated at the level of individual fish in our experiment. However, these tradeoffs do not always occur, and both tolerance and resistance traits can coexist within a single individual host and a population (Hardy et al. 2024; Pagán and García‐Arenal 2018).

Our observed high heritability of susceptibility (and potentially tolerance) has implications for our understanding of how guppy defenses might evolve in natural conditions (Walsman et al. 2022, 2024). High heritability suggests natural selection can act rapidly on host defense traits in populations under high parasite pressure, and indeed host defense appears to have the capacity to evolve rapidly in this system (Dargent et al. 2016; Karl P. Phillips et al. 2018). In some populations, infection prevalence is very high (> 90%; Clark et al. 2023; Stephenson et al. 2015), potentially selecting for increased infection tolerance and decreased resistance, a phenomenon known as the “resistance is futile” effect (Walsman et al. 2023). Importantly, here we deliberately excluded any potential role of parental effects of infection on offspring phenotypes, as is the case for almost all previous studies of host defense traits in this system. However, recent findings suggest that transgenerational immune priming may be important in this system, which would reduce the strength of selection on these traits, particularly when prevalence is high (Weiler et al. 2025).

We found that bacterial community type predicted infection severity (but interestingly not tolerance) independently, and to a similar degree as host genetics. Despite the fact that host characteristics such as the immune system, metabolism, and growth rate have been shown to interact with the microbiome (Speer 2022; Zheng et al. 2020), we found no difference between our lines in skin microbiome composition, alpha diversity, or bacterial community types. Potentially this is because fish skin microbiomes are highly variable and may not reflect the host's genetics but instead drift stochastically (Jones et al. 2022; Wang et al. 2023). Additionally, while we implemented several key controls to account for environmental microbes, that is, washing fish in sterile water, there could still be some transient environmental microbes within our samples. However, many of the top contributing taxa to the DMM bacterial community types (Flavobacterium, Limnobacter, Aeromonas, Acinetobacter, etc.) are commonly found on fish skin and mucus (Takeuchi et al. 2021; Hu et al. 2021; Bruno et al. 2023). More fundamentally, given that this ectoparasite lives at the host‐water interface and is directly exposed to the aquatic microbiome, environmental microbes present on host skin could represent biologically relevant interactions between parasite and water microbes.

Recently it has come to light that different organs may provide unique insights into host health, that is, eye microbiome predicted disease outcome but gut microbiome did not in Island foxes (DeCandia et al. 2025). Therefore, it is possible that microbiomes associated with other tissue could have been modified due to our selection, though we did not test for this. In contrast to our study, artificial selection for body size in gilthead sea bream ( Sparus aurata ) significantly altered the host‐associated microbiome, and the microbiome in the fast‐growing line had greater parasite colonization resistance (Piazzon et al. 2020). Similarly, size selection experiments indirectly affect guppy resistance to G. turnbulli (Bartuseviciute et al. 2022), potentially mediated through the microbiome. These results suggest that selection on host growth does influence the associated microbiome composition, potentially through effects on host metabolism, with downstream effects on parasite resistance (Bartuseviciute et al. 2022; Piazzon et al. 2020), whereas selection on parasite resistance appears not to have direct effects on skin‐associated microbial communities (this study).

Our findings thus contribute to the growing literature linking host‐associated microbiomes and infection severity, pointing to important interactions between the skin microbiome and parasites (Scheifler et al. 2022, 2023; Speer 2022; Trevelline et al. 2020). The bacterial taxon showing the biggest difference in relative abundance between the bacterial community types was the genus Limnobacter, with a lower relative abundance in bacterial community type A. Limnobacter could be of particular interest as it is linked to mucus secretions in fish skin (Hu et al. 2021). Gyrodactlids feed on epidermal cells and mucus (Bakke et al. 2007), therefore, the bacteria could be aiding in parasite colonization by providing resources. Limnobacter is also known for its ability to oxidize thiosulfate, potentially contributing to the chemical cues of infected guppies that uninfected fish avoid (Stephenson et al. 2018). More broadly, the microbiome's role in determining resistance has been seen in several taxa, where pre‐infection microbiome composition, alpha diversity, or microbiome function is a powerful tool to predict individuals that are most susceptible to parasites within a population (Bates et al. 2022; Lee et al. 2019). Intriguingly, we found five of the eight ASVs reported to be indicator taxa of epizootic dynamics in a recent study of skin microbial metabolism function in amphibian chytridiomycosis (Bates et al. 2022): Burkholderiales, Pseudomonadales, Sphingomonadales, Rhizobiales, Flavobacteriales (see Figure S2 for contributing species). While beyond the scope of our study, it has been shown that functional profiles of the microbiome likely play a critical role in the microbiome's ability to defend against parasitic infection (Bates et al. 2022).

Sex differences are pervasive throughout our results. Sexual dimorphism in immunocompetence has been reported in many taxa, including Trinidadian guppies. Our results reveal additional evidence of sex‐specificity in defenses against the specialist parasite G. turnbulli: resistance responded more strongly to selection among female than male guppies, as has been found previously (Dargent et al. 2016). Potentially, this result reflects the fact that males in nature and in our experiment were subject to sexual selection in addition to direct selection for parasite resistance. Females in many populations prefer males with larger orange ornaments (Endler et al. 2001; Houde 1988, 1997), with implications for disease dynamics (Rovenolt et al. 2024). Larger orange ornaments indicate male resistance to this parasite (Stephenson et al. 2020). Compared to females, males tend to develop lower infection loads in general in the lab (Johnson et al. 2011), are less likely, and less heavily infected in surveys of natural populations (Clark et al. 2023; Gotanda et al. 2013; Stephenson et al. 2015). Thus, if female preference has already eroded genetic variation in parasite resistance among males, and females across our lines selected the most resistant males as mates, this may explain the overall high and invariable resistance we observe in males. By contrast, and again in line with previous results from field surveys (Stephenson et al. 2015), females—particularly large ones—appear to be more tolerant of infection than males. Notably, we found evidence for a tradeoff between resistance and tolerance at the line level among females, but not among males, perhaps indicating different trade‐offs for each sex, though there was no sex difference in the tradeoff at the level of individual fish.

The strong sex‐specificity in host defense traits we observe here highlights the need to investigate aspects of guppy response to Gyrodactylus turnbulli beyond the MHC. Specifically, if the primary mechanism of resistance were driven by MHC parasite resistance, we would not see such pronounced sex differences in response to selection: MHC resistance is a mixture of both mother and father MHC genes. While our data indicate that the microbiome is not driving this sex difference in response to selection, the microbiome clearly contributes to parasite resistance in this system and appears to differ between the sexes: they differed in alpha diversity metrics of bacterial communities, as reported before (Kramp et al. 2022). This finding is similar to the sex differences in the skin microbiome in other fish, octopus and mammals (Coetzer et al. 2021; Rodríguez‐Barreto et al. 2023; Park et al. 2022). Recently, it was proposed that the sex difference in the microbiome of the lungs could be driving the severity difference to respiratory illness in humans (Beauruelle et al. 2021). These possibilities will hopefully inspire exciting and novel exploration of the evolutionary dynamics of the microbial communities' contributions to resistance and tolerance between the sexes.

Author Contributions

Rachael D. Kramp: conceptualization (equal), data curation (supporting), formal analysis (lead), funding acquisition (supporting), validation (lead), visualization (lead), writing – original draft (lead), writing – review and editing (lead). Mary J. Janecka: data curation (supporting), formal analysis (supporting), writing – original draft (supporting), writing – review and editing (supporting). Nadine Tardent: data curation (supporting), writing – review and editing (supporting). Jukka Jokela: funding acquisition (supporting), project administration (supporting), writing – review and editing (supporting). Kevin D. Kohl: formal analysis (supporting), writing – review and editing (supporting). Jessica F. Stephenson: conceptualization (lead), data curation (lead), formal analysis (supporting), funding acquisition (lead), investigation (lead), methodology (lead), project administration (lead), supervision (lead), writing – original draft (supporting), writing – review and editing (equal).

Funding

This work was supported by the Fisheries Society of the British Isles (studentship to J.F.S.), the Center for Adaptation to a Changing Environment at ETH Zürich (fellowship to J.F.S.), the University of Pittsburgh, and Howard Hughes Medical Institute (Gilliam Fellowship for Advanced Study to R.D.K.).

Ethics Statement

The permit to collect and export guppies was granted by the Director of Fisheries in the Ministry of Agriculture, Land and Fisheries Division, Aquaculture Unit of the Republic of Trinidad and Tobago. Import to Cardiff University was under Cefas APB authorisation number CW054‐D—187A. Fish husbandry and experimental infections of the F1 generation at Cardiff were under UK Home Office license PPL 302876, with approval by the Cardiff University Animal Ethics Committee. Guppies were shipped to Zürich, and their husbandry and experimental infections of the F3–6 generations were approved by the Veterinäramt of Kanton Zürich (ZH177/15).

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Appendix S1: ece372923‐sup‐0001‐AppendixS1.pdf.

Appendix S2: ece372923‐sup‐0002‐AppendixS2.pdf.