A globally ubiquitous symbiont can drive experimental host evolution

Organisms harbour myriad microbes which can be parasitic or protective against harm. The costs and benefits resulting from these symbiotic relationships can be context‐dependent, but the evolutionary consequences to hosts of these transitions remain unclear. Here, we mapped the Leucobacter genus across 13,715 microbiome samples (163 studies) to reveal a global distribution as a free‐living microbe or a symbiont of animals and plants. We showed that across geographically distant locations (South Africa, France, Cape Verde), Leucobacter isolates vary substantially in their virulence to an associated animal host, Caenorhabditis nematodes. We further found that multiple Leucobacter sequence variants co‐occur in wild Caenorhabditis spp. which combined with natural variation in virulence provides real‐world potential for Leucobacter community composition to influence host fitness. We examined this by competing C. elegans genotypes that differed in susceptibility to different Leucobacter species in an evolution experiment. One Leucobacter species was found to be host‐protective against another, virulent parasitic species. We tested the impact of host genetic background and Leucobacter community composition on patterns of host‐based defence evolution. We found host genotypes conferring defence against the parasitic species were maintained during infection. However, when hosts were protected during coinfection, host‐based defences were nearly lost from the population. Overall, our results provide insight into the role of community context in shaping host evolution during symbioses.

its host and its positioning along the mutualist-parasite continuum can change due to host genetics (Tianero et al., 2015;Wollenberg & Ruby, 2012), environment (Baker et al., 2018;Holt et al., 2020;Khojandi et al., 2019) or symbiont conflict/cooperation in the microbiota (Fraune et al., 2015). Defensive symbioses, whereby a symbiont protects its host from enemy attack, can be particularly dynamic along the continuum. Parasites might shift to being protectors, or vice versa depending upon the symbiont community (King et al., 2016), whilst the effectiveness of protection can vary depending upon host genetic background (Rafaluk-Mohr et al., 2018;Vorburger & Gouskov, 2011).
Some beneficial microbial associations can lead to the evolution of host dependence on the relationship (Fisher et al., 2017).
However, the extent to which context-dependent transitions along the parasite-mutualist continuum can alter patterns of host evolution in defensive symbioses are unclear. Whilst parasites should favour the spread of host genotypes defended against infection, this evolutionary outcome could be hampered by a coinfecting protective symbiont. Symbiont-mediated protection could cause hosts to plastically divest in costly host defences (Ford & King, 2021) or reduce their selective advantage (Martinez et al., 2016) leading to host dependence (Metcalf & Koskella, 2019). Alternatively, if protection is brought about by a joint response between the host immune system and symbiosis (Johnston & Rolff, 2015), then host-encoded defences in the population may be maintained (Hrček et al., 2018).

Moreover, interactions between symbiont species and different
host genetic backgrounds have the potential to alter host fitness landscapes (Clay et al., 2005;Jaenike et al., 2010;Parker et al., 2017;Weldon et al., 2019). By generating novel asymmetries in host competition, symbionts could maintain host genotypes in the population that would otherwise be lost in the absence of symbiosis (Murfin et al., 2019;Park, 1948). Understanding these evolutionary interactions in the ecological noise of wild symbioses remains a challenge given the heterogeneity in both host genetics and microbial symbiont communities.
Here, we use Caenorhabditis nematodes and their natural bacterial symbionts of the Leucobacter genus to examine phenotypic variation in symbioses and the evolutionary consequences for hosts. Leucobacter spp. have been shown to colonise Caenorhabditis spp. in Cape Verde, South Africa, Japan, France, India, the USA and have been identified from metagenomic sequences (Clark & Hodgkin, 2015;Hodgkin et al., 2013;Johnke et al., 2020;Muir & Tan, 2008;Percudani, 2013;Samuel et al., 2016;Somvanshi et al., 2007). Leucobacter species can be parasitic or protect hosts from parasitic infection (Hodgkin et al., 2013). We first sought to establish the global distribution and host range of Leucobacter to better understand whether this genus is predominantly specialised to colonising nematodes, or displays a more generalist habit. Using data from the publicly available Earth Microbiome Project (EMP) database (Thompson et al., 2017), we reveal that Leucobacter is globally distributed and often found in symbiosis with animals and plants.
We then assessed the range of Leucobacter interactions with wild C. elegans (our focal host taxon) isolates that span our mapped symbiont distribution. Next, we analysed 16S amplicon sequencing data to show that coinfection of wild Caenorhabditis nematodes by multiple Leucobacter sequence variants is common. Based on our findings of heterogeneity in virulence of Leucobacter isolates and the occurrence of coinfection in nature, we hypothesised that Leucobacter community composition may have implications for host fitness and evolution. We used experimental evolution to test how the ability of Leucobacter species to be both parasitic and protective can shape the speed and trajectory of C. elegans genetic composition and defence evolution. Finally, we assessed how host genotype can impact the epidemiology of Leucobacter in terms of certain genotypes super-spreading to more susceptible hosts. We discuss how such ecological processes may be linked to observed host evolutionary trends.

| Leucobacter global distribution
The EMP tool Redbiom (McDonald et al., 2019) was used to search for the presence of operational taxonomic units (OTUs) assigned to the Leucobacter genus in the Qiita database (accessed 06.05.2020). The commands "redbiom search taxon" and "redbiom search features" were used to search 16S rRNA libraries from the EMP sequencing context "Pick_closed-reference_OTUs-Greengenes-Illumina-16S-V4-150nt". Samples with assigned EMP ontology (EMPO) metadata and geographic coordinates were selected for downstream analysis.
Overall, Leucobacter OTUs were assigned to 13,992 EMP samples of which 13,715 samples met criteria (sufficient metadata) for downstream analysis. Samples that did not have EMPO annotation or had a biologically irrelevant annotation were manually annotated. A map of the distribution of sampling locations of EMP samples was plotted using R v 3.6.0 (http://www.r-proje ct.org/).

Previous microbiome studies have found Leucobacter in both wild
Caenorhabditis and their natural environment (Johnke et al., 2020;Samuel et al., 2016); however, Leucobacter prevalence and colonisation dynamics are poorly characterised. To better assess the prevalence and occurrence of coinfection of Leucobacter in wild Caenorhabditis, we analysed 16S amplicon sequencing data from a large published study that collected microbiome samples from Caenorhabditis populations in France, Germany and Portugal (Dirksen et al., 2016). Sequence data was analysed using the R package dada2 (Callahan et al., 2016). Primers were removed and forward reads were trimmed to 240 bases and reverse reads were trimmed to 220 bases. Reads were quality scored and trimmed at the first appearance of a base with quality score of two or lower. Reads containing nonassigned bases (N) and bases with an expected error rate higher than two were removed. Reads matching the PhiX sequencing standard genome were also removed. DADA2 partitions sequencing error in each read from true genetic diversity and categorises reads together into "amplicon sequence variants" (ASVs). ASVs in our data set were computed and paired end reads were merged into single consensus reads. Chimeric sequences were removed from the data set. Taxonomy was assigned to ASVs using the naïve Bayesian classified algorithm (Wang et al., 2007) with the Silva training set (version 138) (Quast et al., 2012). Downstream analysis was conducted using the R package phyloseq (McMurdie & Holmes, 2013). We calculated the sample prevalence of Leuobacter by counting the occurrence of at least one Leucobacter ASV in a single host. We further examined the potential for coinfection by Leucobacter genotypes or species by computing the number of individual hosts colonised by more than one Leucobacter ASV. Finally, we summarised the occurrence of Leucobacter ASVs by nematode species to better characterise niche breadth and colonisation dynamics in wild nematodes.

| Bacteria and nematode maintenance
We used Leucobacter musarum sp. nov. subsp. musarum subsp. nov strain CBX152 T (herein Verde2) and Leucobacter celer subsp. astrifaciens subsp. nov strain CBX151 T (herein Verde1) which were isolated from wild Caenorhabditis found in rotting banana stem on the island of San Antao, Cape Verde (Republic of Cape Verde) (Hodgkin et al., 2013). Under standard culture conditions, Verde1 colonises the nematode cuticle where it may cause slight impairment to host growth, but does not prevent successful reproduction. Meanwhile Verde2 is highly virulent, invading the rectal region of nematodes where it goes onto induce vacuole formation in the body cavity and subsequent lysis of host internal structures (Hodgkin et al., 2013). In addition, we used Leucobacter isolates collected from Caenorhabditis sp. in South Africa (JUB18) and France (JUB111). Little is known of the effects of JUB18 and JUB111 on Caenorhabditis health. These isolates were included primarily to compare their virulence with those of Verde1 and Verde2 across C. elegans genotypes. Each Leucobacter isolate used in assays was grown from a single colony to stationary phase in Luria-Bertani (LB) broth (10 g tryptone, 5 g yeast extract, 5 g NaCl, 1 L H2O) at 30°C with constant shaking at 200 rpm (Hodgkin et al., 2013). Bacteria were then standardised to 10 7 colony forming units per ml to allow for doses to be maintained consistently within and across experiments.
We used the hermaphroditic, wild-type C. elegans strain N2 and the following wild caught isolates: ED3040 (South Africa), JU1240 (France), and JU1088 (Japan). For the evolution experiment, we also used the N2-mutant srf-2 (CB6850 = ccIs4251 srf-2(yj262)). This mutant has altered surface glycosylation rendering it defended against Verde2, but hyper-susceptible to Verde1. N2 worms show the opposite pattern of susceptibility (Hodgkin et al., 2013). The srf-2 mutant had a green fluorescent protein reporter enabling it to be distinguished from the unlabelled N2 genotype. The GFP insertion has previously been demonstrated to not carry host fitness costs (Hodgkin et al., 2001). Nematodes were sourced from the Caenorhabditis Genetics Centre (University of Minnesota, Minneapolis, MN) and the laboratories of Marie-Anne Felix and Jonathan Hodgkin.
Nematodes were maintained on nematode growth medium (NGM) (3 g NaCl, 2.5 g Bacto-Peptone, 17 g agar, 1 ml cholesterol, 975 ml H 2 O) seeded with Escherichia coli OP50 (Caenorhabditis Genetics Centre, University of Minnesota, Minneapolis, MN). This food bacterium was grown at 30°C shaking at 200 rpm overnight in LB broth, with 100 µl of culture subsequently spread onto each NGM plate and incubated overnight at 30°C. Nematodes were treated with bleach (NaClO) and sodium hydroxide solution to synchronise for life-stage and collect sterile unhatched eggs (Stiernagle, 2006).

| Caenorhabditis elegans-Leucobacter interactions
We tested for variation in host susceptibility to Leucobacter and in virulence of Leucobacter isolates. We exposed three wild-caught C. elegans isolates originating from South Africa (ED3040), France (JU1240) and Japan (JU1088), with the laboratory type strain (N2), to four Leucobacter isolates collected from wild Caenorhabditis in Cape Verde (Verde1 and Verde2), South Africa (JUB18), and France (JUB111) ( Table S1). Host genotypes were also maintained on OP50 as a control group. Each host-symbiont combination was replicated four times. For symbiont exposures, microbial lawns were grown by streaking 74 µl of inoculum containing 20% Leucobacter and 80% OP50 onto 5.5 cm NGM plates and incubating at 25°C for 24 h.
Twenty C. elegans larval stage 4 (L4) nematodes were then picked onto their respective NGM symbiont exposure treatment and left for 24 h at 25°C. After 24 h, the number of live and dead C. elegans were counted for each plate. Nematodes were considered dead when they did not respond to touch with a platinum wire.

| Variation in host mortality with symbiont community
Exposure plates were made following the same procedure, with the exception that Leucobacter coinfection consisted of 10% Verde1, 10% Verde2 and 80% OP50. Four replicate populations were used for eight treatment groups (N2 and srf-2 each exposed to: OP50, OP50+Verde1, OP50+Verde2, OP50+Verde1+Verde2). Twenty L4 N2 or srf-2 nematodes were placed onto symbiont exposure plates and incubated at 25°C. After 24 h, the total number of live and dead nematodes were counted on each plate. The proportion of dead nematodes were calculated as a measure of host susceptibility.

| Experimental evolution
We hypothesised that variation in mortality in the host prereproductive period driven by the symbiont community would impose diversifying selection. To test this hypothesis, we conducted an evolution experiment in which populations of N2 and srf-2 L4 stage nematodes were exposed for 24 h to Verde1, Verde2, or both simultaneously on 9 cm NGM plates. For the single Leucobacter exposure, 160 µl of OP50 was mixed with 40 µl of Verde1 or Verde2 and spread on a 9 cm NGM plate for each replicate. For the OP50+Verde1+Verde2 treatment, 160 µl of OP50 with 10% (20 µl) Verde1 and 10% (20 µl) Verde2 was spread onto NGM plates. The control group consisted of only 200 μl of OP50.
Inoculated plates were incubated for 24 h at 25°C to grow bacterial lawns. The experiment was started with approximately 500 of the N2 and srf-2 genotype (n per plate = 1000) per replicate population.
Nematode-symbiont exposures were incubated for 24 h at 25°C.
After 24 h of symbiont exposure, unhatched sterile eggs were collected from nematode populations after bleaching. After 12 h incubation in M9 buffer, population genotype frequency estimates were determined by counting hatched L1 nematodes for each treatment in eight 2 µl drops of M9 buffer under a fluorescent microscope. Approximately 1000 nematodes from each replicate were then plated onto 200 µl OP50 plates and raised to L4 stage before being re-exposed to their respective symbiont treatments. This procedure was performed for each of the three host passages and was replicated five times.

| Contact-driven transmission
We hypothesized that the maintenance of N2 genotype during Verde2 exposure may be linked to host dilution effects (Rudolf & Antonovics, 2005). The high population frequency of srf-2 worms could limit frequency-dependent transmission of Verde2 in the susceptible N2 host genotype. We further predicted that the rapid extinction of srf-2 worms in the Verde1 treatment may be facilitated by increased parasite transmission ("super-spreading") from defended N2 worms. We tested these predictions by introducing Verde1 or Verde2 exposed N2 or srf-2 worms into naïve populations of both genotypes, thus minimising environmental symbiont to host transmission. Infected "source" populations of monoculture N2 and srf-2 C. elegans were generated by exposing L4 stage worms to Verde1, Verde2 or an E. coli OP50 control following the same procedure as for the "experimental evolution" (above). In addition, unexposed "sink" L4 populations of monoculture N2 and srf-2 worms were prepared on 5.5 cm NGM +OP50 plates with 100 worms per plate.
After 10 h Leucobacter exposure, five worms from source populations were picked onto sink population plates and the total mortality on each plate was counted over a 72 h period. This experiment was replicated three times.

| Statistical analysis
All statistical analyses were conducted in R version 3.6.0 (http:// www.r-proje ct.org/). To compare the proportion of host mortality among natural isolate interactions, we used a Kruskal Wallis test followed by a pairwise Wilcoxon rank sum test (function "pairwise.wilcox. test") with multiple comparison correction using the false discovery rate (FDR) procedure. We tested for differences in the proportion of host mortality between C. elegans genotypes exposed to different symbiont communities using a Kruskal Wallis test followed by a pairwise Wilcoxon rank sum test to test for differences in relative fitness between experimental groups with FDR correction of p-values.
For the evolution experiment, we calculated the relative fitness (ln w) (Hartl & Clark, 1997) of N2 compared to srf-2. We considered the initial frequency of each at 50% and the estimated frequencies of both genotypes at the end of the experiment. Relative fitness (ln w) was calculated as a function of the change in genotype frequency using the following equation: ln(w) = ln(f N2/f srf-2) where f represents genotype frequencies. In the few cases where one genotype had a frequency of zero, we assumed a frequency of 1/(n+1) for that genotype (Zbinden et al., 2008), where n is the total number of individuals genotyped in that replicate. We used a Kruskal-Wallis test followed by a pairwise Wilcoxon rank sum test to test for differences in relative fitness between experimental groups with FDR correction of p-values for multiple comparisons.
To compare the proportion of mortality in different treatment groups for the transmission experiment, we used the R package brglm (Kosmidis & Firth, 2021) with treatment group as a predictor and mortality fitted with a binomial distribution as the response variable. The overall significance of model predictors was assessed using a Wald Chi square test (Anova type=II command, car package).
Post hoc comparisons between treatment groups were performed using the multcomp package (Hothorn et al., 2008).   (3) = 1.52, p = .68, Figure 1b). Given our findings that Leucobacter isolates vary in virulence and that coinfection of Leucobacter ASVs occurs in wild nematodes, we hypothesised that Leucobacter community composition may impact host fitness. We first confirmed and then further quantified findings from a previous study to show that coinfection by the Cape Verde Leucobacter species in individual nematodes shifts the symbiosis along the parasite-mutualist continuum (Hodgkin et al., 2013). We singly and doubly infected L4 N2 nematodes -naturally defended against Verde1 and susceptible to Verde2 -as well as srf-2 mutants with altered surface antigenicity conferring the opposite susceptibility pattern. After 24 h, Verde2 was more virulent to N2s, and Verde1 to srf-2 nematodes (Kruskal-Wallis: χ 2 (7) = 27.83, p < .001, pairwise Wilcoxon test, p = .05, Figure 2a, Table S2). Coinfection in N2s resulted in symbiont-mediated protection, with a 31.4% average reduction in mortality compared to Verde2 single infection (pairwise Wilcoxon test, p = .05, Figure 2a, Table S2). This outcome was not replicated in srf-2 nematodes (pairwise Wilcoxon test, p > .05, Table   S2). This result suggests that strong genetic-based defence against a single parasite has the potential to limit broader aspects of immunity.

| RE SULTS
Using experimental evolution, we showed that after only three passages, genetic-based resistance in srf-2 to Verde2 rapidly was se- Verde2 p = .210, Verde1 -Verde2 p = .015, Verde1 -Verde1/Verde2 p = .015, Verde2 -Verde1/Verde2 p = .015, Figure 2b,c). The hypersensitivity of srf-2 to Verde1 led to swift extinction of srf-2 and a delayed, but similar trajectory under coinfection, where N2 gained a heightened competitive advantage through Verde1-mediated protection (Figure 2b,c). We did not observe dauer larvae over the course of the evolution experiment.
We hypothesized that host genotype competition in the evolution experiment could be facilitated by either dilution effects or super-spreading by defended genotypes. Our findings broadly recapitulate the results of our evolution experiment, specifically, results of a binomial glm show significant differences in mortality of our treatment groups (χ 2 (11) = 260.33, p < .001), with high mortality when N2 and srf-2 worms were exposed to worms carrying Verde2 and Verde1, respectively (Figure 3, Table S3). We did not find evidence for differential transmission of Verde2 from exposed N2 and srf-2 worms. However, we did find host genotype effects for Verde1 epidemiology, whereby srf-2 mortality is greater when Verde1-exposed N2 hosts are introduced compared to srf-2 (p < .05).

| DISCUSS ION
Microbial symbionts can shape host evolution in a diversity of ways (Brucker & Bordenstein, 2012, 2013Gould et al., 2018;Shapira, 2016;Sharon et al., 2010). The extent to which the strength and outcomes of selection could be impacted when symbionts transition along the mutualist-parasite continuum is unclear (Drew et al., 2021). The global distribution of Leucobacter among plants, animals and free-living environments demonstrates extraordinary niche diversity indicative of a generalist. This broad niche diversity is not uncommon among bacterial genera with similarly wide-ranging host-associations, such as described among Burkholderia (Eberl & Vandamme, 2016) and Pantoea (Walterson & Stavrinides, 2015). We show that within-genus variation in virulence could play a key role in shaping host fitness in C. elegans. Whilst most Leucobacter isolates tested did not have a major killing effect on C. elegans, Verde2 emerged as virulent across a range of wild host genotypes. We further show that coinfection by multiple Leucobacter sequence variants occurs across different nematode host species (Dirksen et al., 2016). The relatively common occurrence of Leucobacter in the wild Caenorhabditis microbiome suggests a consistent association with this host genus in nature. Future work could examine the functional role of different Leucobacter isolates in nematodes to determine how microbial determinants can shape host health.
The combination of variation in virulence and natural coinfection led us to hypothesise that Leucobacter community composition may impact selection of host genotypes. Variation in mortality caused by the different symbiont species, singly and in coinfection, imposed diversifying selection in host populations. Despite the high virulence of Verde2 to wild-type genotypes, the nematode genotype conferring defence (srf-2) does not fix. This result is not due to host dilution effects (Rudolf & Antonovics, 2005) as we found srf-2 does not shield the N2 genotype from parasite acquisition. Conversely, this result could be explained by the cost of the mutant phenotype. Parasite defences are well-established as costly for hosts (Antonovics & Thrall, 1994;Bartlett et al., 2018;Graham et al., 2005;Schmid-Hempel, 2003;Schwenke et al., 2016), although evidence exists for non-costly defence in some instances (Penley et al., 2018).
Low relative srf-2 fitness might be driven by negative covariance between defences against Verde2 and other key physiological components leading to strong directional selection. Specifically, changes to the cuticle structure (altered surface antigenicity) of srf-2 may impact locomotory efficiency and consequently the ability to forage or behaviourally avoid parasites. Similar negative covariance dynamics have been observed between genetic-based immunity and other traits in invertebrate models such as Drosophila, where increased resistance was linked to an inferior competitive ability (Kraaijeveld & Godfray, 1997;Luong & Polak, 2007).
A major challenge associated with evolving genetic-based defences is the cost in the absence of parasite attack (Lochmiller & Deerenberg, 2000;McKean et al., 2008). N2 defences against Verde1 clearly carried lesser cost than srf-2 resistance to Verde2, since N2 rapidly spread to fixation. The purging of srf-2 susceptible genotypes in the Verde1 treatment was probably due to our finding that N2 is super-spreading the parasite within the experimental arena (Lloyd-Smith et al., 2005;Stein, 2011), combined with an overall weaker competitive ability against the N2 genotype.
Superspreading of Verde1 by the N2 genotype may also be an example of symbiont mediated competition, whereby a symbiont that has little cost to one host lineage imposes high costs/death on another, thus benefiting the nonsusceptible symbiont carrier (Murfin et al., 2019). Ultimately, the variable symbiont phenotypes, as well as the fitness constraints and limitations of host genetic-based defences, determined host evolutionary trajectories. Symbiont coinfection can degrade (Ademola & Odeniran, 2016;Leclair et al., 2017) or improve (Bazzone et al., 2008;King et al., 2016) host health relative to single infection. Here, we confirmed that coinfection may help nematodes survive virulent Leucobacter in the environment (Hodgkin et al., 2013). Verde1-mediated protection against Verde2 parasites demonstrates how co-opting protective symbionts could offer an effective and perhaps less costly form of defence (King, 2019;Martinez et al., 2016). Coinfection by Leucobacter species results in the reduction of srf-2 frequency over evolutionary time. Conversely, this defended genotype spreads during exposure to Verde2 parasites. The results suggest that the benefits of the host-based defences in srf-2 are made redundant because of the protection conferred to N2 by coinfection. Martinez et al. (2016) similarly found that symbiont protection caused a host gene conferring resistance to parasites to be at a low frequency in the host population. The decrease in srf-2 frequency during coinfection occurs at a slower rate than during single Verde1 parasite exposure. The speed of loss could be affected by the reduced difference in competitive ability than single parasite exposures (since some N2 hosts succumb to Verde2 killing in coinfection), as well as by a possible dose effect during coinfection.
Divestment of host-based defences has been observed in defensive symbioses across ecological and evolutionary timescales (Ford & King, 2021;Martinez et al., 2016). Although we observed Leucobacter to be widely host-associated in nature, it is often free-living, and the species used herein are probably components of the Caenorhabditis environmental microbiome. As the nematodes disperse or as the environmental microbial community changes, the risks of parasitism may vary across a nematode's lifetime or generations. Forming temporary associations with protective species may therefore subvert the long-term fitness costs of host-encoded defences. This strategy may also minimise costs of symbiont maintenance when infections have been cleared or are absent from the community (Oliver et al., 2008). Such dynamics have been shown experimentally and in the field in the aphid, Acyrthosiphon pisum, whereby abundance of the protective symbiont Hamiltonella defensa has been reported to be associated with seasonal changes in parasitism pressure (Oliver et al., 2008;Smith et al., 2015). Further examination of Leucobactermediated host protection in nature and across a wider range of host genotypes or species will therefore be valuable in understanding the prevalence and evolutionary importance of this symbiosis.
In conclusion, our results show that the variation present in nematode-Leucobacter interactions across the mutualist-parasite continuum affected the speed and trajectory of host evolution.
Given these findings and the expansive geographic and host taxa range of Leucobacter, this genus warrants further consideration as a symbiont that may naturally impact host eco-evolutionary dynamics on a global scale.

ACK N OWLED G EM ENTS
We thank Jonathan Hodgkin for his guidance and for the gifts of Cape Verde Leucobacter isolates. We also thank Marie-Anne Felix for the gifts of other Leucobacter isolates. We are grateful for funding from the Biotechnology and Biological Sciences Research Council