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.

Raw Data

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, 22deg C day and 18deg 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 experimentalplants and 20 contamination control plants).

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 4deg C for counting and collecting rhizobial nodules. 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. We then dried and weighed the belowground portion to estimate belowground biomass.

Data 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 functionfrom the lme4 package (Bates et al., 2015). Untransformed data did not meet parametric tests for normality, homoscedasticity, and linearity ("checkAssumptions.R", 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 chi-squared test statistic to compare models with and without the effect of interest. We halved the pvalues for log-likelihood ratio tests because they are one-tailed tests ofwhether a variance is greater than zero. We analyzed genetic correlations among traits and across treatments us-ing the rcorr function from the Hmisc package (version 4.4-0,Harrell,2020) to calculate a matrix of Pearson correlation coefficients andcorresponding p values. To assess whether significant GxE interactions were due to changes in rank order versuschanges in the magnitude of genetic variance, we used an equation originally described by Cockerham (1963) and applied by Batstone et al.(2020). To obtain genotypic variance components (Vg), we used LMMs of log-transformed data withineach treatment with genotype as the main random effect, as well as arandom effect of block. We obtained Pearson correlation coefficients be-tween environments (rij) using the cor function from the stats package(version 3.6.2,R Core Team, 2019), with raw line means that were alsolog-transformed to remain consistent across all G3E analyses. We determined genetic correlations across environments for all traits (see Supplemental Figure 1), but we only usedrijfrom traits that showedsignificant G3E for the rank versus variance analysis.

Poor germination and establishment led to a final sample of 666 plants that survived until the end of the experiment. Only the surviving individuals are included in the dataset.