Local parasite pressures and host genotype modulate epigenetic diversity in a mixed‐mating fish

Abstract Parasite‐mediated selection is one of the main drivers of genetic variation in natural populations. The persistence of long‐term self‐fertilization, however, challenges the notion that low genetic variation and inbreeding compromise the host's ability to respond to pathogens. DNA methylation represents a potential mechanism for generating additional adaptive variation under low genetic diversity. We compared genetic diversity (microsatellites and AFLPs), variation in DNA methylation (MS‐AFLPs), and parasite loads in three populations of Kryptolebias hermaphroditus, a predomintanly self‐fertilizing fish, to analyze the potential adaptive value of DNA methylation in relation to genetic diversity and parasite loads. We found strong genetic population structuring, as well as differences in parasite loads and methylation levels among sampling sites and selfing lineages. Globally, the interaction between parasites and inbreeding with selfing lineages influenced DNA methylation, but parasites seemed more important in determining methylation levels at the local scale.


| 8737
BERBEL-FILHO Et aL. a wide range of environments (Zhang, Zhang, & Barrett, 2010), suggesting that their adaptation and long-term survival could be facilitated by other factors in addition to genetic variability (Verhoeven & Preite, 2014).
Negative associations between genetic diversity and parasite loads have been previously observed in mixed-mating animals (Ellison, Cable, & Consuegra, 2011;Lively & Morran, 2014), with inbred individuals usually harboring more parasites. The relationship between epigenetic variation, parasites, and mixed-mating, however, has not been explored.
Here, we compared genetic diversity, variation in DNA methylation, and parasite loads in three natural populations of the mixed-mating mangrove killifish Kryptolebias hermaphroditus distributed along the Brazilian coast (Tatarenkov et al., 2017). The genus Kryptolebias contains the only known mixed-mating vertebrate species (K. marmoratus and K. hermaphroditus), characterized by variable rates of selfing and outcrossing (Tatarenkov et al., 2017). Populations of both species consist mainly of self-fertilizing hermaphrodites and varying levels of males at low frequencies (Berbel-Filho, Espirito-Santo, & Lima, 2016;Tatarenkov et al., 2017), and exhibit high levels of homozygosity (Tatarenkov et al., 2017;Tatarenkov, Lima, Taylor, & Avise, 2009), suggesting that self-fertilization is the most common mode of reproduction .
We analyzed microsatellites (previously shown to correlate with parasite loads in the closely related K. marmoratus, see Ellison et al., 2011) and genome-wide methylation based on identification of anonymous CpG by methylation-sensitive AFLP (MS-AFLPs, previously used in nonmodel organisms) to identify epigenetic variation associated with parasite loads (Wenzel & Piertney, 2014

| Study system, field sampling, and parasite screening
A total of 128 specimens of K. hermaphroditus were collected using hand-nets from three sampling sites on isolated mangroves on the northeastern coast of Brazil between January and September 2015: Ceará-Mirim River-Site 1; Curimataú River-Site 2; Ipojuca River-Site 3 (Figure 1). K. hermaphroditus is distributed along the Brazilian coast (Tatarenkov et al., 2017) and is typically found in shallow pools of high salinity levels (>30 ppt), clear waters, and muddy substrates, where there are few other sympatric fish (Berbel-Filho et al., 2016;Lira, Paiva, Ramos, & Lima, 2015). All specimens displayed the common hermaphrodite phenotype (dark color with well-defined ocellus on the caudal fin; Costa, 2011). Fish were euthanized using an overdose of tricaine methanesulfonate (MS-222) following UK Home Office Schedule 1 (Hollands, 1986), standard length was measured using a digital calliper (mm), and the whole fish were preserved in 95% ethanol at −20°C for parasite screening and DNA extraction.
In the laboratory, fish were dissected and screened for both external and internal parasite infections using a dissecting microscope following the methods of Ellison et al. (2011). Macroscopic parasite analyses focused on the three most common types of parasites identified. To assess the reliability of parasite screening, a subsample of five fish was examined by a different observer and the agreement was 100%. We defined parasite loads using a scaled measure of parasite abundance, where for each parasite morphotype (i), the number of parasites per individual (Ni) was divided by the maximum number found across all individuals (Nimax). The final value of the scaled parasite load represents the sum of scaled parasite loads across all parasite types. Given their uneven abundance (Table 1), this approach minimizes bias when parasite loads are heavily influenced by a very abundant parasite type (in our case, bacterial cysts) (Bolnick & Stutz, 2017).
The inbreeding coefficient (F IS ) was calculated in GENEPOP. Global heterozygosity for individual fish was estimated using the homozygosity by locus index (HL) implemented in the Excel macro Cernicalin v 1.3 (Aparicio, Ortego, & Cordero, 2006).
We also used the Bayesian clustering algorithm INSTRUCT (Gao, Williamson, & Bustamante, 2007) to estimate the optimal number of selfing lineages (k) with four simultaneous chains of 2,000,000 MCMC runs, 10 as thinning, and 100,000 of burn-in period, resulting in 100,000 interactions for each chain. The potential number of k tested ranged from 2 to 12. We used the individual q-values (the likelihood of membership to a particular genetic cluster or selfing lineage) from INSTRUCT to classify individuals as either selfed or outcrossed (Vähä & Primmer, 2006). A threshold of q-value ≥ 0.9 was used to classify selfed individuals, while <0.9 represented hybrids between two different selfing lineages, suggesting an outcrossing event (Ellison et al., 2011;Vähä & Primmer, 2006). Pairwise F ST values among sampling sites and selfing lineages were estimated with Arlequin v. 3.5.2.2 (Excoffier & Lischer, 2010) using 10,000 permutations. We used hierarchical analysis of molecular variance (AMOVA) to investigate population structuring among sampling sites and selfing lineages (according to individual q-values) using 10,000 randomizations. Differences between selfed and outcrossed groups in the total number of parasites and homozygosity by locus (microsatellites) were analyzed using median Mann-Whitney rank tests in R v. 3.3.
TA B L E 1 Genetic diversity (at 27 microsatellite loci), mean parasite number (standard deviation in brackets), and parasite prevalence in Kryptolebias hermaphroditus at sampling sites in northeastern Brazil

| Epigenetic analysis
We used methylation-sensitive amplified fragment length polymorphisms (MS-AFLPs) to assess genome-wide DNA methylation patterns . Epigenetic (MSL) and genetic (NML) differentiation at AFLPs among sampling sites, selfing lineages, and between outcrossed and selfed groups was assessed by AMOVA with 10,000 randomizations.
Epigenetic (MSL) and genetic (AFLP and microsatellites) differentiation among sampling sites, selfing lineages, and inbreeding status was visualized by principal coordinates analysis (PCoA). Mantel tests based on distance matrices (Mantel, 1967) were used to test for potential correlations between epigenetic and genetic data for MSL, NML, and microsatellites using GENALEX v. 6.5 with and 10,000 permutations. To identify disproportionately differentiated methylation states, we used a F ST outlier approach implemented in BayeScan 2.1 (Foll & Gaggiotti, 2008;Perez-Figueroa et al., 2010), with 2 × 10 6 iterations (thinning interval 20 after 20 pilot runs of 10 4 iterations each) and a burn-in of 5 × 10 5 . We tested for outliers based on the MSL data generated on the comparisons among sampling sites, selfing lineages, and between inbreeding status (inbred or outbred).

| Statistical analyses
A Kruskal-Wallis test was used to examine the differences on Model selection was conducted using the multimodel averaging approach implemented in the R package glmulti v 1.0.7 (Calcagno & de Mazancourt, 2010). We chose the minimal adequate models based on the lowest AICc values (Akaike information criterion corrected for small sample size), Akaike weight (W i ), and evidence ratios (Burnham & Anderson, 2004). Models (within 2 AIC units) were also reported. Predictors were checked for collinearity using pair.panels function in R package psych (Revelle et al., 2019). Model residuals were checked and assumptions validated.
To disentangle potential confounding effects arising from the unequal distribution of selfing lineages among sampling sites (i.e., five lineages are exclusive to a particular sampling site, Table S1), we repeated the analyses (AMOVA, Mantel test, PCoA, and GLMs) for both genetic (microsatellites and AFLPs) and epigenetic (MSL) data using only individuals from Site 1 (68 individuals for microsatellites and 62 for MS-AFLPs), as this was the only site with more than two selfing lineages (Table S1).

| Parasite screening
Bacterial cysts were present on the gills and consisted of white to yellow spherical cysts circumscribed by a capsule, which resulted in hypertrophied gill filaments. They were the most common type of pathogen appearing in 83.6% of the individuals screened, with a prevalence ranging from 1 to 19 (mean = 2.73, SD = 2.99), and were more prevalent in Site 1 (mean = 3.16, SD = 3.16), followed by Site 2 (mean = 2.66, SD = 3.10) and Site 3 (mean = 1.27, SD = 0.80). The second most common macroscopic parasites were protozoan cysts, which consisted of small dark oval cysts over the gill arch and filaments. In total, 19.53% of the total number of individuals were infected with these cysts, ranging from 1 to 6 (mean = 0.54, SD = 1.26).

| Genetic diversity and population structuring based on microsatellites
No linkage disequilibrium was detected between any pair of microsatellite loci. As expected from the high levels of self-fertilization of the species, no loci were found to be in Hardy-Weinberg equilibrium, and all 27 microsatellite loci showed an excess of homozygotes. The global homozygosity index (HL) was very high (mean = 0.95), as well the estimated selfing rates (Table 1) Site 3) (Figures 1 and 2; Table S1). High F ST values were found both among sampling sites (mean = 0.28, SD = 0.02) and selfing lineages (mean = 0.32, SD = 0.05). All pairwise comparisons were highly significant (Table S2).  (Table 2).
These patterns were also seen on PCoA, with individuals generally clustering by selfing lineages in the microsatellites data (25.84% of overall variation), with individuals from lineage 4 being the most differentiated from the other lineages on Site 1. In this site, substantial overlap was found among selfing lineages and between selfed and outcrossed, despite its significant differences (F ST = 0.03, p = 0.001) (Table S4; Figure S3).

| Genetic and epigenetic variability and population structuring based on MS-AFLPs
The epigenetic analysis identified 381 MS-AFLP loci, of which 267 (70.07%) were methylation-susceptible loci (MSL) and 106 F I G U R E 2 Genetic assignment of Kryptolebias hermaphroditus to six selfing lineages using INSTRUCT. Each individual is represented by a bar, which represents the likelihood of the individual to belong to a specific genetic cluster (color) which is within the normal reproducibility range for AFLPs genotyping (Bonin et al., 2004). AMOVA for reproducibility also revealed no significant differences between methylation and AFLP variation patterns between original and replicated set of individuals (Table S3). Average methylation ranged from 47.51% on lineage 2 to 38.17% on lineage 5, and was 44.82% for inbred and 45.77% for outbred individuals.  (Table 3). As with microsatellites, no clear genetic (at AFLPs) or epigenetic differentiation was found between selfed and outcrossed individuals ( Figure S2).
No MSL epiloci were identified as an F ST outlier in any of the comparisons.  (Table S4). In the PCoA, substantial overlap was found among selfing lineages and between selfed and outcrossed individuals ( Figure S3). Mantel tests between genetic and epigenetic data indicated a significant positive association between AFLPs and MSL data (r = 0.21; p < 0.001), but not between microsatellites and MSL (r = −0.005; p = 0.45).

| Parasite loads, genetic variation, and epigenetic variation
According to a multimodel testing approach, the most plausible  (Tables S5-S6).
Overall, the results of the single-taxa models (using number bacterial cysts) were very similar to those for scaled parasite loads.
The best model to explain the proportion of methylated loci included selfing lineage, and the interactions between selfing lineage and bacterial cysts, and selfing lineage and inbreeding (Table S7).
When using only individuals from Site 1 (to remove any potential confounding effect between sampling site and selfing lineages) for the proportion of methylated loci, the model with the lowest AIC indicated that selfing lineage, inbreeding, and the interactions between inbreeding and selfing lineage and inbreeding and scaled parasite loads were all significant predictors (Table S8). However, the second best-fitting model (ΔAICc = 0.02) explained the same amount of variation (weight = 0.39) and the evidence ratio (−0.66) suggested that it was more likely (evidence ratio of 1.50) than the first model. This second model indicated that overall, the proportion of methylated DNA significantly increased with scaled parasite loads (estimate = 0.43, SE = 0.11, p = 0.03) and that DNA methylation levels were also affected by the interaction between scaled parasite loads and inbreeding (estimate = −1.29, SE = 0.38, p < 0.001), with inbred individuals having increased methylation levels with increased parasite loads, while outbred individuals had decreased methylation levels with increased parasite loads (Figure 3d; Table 4).

| D ISCUSS I ON
DNA methylation could play an important adaptive role in organisms with low heterozygosity, including self-fertilizing species, potentially increasing their plasticity capacity to cope with environmental change (Douhovnikoff & Dodd, 2015;Verhoeven & Pretie, 2014). However, our results did not indicate significant differences in genome-wide DNA methylation variation between selfed and outcrossed individuals, and our models only identified inbreeding status (defined as originating from selfing or outcrossing) significantly related to DNA methylation via its interaction with selfing lineage (all sampling sites) and parasites (at the local scale in Site 1). Higher variation in DNA methylation has been reported for clonal and inbred individuals Massicotte & Angers, 2012;Nakamura & Hosaka, 2010;Richards et al., 2012;Veerger et al., 2012) and has been interpreted as an adaptive mechanism to compensate for low genetic variation , or as a potential consequence of inbreeding (as in Vergeer, Wagemaker, & Ouborg, 2012) responsible, at least in part, for inbreeding depression (Nakamura & Hosaka, 2010;Vergeer et al., 2012). Yet, our results suggest that, at least in this species, either inbreeding does not affect genome-wide DNA methylation variation or it does in a gene-specific manner (Venney, Johansson, & Heath, 2016), although further research would be needed to address this question.
We found that the different selfing lineages of Kryptolebias hermaphroditus distributed in three sampling sites of northeastern Brazil differed significantly in parasite loads, genetic composition, and DNA methylation patterns, which might indicate specific interactions between host genotypes, epigenotypes, and parasites (Dybdhal & Lively, 1998;Ebert, 2008). Previous studies on mangrove killifishes had identified extensive genetic structuring both between (Tatarenkov et al., , 2017 and within mangrove systems even at close geographical proximity (Ellison et al., 2012;Tatarenkov, Earley, Taylor, & Avise, 2012;Tatarenkov et al., 2007), as a consequence of the self-fertilizing nature of these fish. We found strong evidence of genetic structuring between sampling sites and selfing lineages using microsatellites, but lower differentiation for AFLP genetic markers (likely due to the different mutation rate of the markers) and epigenetic markers (MS-AFLPs).
Overall, more inbred individuals harbored higher parasite loads than their outcrossed counterparts, supporting the prediction that low heterozygosity due to self-fertilization may reduce fitness (considering parasite loads as a proxy for pathogen pressure), as for other mixed-mating species (Ellison et al., 2011;King, Jokela, & Lively, 2011;Lively & Morran, 2014). Extensive periods of self-fertilization can reduce offspring fitness due to the accumulation of deleterious alleles and inbreeding depression (Charlesworth et al. 1993 Schmidt, Gelarden, Parrish, & Lively, 2011), which can generate genetic diversity to face natural enemies, such as parasites ). Here, the relationship between parasites and inbreeding status (selfed or outcrossed) suggests that outcrossing might confer a fitness advantage (in terms of parasite loads), even when it occurs at very low frequencies (Ellison et al., 2011). However, despite the adaptive potential of outcrossing, the main reproductive mode of K. hermaphroditus seems to be self-fertilization (Tatarenkov et al., 2017). This suggests that other evolutionary mechanisms may be balancing the harmful effects of parasite infections or that parasite selection is of low , as theory predicts that low selection levels imposed by natural enemies are consistent with the maintenance of asexual reproduction (Judson, 1997;Ladle, Johnstone, & Judson, 1993). For example, in the mixed-mating Potamopyrgus snails, the oldest asexual lineages are restricted to populations where parasites are rare (Neiman, Jokela, & Lively, 2005 (Teh et al., 2014). Using data from all sampling sites, we found that genome-wide DNA methylation was strongly influenced by selfing lineage, and only at a smaller scale by inbreeding through its interaction with selfing lineage (Bell et al., 2011;Bjornsson et al., 2008;Dubin et al., 2015;Gertz et al., 2011). Strong epigenetic differences between selfing lines had been identified previously in K. marmoratus Ellison et al., 2015), indicating an important role of the genetic background in the epigenetic variation of mangrove killifishes. In addition, we also found a significant correlation between DNA methylation and genetic variation (at both AFLP and microsatellites data), suggesting that autonomous variation in DNA methylation may be limited (Dubin et al., 2015).
Several abiotic and biotic factors, including parasites (Hu et al., 2018;Norouzitallab et al., 2014) as well as stochastic epimutations (Leung et al., 2016), are known to influence DNA methylation variation. Our results showed that genome-wide DNA methylation levels for all sampling sites were significantly influenced by parasite loads through the interaction with selfing lineage, suggesting a potential genotype-by-environment interaction on parasite responses. Yet, as most of the selfing lineages were exclusive to specific sampling sites, we could completely discard confounding effects between both variables. In fact, selfing lineage did not affect genome-wide DNA methylation levels in Site 1, but only parasites and their interaction with inbreeding status. Increasing evidence has been showing that DNA methylation is involved in the modulation of host-pathogen interactions (Gómez-Díaz et al., 2012). The bacterial parasite Wolbachia, for example, alters host genome-wide DNA methylation patterns resulting in the feminization of infected leafhoppers (Zyginidia pullulan) to increase its transmission (Negri et al., 2009). Experiments in plants with both hyper-and hypomethylated mutants indicate that genome-wide DNA demethylation enhances immune responses to both bacterial (Dowen et al., 2012;Yu et al., 2013) and fungal (Le et al., 2014) infections. Additionally, almost half of the resistance genes in the The relationship between parasite loads and outcrossing seems to be common to several mixed-mating species (Ellison et al., 2011;King et al., 2011;Steets, Wolf, Auld, & Ashman, 2007) in addition to K. hermaphroditus, suggesting that the influence of parasites in the regulation of mixed-mating could be generalized. The extent of this relationship, however, may depend on the severity of the selection imposed by coevolving parasites (Lively & Morran, 2014).
Our results indicate that genotype composition (and its interaction with inbreeding) may be important in DNA methylation responses to environmental variation in wild populations, and that, if DNA methylation responded in a genotypic-specific manner to parasites pressures, it could contribute to local adaptation (Foust et al., 2016;Smith, Mártin, Nguyen, & Mendelson, 2016). The mangrove killifish, with its naturally inbred populations and marked methylation differences between populations and genotypes, represents an ideal model to analyze the relative roles of genetic and epigenetic diversity in modulating local adaptation.