The Potential for Genotype-by-Environment Interactions to Maintain Genetic Variation in a Model Legume-Rhizobia Mutualism.

The maintenance of genetic variation in mutualism-related traits is key for understanding mutualism evolution, yet the mechanisms maintaining variation remain unclear. We asked whether genotype-by-environment (G×E) interaction is a potential mechanism maintaining variation in the model legume-rhizobia system, Medicago truncatula-Ensifer meliloti. We planted 50 legume genotypes in a greenhouse under ambient light and shade to reflect reduced carbon availability for plants. We found an expected reduction under shaded conditions for plant performance traits, such as leaf number, aboveground and belowground biomass, and a mutualism-related trait, nodule number. We also found G×E for nodule number, with ∼83% of this interaction due to shifts in genotype fitness rank order across light environments, coupled with strong positive directional selection on nodule number regardless of light environment. Our results suggest that G×E can maintain genetic variation in a mutualism-related trait that is under consistent positive directional selection across light environments.

Introduction

Characterizing the forces that maintain genetic variation in the face of ongoing selection is a long-standing, classic question in evolutionary biology (Walsh and Blows, 2009). For mutualisms—interspecific interactions that benefit both species—answering this question is especially challenging. By definition, the traits that promote mutualisms are positively related to fitness, and thus variation for fitness-related traits should be reduced (Charlesworth, 1987; Kruuk et al., 2000; Heath and Stinchcombe, 2014). Moreover, the mechanisms thought to promote stability and maintenance in mutualisms are themselves predicted to reduce variation in mutualist quality (Heath and Stinchcombe, 2014). One straightforward mechanism capable of maintaining variation in mutualisms, like all traits, is genotype-by-environment (G×E) interaction (Gillespie and Turelli, 1989), although we have few examples of whether G×E maintains variation in mutualism-related traits.

Genotype-by-environment interactions can occur either through changes in the rank order of genotype fitness or in the magnitude of genetic variance across environments (Bowman, 1972; Muir et al., 1992; Hühn et al., 1993; Des Marais et al., 2013). Only changes in the rank order of genotype fitness can contribute to the maintenance of genetic variation, because genotypes that perform well in one environment are not the same genotypes that perform well in another (Mitchell-Olds, 1992; Johnson, 2007). In contrast, G×E interactions that are mostly due to changes in genetic variance might affect the response to selection across environments (Fisher, 1958) but will not effectively contribute to the maintenance of genetic variation. Several studies have tested for G×E interactions in a variety of symbiotic systems, such as plant–fungi systems (Ahlholm et al., 2002; Werner et al., 2018), the wasp–Wolbachia symbiosis (Mouton et al., 2007), ant–plant mutualisms (Abdala-Roberts et al., 2012; Abdala-Roberts and Mooney, 2013), and legume–rhizobia mutualisms (Heath and Tiffin, 2007; Heath et al. 2010; Ossler and Heath, 2018; Batstone et al., 2020; Heath et al., 2020; Magnoli and Lau, 2020 ). Most of these studies concluded that G×E could contribute to genetic variation in mutualism-related traits and used crossing reaction norms as qualitative support. However, only two of these studies have explicitly tested the maintenance of genetic variation by including a quantitative assessment of the relative contribution of genotype rank order shifts versus changes in genetic variance across abiotic environments (e.g., Batstone et al., 2020; Heath et al., 2020).

Genetic variation in mutualisms could also be maintained through environmental heterogeneity in selection (Mitchell-Olds et al., 2007; Heath and Stinchcombe, 2014). Shifts in resource availability can alter the costs, benefits, and ultimately the strength of the mutualistic interaction, which could potentially cause changes in the magnitude and/or direction of selection on mutualism-related traits across environments (Thrall et al., 2007; Batstone et al., 2018). Changes in nutrient availability have been shown to reduce the net benefits for hosts in a variety of mutualistic systems, including coral, lichen, plant–fungi, and legume–rhizobia (Shantz et al., 2016). For example, in some legume–rhizobia systems, nitrogen enrichment decreased the net benefits provided by rhizobia to plant hosts, causing reduced investment by hosts in the mutualism (Regus et al., 2015, 2017; Weese et al., 2015). Such a reduction in mutualism investment could lead to relaxed selection by hosts for mutualistic partners, resulting in increased mutualism variation (or slowing the loss of variation) over time. Van Cauwenberghe et al. (2015) showed that long-term nitrogen addition led to increased rhizobial genetic diversity in the legume–rhizobia system Vicia cracca–Rhizobium leguminosarum. In addition, Simonsen and Stinchcombe (2014b) found a reduction in fitness costs for legumes that hosted ineffective rhizobia in the presence of herbivory, which could directly lead to reduced selection against ineffective partners and maintain mutualism variation. Batstone et al. (2020) explicitly tested for environmental heterogeneity in selection for a legume–rhizobia mutualism through a selection-by-environment (S×E) interaction analysis and found that the magnitude of selection for a mutualism-related trait varied across greenhouse and field plot conditions.

A common and useful model for studying the effects of genetics and the environment on the maintenance of mutualism variation is the legume–rhizobia mutualism. It is a resource-based mutualism, with leguminous plants providing carbon and shelter in root nodules, and mutualistic root-dwelling bacteria providing fixed nitrogen (Kiers et al., 2003; Gorton et al., 2012). There is a vast amount of information on the mechanisms regulating legume–rhizobia mutualisms—from genomic, cellular, to the individual level—making it a powerful system for studying mutualism evolution (Jones et al., 2007). In addition, abiotic and biotic factors can be effectively manipulated to shift the costs and benefits of the interaction due to resource availability (Bronstein, 1994; Heath and Stinchcombe, 2014). However, environmental manipulations of the legume–rhizobia mutualism tend to be largely one-sided, with most studies focusing on effects of soil nitrogen (N) availability due to its ecological and agricultural importance in artificial fertilizer use (e.g., Heath, 2010). The results of these studies tend to show variable outcomes for G×E, with some finding a significant G×E interaction with varying N availability (Heath and Tiffin, 2007; Heath, 2010), and others showing that plant hosts either consistently associated with (Grillo et al., 2016) or discriminated against (Regus et al., 2014) rhizobial strains regardless of N environment. The availability of carbon (C) resources is also equally important for understanding the relationship between legumes and rhizobia, but literature addressing this is scarce. Elevated CO2 (eCO2) has gained recent interest, but its effects on legume–rhizobia mutualisms also show inconsistent results depending on host species or genotype (Lüscher et al., 2000; West et al., 2005), rhizobial genotype (Bertrand et al., 2007), duration of eCO2 exposure (Hungate et al., 2004; Rogers et al., 2006), or the surrounding environment (West et al., 2005). However, the results on eCO2 do indicate that partner genotype and C availability have the potential to influence mutualism efficiency and genetic variation for mutualism-related traits.

An alternative way to test for the influence of C availability on mutualisms is through shifts in light availability. The physiological effects of light condition on plants are well documented (e.g., Evans, 1989; Cronin and Lodge, 2003), and individuals of most plant species show variation in their plastic responses to light (Winn and Evans, 1991; Sultan and Bazzaz, 1993). In addition, both G×E and adaptive plastic responses to light quality and quantity are well known (Donohue et al., 2000; Schmitt et al., 2003; Heschel et al., 2004; Stinchcombe et al., 2010). If legumes respond differently to varying light conditions like most plant species, it could affect how legume hosts allocate C resources and ultimately associate with rhizobia under varying light conditions (Schwartz and Hoeksema, 1998). A few studies have tested the effects of light availability on legume–rhizobia mutualisms using one-to-one pairings, but with variable results. Some of these studies suggested that maintaining the association becomes costly for plant hosts under low light, resulting in a shift from mutualism to parasitism (Lau et al., 2012; Ballhorn et al., 2016), while another suggested that the mutualistic association becomes commensal (Friesen and Friel, 2019). Heath et al. (2020) explicitly tested for G×E effects on legume–rhizobia mutualisms using legumes of a single genotype inoculated with one of 11 strains of rhizobia and two light environments. They showed a rhizobia genotype-by-light environment interaction for plant performance traits, but not for mutualism-related traits such as nodule number. Using several legume genotypes might more directly test for genetic variation in resource allocation and mutualism association under varying light conditions to effectively determine whether G×E interactions maintain genetic variation in mutualism-related traits.

We manipulated light availability, and hence the amount of C available for plants to provide to rhizobia, in an experiment testing for G×E in the model legume–rhizobia system, Medicago truncatula–Ensifer meliloti. Specifically, we asked the following questions: (1) Do legumes form fewer mutualistic associations with bacteria in low light (low C) environments? (2) Is there G×E for plant performance and mutualism-related traits in response to low light, and are these interactions capable of maintaining genetic variation? And (3) How does natural selection act on morphological and mutualistic traits in contrasting light environments? Here, we present the results of an experiment testing for G×E interactions in a plant–bacteria mutualism, quantitatively examining their potential to maintain variation in mutualism-related traits.

Results

Poor germination and establishment led to a final sample of 666 plants that survived until the end of the experiment. While this germination and establishment fraction was lower than our past work with Medicago (and these lines), we verified that the total fraction of germinants was not correlated with the focal traits analyzed below (|r| < 0.19, P > 0.17), suggesting that germination and establishment variation did not bias our estimates of plant performance or nodulation. The experiment was terminated at 9 weeks to prevent nodule senescence (Puppo et al., 2005; El Msehli et al., 2019), at which point 409 (61%) plants had flowered. We found 25 ± 6.50 (mean ± S.E.) nodules in the uninoculated control plants, which suggests some rhizobial transfer among individuals within each block. However, this will not affect the results since all experimental plants received the same inoculation treatment. We excluded uninoculated control plants in the analysis of our data.

Plant Genotype and Light Environment Affect Plant Growth

The productivity and growth of plant hosts were estimated via leaf number, and aboveground and belowground biomass. Plants growing in shaded conditions produced marginally fewer number of leaves than ambient plants, showing a ∼30% decrease in leaf number from ambient to shade conditions (ambient 7.39 ± 0.44 versus shade 5.18 ± 0.17, Table 1 and Figure 1A). However, plant genotype had a significant effect on leaf number (Table 1), with the average number of leaves at harvest ranging from three to ten leaves across genotypes. There was no significant G×E for leaf number (Table 1 and Figure 1B).

Table 1

Summary Statistics for the Effects of Light Treatment and Plant Genotype on Plant Performance and Mutualism-Related Traits.

	Leaf Number	Aboveground Biomass	Belowground Biomass	Aboveground:Belowground Ratio	Nodule number
Fixed effects	F(1,9.2944)	F(1,13.209)	F(1,13.281)	F(1,9.0881)	F(1,11.67)
Light treatment	3.75	18.33∗∗∗	19.16∗∗∗	19.42∗∗	23.19∗∗
Random effects	χ2df=1	χ2df=1	χ2df=1	χ2df=1	χ2df=1
Plant genotype	14.27∗∗∗	5.16∗	3.88∗	2.76∗	8.13∗∗
GxE	0.63	13.22∗∗∗	32.45∗∗∗	0.20	5.42∗

Linear mixed models included light treatment as a fixed effect, and plant genotype and genotype-by-treatment interaction (G×E) as random effects. Analyses included ANOVA F(num df, denom df) for fixed effects and log-likelihood ratio tests (χ2) for random effects. Significant values are in bold and coded as: ∗∗∗p < 0.001, ∗∗p < 0.01, ∗p < 0.05.

Figure 1

Effects of Light Treatment and Genotype-by-Treatment Interaction on Plant Performance and Mutualism-Related Traits.

(A, C, E, G, and I) The main effect of light treatment (shade and ambient), with points representing raw treatment means ± 95% confidence intervals (CIs).(B, D, F, H, and J) The genotype-by-treatment interaction effect, with points representing least-squares line means ± 95% CIs.

Plants that received full light were significantly larger on average than plants that were shaded. Aboveground biomass decreased by ∼61.5% from ambient to shaded conditions (ambient 0.22 ± 0.024 g versus shade 0.085 ± 0.0086 g, Table 1 and Figure 1C), but was significantly affected by plant genotype and G×E (Table 1 and Figure 1D). Belowground biomass also decreased by ∼72% from ambient to shade (ambient 0.083 ± 0.0084 g versus shade 0.023 ± 0.0012 g, Table 1 and Figure 1E), and showed a significant genotype and G×E effect (Table 1 and Figure 1F). The relative contribution of genotype rank order shifts in the G×E interaction was ∼31.6% for aboveground biomass and ∼12.5% for belowground biomass (Supplemental Table 1), suggesting that a majority of the G×E was due to changes in the magnitude of genetic variance between environments; specifically, an increase in variance in the ambient treatment compared with the shade treatment (Figure 1D and 1F).

We also analyzed the ratio of aboveground to belowground biomass as an index of plant investment into different components of growth. Light availability significantly affected this ratio, but the response was in the opposite direction of the other traits. Aboveground:belowground ratio increased by ∼50.3% from ambient to shade (ambient 2.470.11 g versus shade 3.710.24 g, Table 1 and Figure 1G), indicating increased allocation to aboveground growth in the light-limited conditions. There was a marginal effect of plant genotype and no significant treatment-by-genotype interaction on aboveground:belowground biomass (Table 1 and Figure 1H).

Plant Genotype–Light Environment Interaction Affects Nodule Number

Nodule number showed a significant main effect of light treatment and genotype, as well as a significant G×E interaction effect (Table 1). Plants in the shade had about 20 fewer nodules on average than in full light (∼63% decrease, ambient 32.05 ± 2.51 versus shade 11.85 ± 0.73, Figure 1I). Approximately 83.3% of the significant G×E interaction was explained by a change in genotype rank order between treatments (Supplemental Table 1). Most, but not all, plant genotypes showed higher average nodule number in the ambient treatment than shade, and the magnitude of this plastic response was variable among genotypes (Figure 1J).

Genetic Correlations between Traits and across Environments

There was a significant positive genetic correlation between most traits within each light treatment (Figure 2). Our focal mutualism-related trait, nodule number, was significantly positively correlated with all plant performance traits in both ambient and shade environments, except the aboveground:belowground ratio, which showed no significant correlation. When comparing across light treatments, all traits showed a significant positive correlation between the shade and ambient environments ranging from 0.34 for aboveground:belowground ratio to 0.76 for belowground biomass (Figure 2).

Figure 2

Genetic Correlation Matrix between Traits and across Environments.

Correlations are of raw line means for each trait within the ambient (gold) and shade (green) treatments. Dark bordered boxes show between treatment correlations for each trait. Numbers present correlation coefficients, which are also represented by the color legend on the right. Insignificant correlations (p > 0.05) are indicated by the absence of an ellipse.

Selection for Nodule Number Does Not Differ between Light Environments

All three methods of our selection analyses (absolute, global, and local) revealed a significant S×E interaction for leaf number, while showing that selection for nodule number did not differ between light environments (Table 2 and Figure 3). The global and local analyses showed inconsistencies for belowground biomass, with a non-significant S×E interaction in the global analysis and significant S×E in the local analysis (Table 2 and Figure 3F and 3I, respectively). Estimates of the selection gradient (β) from the local analysis indicated positive directional selection for leaf number in both environments, but the magnitude of selection was significantly greater in the shade than ambient (Table 3 and Figure 3H). For nodule number, both the magnitude and direction of selection remained similar across light environments (Table 3 and Figure 3G).

Table 2

Selection-by-Environment Interactions Using Three Methods of Scaling Fitness (Aboveground Biomass) and Traits.

Model	Term	Sum of Squares	F(1,92)
Absolute	Intercept	0. 018206	17.7636∗∗∗
	Treatment	0. 002129	2.0772
	Nodule number	0. 006355	6.2011∗
	Leaf number	0. 014393	14.0437∗∗∗
	Belowground biomass	0. 065165	63.5825∗∗∗
	Nodule × treatment	0. 000000	0.0002
	Leaf × treatment	0. 006106	5.9578∗
	Belowground × treatment	0. 000700	0.6834
Global	(Intercept)	23.0933	539.2206∗∗∗
	Treatment	0.8351	19.4995∗∗∗
	Nodule number	0.2656	6.2011∗
	Leaf number	0.6015	14.0437∗∗∗
	Belowground biomass	2.7231	63.5825∗∗∗
	Nodule × treatment	0.0000	0.0002
	Leaf × treatment	0.2552	5.9578∗
	Belowground × treatment	0.0293	0.6834
Local	(Intercept)	100.000	1456.8018∗∗∗
	Treatment	0.000	0.0000
	Nodule number	0.571	8.3197∗∗
	Leaf number	0.982	14.3004∗∗∗
	Belowground biomass	8.684	126.5117∗∗∗
	Nodule × treatment	0.019	0.2830
	Leaf × treatment	0.449	6.5364∗
	Belowground × treatment	0.854	12.4383∗∗∗

ANCOVA F(num df, denom df) are presented for three linear models: absolute fitness, globally (across light treatments), and locally (within light treatment) relativized fitness and standardized traits. Significant values are in bold and coded as: ∗∗∗p < 0.001, ∗∗p < 0.01, ∗p < 0.05.

Figure 3

Partial Regression Plots of Selection-by-Environment Analyses for Plant Performance and Mutualism-Related Traits.

Aboveground biomass is used as the fitness proxy, and three scaling methods were applied: absolute (A–C), global (among treatments, D–F), and local (within treatments, G–I). The y axes of all plots were transformed back to their original scale by adding mean relative fitness, and the x axes of (A)–(C) were similarly transformed by adding mean trait values to better reflect the ANCOVA results in Table 2. Points represent raw line means in the ambient (gold, triangle) and shade (green, circle) treatments that were scaled based on the three scaling methods. Solid lines and shading represent selection gradients and ±95% CIs, respectively. Corresponding within-treatment β estimates are provided for all traits, with bolded values representing traits that showed a significant (p < 0.05) selection-by-environment interaction.

Table 3

Summary Statistics of Local Scale (within Light Treatment) Selection Analysis.

Treatment	Term	β	SE	t value	Sum of Squares	F(1,46)
Ambient	Nodule number	0.143	0.037	3.85∗∗∗	0.352	14.81∗∗∗
	Leaf number	0.037	0.026	1.40	0.046	1.95
	Belowground biomass	0.623	0.039	16.04∗∗∗	6.114	257.34∗∗∗
Shade	Nodule number	0.099	0.071	1.39	0.220	1.94
	Leaf number	0.191	0.051	3.69∗∗∗	1.548	13.63∗∗∗
	Belowground biomass	0.326	0.067	4.82∗∗∗	2.641	23.27∗∗∗

The summary presents selection gradient (β) estimates with standard errors (SE) and t values, as well as ANCOVA sum of squares and F(num df, denom df). Significant values are in bold and coded as: ∗∗∗p < 0.001, ∗∗p < 0.01, ∗p < 0.05.

Discussion

Explaining the maintenance of variation in mutualisms remains a largely unresolved topic in evolutionary ecology, with several potential mechanisms that require testing (Heath and Stinchcombe, 2014). We used a manipulative experiment to test for G×E in mutualism-related traits, driven by changes in light and hence carbon availability. We found predicted effects of reduced plant performance (leaf number, and aboveground and belowground biomass) and mutualism-related traits (nodule number) under low light conditions, along with G×E interactions for nodule number driven by appreciable changes in rank rather than variance. Our selection analyses suggest that nodule number is under strong positive directional selection in each environment, and this was consistent across all three methods of scaling: absolute, global (across environment), and local (within environment). There was little evidence of environmental heterogeneity in the magnitude or direction of selection for nodule number. Below, we discuss in detail two specific aspects of our main results. First, we evaluate the consequences of G×E being driven by changes in rank order and strong positive selection for nodule number on the evolution of mutualism-related traits. Second, we discuss whether plants produced fewer nodules in low light because they were smaller on a whole-plant level or as an active response against decreased carbon availability, and the implications of this result for the evolution of the mutualism.

G×E and Selection for Nodule Number

Our findings on nodule number illustrate how G×E can contribute to the maintenance of variation in mutualisms. Nodule number showed significant G×E, with the majority being driven by changes in rank order rather than variance. It is important to note, however, that this is not simply a function of the amount of G×E and the across environment correlation. Comparing nodule number and aboveground biomass, nodule number has a higher correlation across environments (i.e., less G×E: Cockerham, 1963), yet a much higher fraction of its G×E is due to changes in rank than variance (83% versus 31%). Nodule number was also under positive selection in each treatment, with little evidence of changes in selection depending on light availability. Consequently, genotypes that produce more nodules are strongly favored in each condition, although which genotypes are most strongly favored differs depending on light availability.

The contribution of genotype rank order shifts in nodule number variation across light environments suggests that some host genotypes invest more in the mutualism in one environment and reduce investment in the other. In addition, there was no significant shift in the magnitude or direction of selection for nodule number between light environments at both the global and local scale (Table 2); there was consistent positive directional selection in both light environments (Table 3 and Figure 3). A lack of S×E interaction coupled with a significant G×E due to genotype rank shifts can maintain genetic variation in nodule number because selection can favor plant genotypes with high nodule production in both light environments, but the genotypes that have high nodule numbers in one environment are not identical to those with high nodule numbers in the other. Furthermore, the genetic correlation for nodule number across environments was strongly positive (but less than 1: Figure 2), which indicates that the groups of genotypes that produce more or less nodules will be similar across environments, but that G×E with rank order shifts can affect which individual genotypes are selected for or against in each environment. Thus, a genotype rank order shift between light environments under positive directional selection could maintain genetic variation in the legume–rhizobia mutualism.

To precisely determine whether G×E is effective at maintaining genetic variation in the legume–rhizobia mutualism, it is also important to consider the frequency of selective environments experienced by natural populations (Gomulkiewicz and Kirkpatrick, 1992). Determining the natural frequency of selective environments proves to be challenging but it is critical for understanding whether the environmental context dependency we found is realistic. The evolutionary response of plastic traits can change drastically with the frequency of selective environments (McDonald and Ayala, 1974; Via and Lande, 1985; Gomulkiewicz and Kirkpatrick, 1992; Arnold and Peterson, 2002; Kelley et al., 2005; Stinchcombe et al., 2010). Our results indicate that plant hosts will experience positive directional selection for nodule number regardless of light environment, but the frequency of light environments will determine which genotypes will be selectively favored, and whether genetic variation will be maintained. For example, if individuals experience full light conditions more frequently than shade, the effects of G×E will no longer be a factor maintaining genetic variation; selection will act on genotypes only as they are ranked in the ambient treatment. Moreover, because there is more genetic variation expressed under ambient light than reduced light, we would expect larger responses to selection in those environments. Natural populations of Medicago truncatula experience a wide range of habitats across the Mediterranean (Batallion and Ronfort, 2006); however, information on the frequency of light conditions experienced by these populations is limited. Future studies must consider how frequently legume hosts experience varying light conditions in nature to determine whether G×E can effectively maintain genetic variation in the legume–rhizobia system.

Reduced Nodule Number under Low Light: Causes and Implications

We found a significant average decrease in the number of nodules—the sites of resource exchange in the mutualism and our proxy for the prevalence of the mutualistic interaction (Heath and Tiffin, 2009). There was a ∼63% reduction in nodule number from the ambient to shade treatment (Figure 1I), which could either be part of a whole-plant level response to low light availability or an active response by plant hosts to reduce mutualistic interactions. A whole-plant level response is possible since we also found an average reduction in our plant performance traits from the ambient to shade treatment: leaf number, aboveground biomass, and belowground biomass, by approximately 30%, 61.5%, and 72%, respectively (Figure 1). Low light availability affects photosynthetic activity, respiration rates, and resource exchange and allocation, which can drastically influence whole-plant level traits, such as leaf number and biomass (Corré, 1983; Givnish, 1988; Sultan and Bazzaz, 1993; Stuefer and Huber, 1998; Cronin and Lodge, 2003; Casal, 2012). However, our statistical models for nodule number included total biomass as a covariate, and we still found a significant average reduction under low light conditions, indicating that reduced allocation to rhizobia might be an active response to low light and low carbon availability.

A reduction in nodule number from ambient to shade could suggest that plant hosts invested less into the mutualistic interaction in low light conditions. The possibility of reduced mutualistic investment in low light is consistent with classic theoretical work that predicts shifts in mutualistic interactions based on the availability of traded resources (Trivers, 1971; Axelrod and Hamilton, 1981; Schwartz and Hoeksema, 1998; Sachs and Simms, 2006; Nathaniel Holland and DeAngelis, 2009). The carbon resources that plant hosts provide to their rhizobial partners may become less available in shaded conditions due to reduced photosynthetic activity and shifts in resource allocation from belowground to aboveground growth (e.g., Evans, 1989; Cronin and Lodge, 2003). Our experimental legumes showed a greater aboveground:belowground biomass ratio in the shade compared with full light conditions (Figure 1G), suggesting that resources in the shade may have been diverted toward aboveground growth as opposed to belowground, which could lead to reduced nodule production under light-limited conditions. A reduction in nodule number could also be a response by the host to minimize costs of the mutualistic association due to reduced plant performance at low light. Studies on plant–fungal (Zheng et al., 2015) and legume–rhizobial (Regus et al., 2015) systems at varying light conditions showed that, as the mutualistic benefits provided to hosts decreased at low light, so did host investment in the mutualism, which could weaken the interaction under low light availability.

A more precise assessment of whether light availability has a direct effect on the mutualism would require measuring the amount of costly photosynthetic carbon allocated to rhizobia and the relative exchange of carbon and nitrogen resources in shaded versus ambient environments. Nevertheless, our results suggest that legume hosts interact less with their rhizobial partners under low light, regardless of whether this is a direct or indirect response. The prevalence of the mutualism is reduced under low light and carbon conditions, but whether this response has any substantial effect on the evolution of the mutualism depends on how frequently and consistently the system experiences low light conditions.

Methods

Study System

Medicago truncatula Gaertn. is an annual legume species native to the Mediterranean (Bataillon and Ronfort, 2006). It has a short generation time, produces large uniform nodules, and a small genome size, with a growing availability of genomic tools and methods for studying legume–rhizobia mutualisms (Barker et al., 1990; Blondon et al., 1994; Cook, 1999; Stanton-Geddes et al., 2013). It is also a highly selfing species, which makes it possible to grow nearly genetically identical individuals of a single accession in multiple environments, facilitating studies of G×E. We used 50 M. truncatula genotypes supplied by the French National Institute for Agricultural Research in Montpellier, and the United States Department of Agriculture—Agricultural Research Services in Washington. We co-inoculated the legumes with two strains (Em1021 and Em1022) of the mutualistic rhizobial species, Ensifer meliloti. These strains differ in their ability to fix N (Batstone et al., 2017), making the rhizobial environment relatively more complex and realistic than a single strain inoculation. We used lab stocks of Em1021, while Em1022 was supplied by Batstone et al. (2017).

Experimental Design

We conducted a manipulative greenhouse experiment in the Earth Sciences Centre at the University of Toronto. We applied two light treatments: ambient and shade, where plants in the ambient treatment were exposed to normal greenhouse conditions (16:8 h light:dark cycle, 22°C day and 18°C night temperatures), and those in the shade treatment were covered with neutral shade cloth that blocked 70% of light. The experiment was set up as a split-plot randomized design, where each block contained both the ambient and shade treatments separated into two racks. We used one individual per genotype, placed in random locations in each rack, for a total of 50 experimental individuals per rack. To test for contamination among plants, an uninoculated control from an extra genotype was included in each rack. The experiment was replicated using 10 blocks, for a total of N = 1020 plants (1000 experimental plants and 20 contamination control plants).

Medicago truncatula Planting and Data Collection

We prepared M. truncatula seedlings following standard protocols (Garcia et al., 2006; Simonsen and Stinchcombe, 2014a; Wood et al., 2018). We scarified seeds and sterilized them in bleach and 95% ethanol. We then imbibed seeds in distilled water for 30 min before placing them in 1% agar plates for stratification at 4°C for approximately 13 days.

Once stratified, we incubated the seeds in the dark at room temperature (∼20°C) for approximately 24 h to allow for radical elongation, then exposed them to light for another hour to initiate chlorophyll production. The germinated seedlings were then moved to the greenhouse and planted into individual autoclaved 120 ml Cone-tainers plugged with polyester fiber (pillow stuffing) and filled with sand (Wood et al., 2018). Within 2 days of planting, the seedlings were co-inoculated with a mixture of Ensifer meliloti strains, Em1021 and Em1022.

Seedlings were initially misted with distilled water and exposed to full light for the first 2 weeks, regardless of light treatment, to promote seedling establishment. Once the first true leaflets emerged, we switched from top watering to bottom watering and the shade treatment was applied to half of the plants (Lau et al., 2012). Shade-treated plants were covered with a tent made from a 24 × 12 × 7 inch wood frame wrapped in shade cloth and fastened by zip ties. Seedlings that did not successfully germinate or survive were removed from the experiment.

After 6 weeks of growth, we measured the number of leaflets per plant, and by 9 weeks we harvested the plants by removing them from their Cone-tainers and separating the aboveground structures (shoots and leaves) from the belowground structures (roots and nodules). We placed the aboveground portion of plants in a drying oven for about 2 weeks to ensure they were fully dried before weighing them for aboveground biomass. We placed the belowground portion in sealed plastic bags to prevent desiccation and stored them at 4°C for counting and collecting rhizobial nodules. We then dried and weighed the belowground portion to estimate belowground biomass.

Rhizobia Inoculation and Data Collection

We prepared and inoculated Ensifer meliloti cultures following standard protocols (Simonsen and Stinchcombe, 2014a; Batstone et al., 2017). We cultured a sample of each rhizobial stock onto Petri dishes of tryptone yeast (TY) agar medium, then randomly sampled single colonies from the culture and re-plated them onto new plates. We repeated this process four times, taking a single pure colony from the culture each time. Once the rhizobial strains were purified, they were grown on liquid TY medium and diluted to final optical density at 600 nm (OD600) readings of 0.12 and 0.093 for Em1021 and Em1022, respectively. We combined both liquid cultures to create an inoculation broth with a 3:1 volume ratio of Em1021 and Em1022, and then inoculated each plant, excluding controls, with 1 ml of inoculation broth. We used three times as much Em1021 as Em1022 because the latter can easily outcompete the former when they are inoculated at the same density, or even at a 2:1 ratio (Batstone et al., 2017). After harvest, we counted the number of nodules present on plant roots as a proxy for investment in the mutualism by plant hosts and for rhizobial performance (Heath and Tiffin, 2009).

Statistical Analysis

We used R (version 3.6.2, R Core Team, 2019) to conduct analyses for the effects of light environment and plant genotype on the plant performance traits, leaf number, aboveground and belowground biomass, and on the mutualism-related trait, nodule number. We used a linear mixed model (LMM) with log-transformed response variables using the lmer function from the lme4 package (Bates et al., 2015). Untransformed data did not meet parametric tests for normality, homoscedasticity, and linearity (Wood et al., 2018), and model diagnostics for generalized LMMs showed a misfit with the data using the DHARMa package (version 0.2.7, Hartig, 2020). Our LMMs included light treatment as the fixed effect and random effects of genotype, block, and genotype-by-treatment and block-by-treatment interactions. We conducted significance tests using the Anova function from the car package (3rd edition, Fox and Weisberg, 2019) with fixed effects analyzed using the F-statistic and type III sum of squares. We analyzed the significance of random effects with log-likelihood ratio tests using the χ2 test statistic to compare models with and without the effect of interest. We halved the p values for log-likelihood ratio tests because they are one-tailed tests of whether a variance is greater than zero. For analysis of nodule number, we included total biomass as a fixed effect covariate to account for confounding effects of plant size on nodule production. We also assessed belowground biomass as a covariate for nodule number to account for the effects of root size. The results for nodule number were consistent across both these models, so we report and discuss the analysis using total biomass as a covariate.

We analyzed genetic correlations among traits and across treatments using the rcorr function from the Hmisc package (version 4.4-0, Harrell, 2020) to calculate a matrix of Pearson correlation coefficients and corresponding p values. For this analysis, we used line means of raw data to obtain a correlation matrix within the two light environments. Correlation plots were made using the corrplot function in the corrplot package (version 0.84, Wei and Simko, 2017).

G×E interactions can be due to shifts in genotype rank order or in the magnitude of genetic variance across environments; the latter can change the magnitude of the response to selection (Fisher, 1958), but only the former can effectively maintain genetic variation. To assess whether significant G×E interactions were due to changes in rank order versus changes in the magnitude of genetic variance, we used an equation originally described by Cockerham (1963) and applied by Batstone et al. (2020):where V is the genotypic variance component within environments i and j, r is the genetic correlation between environments i and j, and e is the number of environments.

The first half of (Equation 1) describes the G×E variance due to changes in genotype rank order, which can be isolated as:

The second half of (Equation 1) describes how changes in the magnitude of genetic variance contribute to G×E, which can be rewritten as:

We obtained Vrank and Vvariance using Equations 2 and 3, respectively, and estimated the relative portion of total G×E due to changes in rank (Vrank/Vrank + Vvariance), following Batstone et al. (2020). To obtain genotypic variance components (Vg), we used LMMs of log-transformed data within each treatment with genotype as the main random effect, as well as a random effect of block. We obtained Pearson correlation coefficients between environments (r) using the cor function from the stats package (version 3.6.2, R Core Team, 2019), with raw line means that were also log-transformed to remain consistent across all G×E analyses. We determined genetic correlations across environments for all traits (see Supplemental Figure 1), but we only used r from traits that showed significant G×E for the rank versus variance analysis.

Selection Analysis

To test whether there are any differences in how natural selection acts on mutualistic traits across environments (i.e., S×E), we conducted a selection analysis following the logic described by De Lisle and Svensson (2017) and Batstone et al. (2020). We used aboveground biomass as a proxy for plant fitness (most studies show a positive correlation between biomass and fitness components [Younginger et al., 2017]), and nodule number as our focal mutualistic trait. We also included leaf number and belowground biomass to analyze selection on plant performance traits, and to account for the high correlation of these traits with nodule number (see Results, Figure 2). We obtained raw line means to analyze ANCOVAs using linear models with the fitness proxy as a response variable, and focal traits and light treatment as covariate and interaction terms. We used three scaling methods for our S×E analysis: (1) absolute fitness, (2) relativizing fitness and standardizing traits on a global scale (across treatments), and (3) relativizing and standardizing on a local scale (within each treatment). We compared these methods to assess how differences in mean values and standard deviations across and within treatments might affect our S×E analysis (Batstone et al., 2020). We included an analysis using the absolute scale because it is free from any differences that might result from standardization and relativization in the other scaling methods, and it allows an assessment of the relationship between traits and fitness as they were measured. For the global analysis, we relativized aboveground biomass by dividing line means with a global mean across both treatments (line mean/global mean), and we standardized our focal traits by subtracting line means from the global mean and dividing this by a global SD ((line mean – global mean)/global SD). The local analysis was similar except we relativized aboveground biomass using local means from within each treatment (line mean/local mean), and standardized the focal traits using local means and SD ((line mean – local mean)/local SD). Any difference in selection for our focal traits between light environments (i.e., S×E) was indicated by a significant trait-by-treatment interaction. If a significant S×E interaction from the global analysis is not present in the local analysis, it indicates that the interaction is mainly driven by differences in mean fitness because those differences are mathematically eliminated in the local analysis. We used the global and local S×E analyses to determine whether estimates of selection gradients (β) were significantly different between environments (i.e., significant S×E interaction). However, we used β estimates from the local analysis to determine differences in the magnitude and direction of selection between environments, as these are most interpretable for making predictions of evolutionary change (De Lisle and Svensson, 2017).

Funding

We thank Canada for funding that supports our science, Megan Frederickson, Art Weis, and Helen Rodd for comments on our work, and Bill Cole and Thomas Gludovacz for horticultural assistance.

Author Contributions

Conceptualization, P.V. and J.R.S.; Methodology, P.V. and J.R.S.; Investigation, P.V.; Formal Analysis, P.V.; Writing – Original Draft, P.V.; Writing – Review, Editing, & Revising, P.V. and J.R.S.; Resources & Supervision, J.R.S.