Negotiating mutualism: A locus for exploitation by rhizobia has a broad effect size distribution and context‐dependent effects on legume hosts

Abstract In mutualisms, variation at genes determining partner fitness provides the raw material upon which coevolutionary selection acts, setting the dynamics and pace of coevolution. However, we know little about variation in the effects of genes that underlie symbiotic fitness in natural mutualist populations. In some species of legumes that form root nodule symbioses with nitrogen‐fixing rhizobial bacteria, hosts secrete nodule‐specific cysteine‐rich (NCR) peptides that cause rhizobia to differentiate in the nodule environment. However, rhizobia can cleave NCR peptides through the expression of genes like the plasmid‐borne Host range restriction peptidase (hrrP), whose product degrades specific NCR peptides. Although hrrP activity can confer host exploitation by depressing host fitness and enhancing symbiont fitness, the effects of hrrP on symbiosis phenotypes depend strongly on the genotypes of the interacting partners. However, the effects of hrrP have yet to be characterised in a natural population context, so its contribution to variation in wild mutualist populations is unknown. To understand the distribution of effects of hrrP in wild rhizobia, we measured mutualism phenotypes conferred by hrrP in 12 wild Ensifer medicae strains. To evaluate context dependency of hrrP effects, we compared hrrP effects across two Medicago polymorpha host genotypes and across two experimental years for five E. medicae strains. We show for the first time in a natural population context that hrrP has a wide distribution of effect sizes for many mutualism traits, ranging from strongly positive to strongly negative. Furthermore, we show that hrrP effect size varies across host genotypes and experiment years, suggesting that researchers should be cautious about extrapolating the role of genes in natural populations from controlled laboratory studies of single genetic variants.


| INTRODUC TI ON
Mutualisms between hosts and microbes are ubiquitous and play a critical role in spurring evolutionary innovation and powering ecosystem services. However, we still know little about the genetic variants that influence partner fitness in mutualisms, especially compared to antagonisms (Baskett & Schemske, 2015;Stoy et al., 2020). Understanding which genes affect mutualist fitness, how they are transmitted, and when and how they function can help us predict how different mutualism traits will evolve. For instance, genes residing on mobile genetic elements may sweep through microbial populations more rapidly than vertically transmitted genes (Shapiro, 2016), leading to more rapid evolution of mutualism traits. Genes governing early stages of symbiosis, compared to later stages, may increase a symbiont's host range (Radutoiu et al., 2007) and impact its long-term extinction risk (Koh et al., 2004). Genes with pleiotropic effects may experience stronger evolutionary constraints than genes that only affect single traits (Auge et al., 2019), preventing mutualism traits from reaching optimum values for fitness. Mutualism genes have been uncovered by a variety of methods -including mutant screens, association genetics and quantitative trait locus mapping -that associate mutualism phenotypes with the presence of particular genes or variation among alleles (Burghardt et al., 2017;Gorton et al., 2012;Hu et al., 2020;LaPlante et al., 2021;Piculell et al., 2019;Price et al., 2015;Stanton-Geddes et al., 2013;Torkamaneh et al., 2020). Increasingly, these methods are equipped to detect loci that exhibit context dependent phenotypes, such as mutualistic partner-dependent phenotypes (MacPherson et al., 2018;Wang et al., 2018). However, it remains uncertain how much we can extrapolate from highly controlled laboratory studies of single genetic variants to the function of genes in natural populations.
The impact of a gene on rates of phenotypic evolution depends on the effect size distribution of allelic variants, where 'effect size' indicates how much an allele changes a particular trait value (Dittmar et al., 2016). Random mutations generate alleles with a range of effect sizes (Bataillon & Bailey, 2014;Kassen & Bataillon, 2006), and the width of the effect size distribution in a population provides the genetic variation upon which natural selection acts (Salvaudon et al., 2008;Simonsen & Stinchcombe, 2014). Wider effect size distributions (i.e. more genetic variance) can produce faster responses to selection (Li, 1967), although narrow effect size distributions can accelerate evolution over short timescales (Briggs & Goldman, 2006).
Rates of phenotypic evolution will also depend on the amount of context dependency in the effect size of a candidate gene. Context dependency exists when the effect size of a gene varies with environmental conditions (i.e. phenotypic plasticity) or genotypes at other loci (i.e. epistasis; Remold & Lenski, 2004). High context dependency can alter selection on an allele by limiting the contexts in which it confers effects on fitness. Thus, context dependency can prevent an allele from sweeping through a population, even when selection is strong (Chevin, 2019;Höllinger et al., 2019), or conversely, context dependency can accelerate evolution towards a phenotypic optimum (Borenstein et al., 2006). Given the strong impacts that the distribution and context dependency of effect size can have on the tempo of evolution, it is critical to study these parameters for genes important in mutualisms.
The legume-rhizobium symbiosis is a globally important mutualism that shapes the ecology of wild plant communities (van der Heijden et al., 2006) and provides much of the nitrogen needed in agriculture (Goyal et al., 2021). The symbiosis is initiated when soildwelling rhizobial bacteria infect the roots of leguminous plants, forming nodules in which they fix atmospheric nitrogen into a form plants can use for growth (Poole et al., 2018). In a subset of legumes, differentiation of rhizobia into their nitrogen-fixing form is accomplished by host secretion of nodule-specific cysteine-rich (NCR) peptides, which target the organelle-like structures in which rhizobia are sequestered and trigger rhizobia to undergo partial membrane permeabilisation, genome duplication and loss of reproductive viability inside the nodule (Alunni & Gourion, 2016;Ledermann et al., 2021). The genome of the model legume Medicago truncatula encodes an abundant and diverse family of NCR peptides (Montiel et al., 2017), and plant accessions vary substantially in expression levels of individual NCR peptide genes (Nallu et al., 2014). The role of NCR peptides in rhizobial adaptation to the nodule environment is complex, and most have not been functionally studied. On one hand, some NCR peptides have antimicrobial activity and can kill rhizobia in the nodule, depending on the plant genomic background (Yang et al., 2017). On the other hand, certain NCR peptides are required for rhizobia to persist in nodules, and rhizobia die if hosts fail to express these peptides (Horváth et al., 2015;Kim et al., 2015).
Far from being passive recipients of NCR peptide cues, rhizobia can curtail this signalling mechanism by producing peptidases that degrade specific NCR peptides (Benedict et al., 2021;Price et al., 2015). One such agent is Host range restriction peptidase (HrrP), a plasmid-encoded peptidase discovered in Ensifer meliloti, the rhizobial symbiont of M. truncatula (Crook et al., 2012;Price et al., 2015).
hrrP was uncovered in a mutant screen for rhizobia that gain compatibility with novel hosts: disruption of the hrrP locus allows E. meliloti to fix nitrogen on a host with which it is otherwise incompatible. In some cases, hrrP-expressing rhizobia can avoid differentiating and fixing nitrogen inside nodules, which decreases plant fitness but increases rhizobium fitness (Price et al., 2015). However, the effect of hrrP on host and rhizobium fitness is dependent on host genotype and strain genomic background (Price et al., 2015), with some host genotypes appearing totally resistant to the activity of hrrP (i.e. experiencing normal nitrogen fixation from hrrP-expressing rhizobia), and with some strain genotypes exhibiting different phenotypes even when bearing the same hrrP allele (Price et al., 2015). The context dependency of the effects of NCR peptides and rhizobium peptidases on host and symbiont fitness is captured by the 'working balance' hypothesis of peptidase-NCR peptide dynamics , which predicts that rhizobia and hosts benefit from moderate net NCR peptide levels, but show extreme phenotypes when NCR peptides are excessively low or high. Since hrrP serves as a means for rhizobia to tune their host's level of NCR peptide production, the effect size of hrrP is predicted to vary among hosts to the extent that hosts vary in NCR peptide expression, consistent with (Price et al., 2015). Sequence variation in the hrrP locus could also contribute to variation in hrrP effect size, although only one hrrP allele has been empirically studied to date (Price et al., 2015).
Here, we investigate the distribution and context dependency of hrrP effects on mutualism outcomes. We used a wild collection of the rhizobium E. medicae, which can nodulate several Medicago species including M. polymorpha and M. truncatula (Denton et al., 2007). In a previous PCR screen, 12.5% of the E. medicae strain collection was found to bear hrrP (i.e., hrrP+ strains), with the remainder lacking this locus (i.e. hrrP-strains; Wendlandt et al., 2021). We performed targeted gene disruptions in 12 hrrP+ E. medicae strains to generate hrrP-knockout mutants, and measured hrrP effect size for each strain as the relative difference in trait values of hrrP+ and hrrP-strains. Because we used wild E. medicae strains, this measure of hrrP effect size could be influenced by variation among strains in hrrP allelic identity (affecting its specificity and catalytic activity for targeted NCR peptides), hrrP expression level (due to variation in promoter sequence or other modifiers of gene expression), and other loci with epistatic effects. This measure of hrrP effect size is intentionally broad in scope to capture natural phenotypic consequences of disruption to this locus in nature, making it an ecologically relevant way to understand hrrP effect size. To assess how variation among rhizobia strains, plant host genotypes, and environments alters the impact of hrrP on mutualism outcomes, we asked: (

| Rhizobia strains and inocula preparation
In the Knockout Experiment, we used 12 E. medicae strains (Table 1) genotyped as hrrP+ by Wendlandt et al. (2021). These 12 strains span the genetic diversity uncovered for hrrP by Wendlandt et al. (2021) and represent 6 unique hrrP sequences for the partial coding region for which sequence data are available (Table S2). We generated one hrrP-knockout mutant strain from each of the 12 E. medicae strains using homologous recombination insertional mutagenesis. Briefly, a 3 kb non-replicative plasmid encoding a neomycin resistance gene was inserted into the hrrP coding region, and the presence of the insert was verified by testing for neomycin resistance and performing PCR with primers whose product spans the gene-insert junction (Price et al., 2015). Previous tests of mutants made in this way found no pleiotropic effects of neomycin insertion (Paul Price, pers comm).
In the G × G Knockout Experiment, we used five of the wild-type hrrP+ strains used previously as well as their knockout hrrP-derivatives (Table 1).
We prepared rhizobial inocula for the greenhouse experiments by streaking frozen glycerol stocks of each wild-type and knockout strain onto tryptone yeast (TY) agar plates and incubating until single colonies formed. Before preparing inocula, we confirmed that hrrP could be PCR-amplified from each wild-type strain in the upcoming experiment. Single colonies were then used to inoculate 1 ml of aliquots of TY broth, which were incubated at 30°C and 300 rpm for 3 days. Two hundred and fifty microlitres of the 1-mL culture was used to inoculate 4.75 ml of TY broth, which was incubated at 28°C and 300 rpm for 2 days. The OD 600 was measured for each culture to estimate the number of colony-forming units (CFUs) using a conversion factor of CFU ml −1 = 5.8 × 10 7 × OD 600 .
Cells were pelleted, separated from supernatant and resuspended in 0.1X TY broth to concentrations of approximately 10 6 CFU ml −1 . (2009), we assumed that the relationship between OD 600 and CFU was approximately consistent across strains. Since our main goal was to inoculate plants with enough rhizobia (10 6 CFU) that nodule formation would not be limited by the number of rhizobia present, moderate fluctuations in the relationship between OD 600 and CFU among strains should have weak impacts on our findings. See Table S1 for specific methods used in each experiment.

| Plant host genotypes and growth conditions
We generated all seeds in a common garden in greenhouse conditions to minimise maternal effects. In the Knockout Experiment, we used one M. polymorpha genotype (RTM;  (Vincent, 1970) containing 500 µM of NH 4 NO 3 , beginning the week after inoculation. Plants were fertilised weekly for a total of five times throughout the experiment; see Table S1 for details. In total, we measured up to six traits from each experiment. As proxies of plant fitness, we measured leaf count, dry shoot mass and dry shoot mass per nodule (reflecting the balance of benefits to plants versus rhizobia). As proxies of rhizobium fitness, we measured total nodule count, nodule size and CFU per nodule (the latter was log-transformed before analysis). Pairwise correlation coefficients for these responses are reported in Table S3

| Statistical analysis
We analysed data using general linear mixed models implemented with lme4 v.

| The effect size of hrrP varies among strains of wild rhizobia
We found that hrrP had positive effects on many proxies of plant host and rhizobium fitness. On average, hrrP increased leaf count by 8% ( Figure S1), increased shoot mass by 11% (Figure 1a), decreased shoot mass per nodule by 6% (Figure 1b), increased nodule count by 34% (Figure 1c), and increased logCFU per nodule by 7% ( Figure 1d) Figures 1 and S2). Model 2 tested for effects of host genotype on hrrP effect size using 5 E. medicae strains and 2 M. polymorpha host genotypes (G × G Knockout Experiment, Figure 2). Model 3 tested for effects of experiment year on hrrP effect size using 5 E. medicae strains and one M. polymorpha host genotype (pooled Knockout Experiment and G × G Knockout Experiment, Figure 3). For each model and response variable, n indicates the number of plants used in the analysis. ***p < 0.0001, **p < 0.001, *p < 0.05, †p < 0.10.  Figure 2).

| The effect size of hrrP varies among experiment years
For the 5 E. medicae strains tested in two experimental years, mean hrrP effect size tended to be smaller or more negative in 2019 compared to 2018 (Model 3: "Year"; Table 2). On average,  Table 2). hrrP had inconsistent effects on shoot mass between years for three strains (E. medicae AZN234, PEA63 and PEA143; Figure 3a). Since the wild-type E. medicae

| DISCUSS ION
Predicting the evolutionary dynamics of genes involved in mutualism requires that we understand how these genes contribute to standing genetic variation in natural populations and the degree of context dependency of their phenotypic effects. However, this information is lacking for most loci impacting fitness in mutualism, particularly loci that show large effects in controlled laboratory experiments.
Our study reveals that the effects of hrrP, a horizontally transmitted, plasmid-borne locus that can have major effects on the fitness of both mutualist partners, are highly genetically and environmentally context dependent in a set of wild rhizobia strains. We find that: (1) the effect size of hrrP on symbiotic partner fitness can differ in sign and magnitude among wild rhizobia strains, and that hrrP effect size shows context dependency between (2)  the range of what we uncovered in our wider survey of strains. We also saw that hrrP increased nodule count more than shoot mass, such that hrrP decreased shoot mass per nodule and shifted the balance of symbiotic benefits towards rhizobia. This could be a subtle form of exploitation within the constraints of the working balance hypothesis, whereby hrrP evolves to increase benefits to rhizobia more than it increases benefits to plant hosts (Klein et al., 2022).
Although previous research uncovered a single hrrP allele of large effect (Price et al., 2015), we find that hrrP has small or neutral effects on symbiotic traits in many E. medicae strains. This finding suggests that researchers should be cautious about interpreting the role large-effect genes will have in nature until more variants of that gene have been studied. Although in certain contexts, hrrP can be a strong driver of variation in host and symbiont fitness, our data suggest that large-effect hrrP alleles such as the B800 allele (Crook et al., 2012;Price et al., 2015) may not be common in natural populations, where many hrrP alleles have smaller effects on symbiosis traits. One mechanism that could favour small-effect hrrP loci is suggested by the working balance hypothesis, which predicts that hosts and symbionts experience selection towards similar net levels of NCR peptide activity and thus could experience fitness alignment (Friesen, 2012), despite the antagonistic effects of peptidases

| Conclusions
Both mutualists and pathogens can have large fitness effects on their hosts, but we generally know more about the genes underlying pathogen interactions than the genes underlying mutualistic interactions.
Since mutualisms differ from antagonisms in that partners coordinate to exchange a service or resource, the genetic basis of interaction outcomes may be fundamentally more complex for mutualisms than for antagonisms (Stoy et al., 2020). In line with this anticipated complexity, we show that hrrP from a panel of wild E. medicae strains has a wide range of effect sizes on mutualism outcomes for legumes and rhizobia, and that these effect sizes are highly context dependent. These findings are consistent with the "working balance" model of peptidase-NCR peptide activity, in which the fitness value of hrrP depends on host levels of NCR peptide production. The high context dependency of hrrP effects could also contribute to the evolutionary stability of mutualism by preventing genes of large effect from sweeping through symbiont populations. Furthermore, we highlight the importance of measuring effect sizes and context dependency of multiple variants of candidate mutualism genes to understand their probable role in the evolution of wild populations.