No disruption of rhizobial symbiosis during early stages of cowpea domestication

Modern agriculture intensely selects aboveground plant structures, while often neglecting belowground features, and evolutionary tradeoffs between these traits are predicted to disrupt host control over microbiota. Moreover, drift, inbreeding, and relaxed selection for symbiosis in crops might degrade plant mechanisms that support beneficial microbes. We studied the impact of domestication on the nitrogen fixing symbiosis between cowpea and root-nodulating Bradyrhizobium. We combined genome-wide analyses with a greenhouse inoculation study to investigate genomic diversity, heritability, and symbiosis trait variation among wild and early-domesticated cowpea genotypes. Cowpeas experienced modest decreases in genome-wide diversity during early domestication. Nonetheless, domesticated cowpeas responded efficiently to variation in symbiotic effectiveness, by forming more root nodules with nitrogen-fixing rhizobia and sanctioning non-fixing strains. Domesticated populations invested a larger proportion of host tissues into root nodules than wild cowpeas. Unlike soybean and wheat, cowpea showed no compelling evidence for degradation of symbiosis during domestication. Domesticated cowpeas experienced a less severe bottleneck than these crops and the low nutrient conditions in Africa where cowpea landraces were developed likely favored plant genotypes that gain substantial benefits from symbiosis. Breeders have largely neglected symbiosis traits, but artificial selection for improved plant responses to microbiota could increase plant performance and sustainability.

Methods
Data Collection:

Inoculation Experiments — Seeds were surface sterilized in bleach (5% sodium hypochlorite), rinsed in sterile ddH2O, scarified, and planted in bleach-sterilized 1-gallon plastic pots containing an autoclave-sterilized 50:50 mix of silica sand and limestone flour silica sand, which contains negligible nutrients to support plant growth (Regus et al., 2015). Three seeds were planted per pot from 06/13/2018 to 06/15/2018. Greenhouse temperatures averaged 86 °F + 14 °F (standard error, SE) and relative humidity was 55% + 20%. On 6/21/2018 seedlings were thinned to one plant per pot to size match the remaining seedlings among plant lines. One day later rhizobial inoculation followed.

For the clonal strain experiment, Fix+ and Fix- strains were plated on a modified arabinose gluconate medium (MAG; Sachs et al. 2009) and a single colony per strain was spread onto 8-10 plates to generate dense lawns. After 7 days of growth the cells were washed from the plates into liquid MAG media and cell concentrations were quantified by colorimetry. Liquid cultures were centrifuged at ~750g, spent media was removed, and the cells were resuspended in sterile ddH2O at a concentration of 1 x 108 cells ml-1. Plants were inoculated with either 5 ml of the Fix+ or Fix- clonal Bradyrhizobium cells (single inoculation, 5 x 108 cells), 5 ml of a mixture comprising equal concentrations of both strains (co-inoculation, 2.5 x 108 cells of each strain), or 5 ml sterile ddH2O as a control.

To investigate variation in symbiosis traits with other symbionts within a microbial community we performed a soil inoculation experiment Field soil was sampled from the University of California Riverside Agricultural Experiment Station at four sites within a 5-acre field where diverse cowpeas are propagated (coordinates: 33.967, -117.339; Huynh et al. 2018). The field has a history of cultivating cowpea during odd-numbered years, starting 2003. Additionally, the field is intercropped with barley and occasionally with other legume crops such as soybean and pigeonpea.Soil was passed through a sterilized 2mm sieve (6L per site), and apportioned into aliquots of 400g. From each sample, 400mL of sterile water was added, the sieved soil was shaken vigorously, filtered twice through 8 layers of sterile cheesecloth, and the filtered supernatants were pooled into sterile flasks, which were allowed to settle overnight at room temperature to inoculate with microbes in water and avoid sediments (Unkovich and Pate 1998). The supernatant from each flask was divided into two equal portions, one of which was autoclaved and allowed to cool to serve as a negative control, while the other was reserved at room temperature and used for inoculation. Seedlings were inoculated with 10mL of each microbial inoculum (alive or dead) and each one was separately plated (100ul) in MAG and incubated at 29°C for eight days to confirm high densities of slow growing bacteria like Bradyrhizobium.

In both experiments, plants were fertilized weekly by applying 10 ml of Jensen’s solution with 1 g/L K15NO3 (2% 15N by weight), which includes all the necessary micronutrients (Somasegaran and Hoben 1985), and a minimal concentration of nitrogen to support cowpea growth. Plant genotypes and inoculation treatments were randomly arranged within blocks in the greenhouse with five plant replicates per inoculation treatment x plant genotype combination, except for controls that had 3 replicates. The clonal strain experiment had 360 plants, including 300 that were inoculated (20 lines x 3 inoculation treatments x 5 replicates) and 60 control plants (20 lines x 3 replicates). The soil inoculation experiment had 160 plants, including 100 that received the live inoculum (20 lines x 5 replicates) and 60 that received the autoclaved control (20 lines x 3 replicates).

Plant harvest and nodule culturing – All plants were harvested during the course of three weeks, from 7/30/2018 to 8/3/2018 and from 8/13/2018 to 8/23/2018. Plants were removed from pots, washed free of sand, and dissected into root, shoot, and nodule portions. Nodules were counted and photographed. Rhizobia were sub-cultured from nodules of co-inoculated plants to differentiate Fix+ and Fix- strains. Nodules were crushed and streaked on MAG and isolated colonies were subcultured on MAG with rifampicin (500 μg ml-1) and streptomycin (1000 μg ml-1), selecting for Fix+, and YM media, on which Fix- exhibit fast growth and slimy appearance. Five nodules each from three co-inoculated plants per genotype were randomly picked and assessed (~15 nodules per genotype, 268 total). From each nodule, ~50 colonies were counted to estimate the proportion of Fix+ to Fix- strains (11,586 colonies in total).

Leaf 15N ‘atom per cent difference’, a measure of nitrogen fixation (Regus et al. 2014), was estimated as the percentage of 15N atoms over total nitrogen in each sample (Unkovich et al. 2008). The δ15N of each sample was calculated by comparing 15N abundance expressed as parts per thousand relative to atmospheric N2, these values were used to compare among plants inoculated with Fix+ and Fix- strains following the formula:

 δ15N ‰ =  sample atom%15N - 0.36630.3663 x 1000

To calculate these values individual leaves of each plant were oven dried, powdered using steel bead beaters at 14,000 rpm, and 4 mg per plant was transferred into individual tin capsules, including four replicates per genotype for the Fix+, Fix- and two replicates for control inoculation treatments (178 samples total). Isotopic analyses were performed at the UC Davis Stable Isotope Facility.

Data Processing:

Genome-wide variation of Cowpea accessions – To examine genetic variation and admixture between wild and cultivated cowpea we performed a combined analysis of 380 landraces and 58 wild cowpea accessions reported in Huynh et al. (2013) using the 1536-SNP GoldenGate genotyping assay. Huynh et al. (2013) analyzed wild and domesticated genotypes separately, with a focus on geographic origin. To maintain consistency with Huynh et al. 2013, SNPs with a minimum allele frequency (MAF) < 0.05 and with a call rate < 0.90 were discarded, for a final filtered set of 920 SNPs. Genetic differentiation was evaluated using a principal component analysis (PCA) with the package adegenet (Jombart, 2008). Admixture and structure were examined using the R package LEA (Frichot et al. 2014; Frichot and François 2015). One to ten ancestral populations (i.e., entropy criterion; K = 1 to 10) were assumed using 100 repetitions. To test if patterns of genetic diversity differed among populations, a generalized mixed model analysis using SNP loci as our random factor was implemented (Costa et al. 2021; Theodorou et al. n.d.). The GLMM with a Beta distribution and a logit link function was modeled using the package glmmTMB (Brooks et al. 2017; Douma and Weedon 2019). Post-hoc comparisons based on the model were performed with the R package emmeans (Searle et al. 2012). Population statistics were estimated with the R package hierfstat (Goudet 2005). 

To have a more robust estimation of the genomic-level variation and relationships among the twenty focal cowpea lines, we further genotyped the wild accessions using the Illumina Cowpea iSelect Consortium Array, screening 51,128 SNPs across the cowpea genome. Domesticated accessions were previously genotyped with the same array (Muñoz-Amatriaín et al, 2017). SNPs with a MAF < 0.1 and with a call rate < 0.95 were discarded using the R package snpReady (Granato et al. 2018), for a final filtered set of 34,762 SNPs. Pairwise genetic distances were estimated with the R package adegenet (Jombart 2008) and neighbor-joining was used to reconstruct phylogenetic relationships. Branch support values were evaluated by a bootstrap analysis where SNPs were sampled with replacement 100 times using the phylo.boot function of the package ape (Paradis and Schliep 2018).

Trait data analysis – Size comparisons among wild and domesticated populations were performed by calculating scale free measurements to minimize effects of initial seedling size. Investment into symbiosis was calculated by dividing the dry nodule biomass of each plant over the total biomass. Host growth response was calculated by subtracting the mean biomass values (i.e., shoot, root, and nodules) of the control plants within a population from the inoculated plants belonging to the same group, dividing by the control value, and multiplying the quotient by 100 (Regus et al. 2015). Means per population were calculated for plants harvested during the same week to account for variation in days post inoculation.

Host Growth Response %=Total biomass inoculated planti-Mean biomass controlsjMean biomass controlsj  x 100

Where i indicates plant replicate and j indicates population mean value.

Dry nodule biomass values of co-inoculated plants (where a subset of nodules was used for subculturing) were inferred by generating a wet-to-dry nodule weight linear regression (per genotype). To test for post-infection sanctions, a binomial test was used to evaluate whether nodule occupancy of Fix+ deviated from the null expectation of 50% given that the strains were inoculated in equal proportions. Results were analyzed independently for each genotype tested.

Linear mixed models (LMMs) were used to analyze differences in symbiosis traits among the three populations defined by Huynh and colleagues (2013), i.e., Genepool-1, 2, and wild cowpeas (three-population analysis). However, because landraces of Genepool-1 and 2 are each most closely related to wild cowpeas from the same region (Huynh et al., 2013), we also analyzed comparisons that divided the wild cowpeas into southern Africa populations (PI632890, PI632876, PI632892, PI632891; i.e., Wild-1) and northern Africa populations (PI632882, PI632879, PI632880, PI632881; i.e., Wild-2, four-population analysis). Inoculation treatment and population were treated as fixed effects, cowpea genotype and genotype x treatment interactions were treated as random effects, and days post inoculation was used as a covariate. Response variables were transformed if necessary, to improve normality. Analyses were performed using The R project for Statistical Computing version 3.6.1 (R Core Team 2020).

Components of trait variation – Independent linear mixed models were constructed to estimate the components of variation in each symbiosis trait under the clonal inoculation treatments, where genotypic effects could be best isolated. Models of variance-covariance structure were used to test whether the expression of additive genetic variance (s2a) in each symbiosis trait varied among treatments, or among the wild and domesticated populations (three-population analysis), and if the expression of s2a in populations varied among treatments. Because of limited sampling of plant genotypes, it was not practical to conduct this specific analysis using the four-population approach. The variance covariance matrix for the genotype effect known as the additive relationship matrix was estimated from the SNP data with the A.mat function in sommer (Covarrubias-Pazaran, 2016). To test if the additive genetic variance in the trait of interest varies among the levels of the factor of interest (treatment, population, population x treatment), a model where the among-genotype variance was constrained to be the same across levels was compared with a heterogeneous variance structure model (Table S1). Differences in the expression of genetic variance were assessed using log-likelihood tests among models (Shaw, 1991). Breeding values of each genotype were estimated by best linear unbiased prediction (BLUPs) (Bauer et al. 2006; Liu et al. 2008; Piepho et al. 2008), taking into account the additive relationship matrix among genotypes (genomic BLUPS, or GBLUPs). Narrow-sense heritability (h2) was estimated as the proportion of additive variance of two alleles at a locus over the phenotypic variance (h2= VA/VP), (Bernardo 2020). Analyses were performed in the R package sommer (Covarrubias-Pazaran 2016).

Genetic correlations among traits were estimated following Falconer (1952) and implemented by Etterson (2004) and Saxton (2004), where the correlation between any pair of traits i and j, rAij , was estimated as follows, where COVAij is the covariance between an individual’s breeding value for one trait and its breeding value for the other trait:

rAij= COVAijVAiVAj 

VAi is the genetic variance of trait i and VAj is the genetic variance of trait j. To estimate the genetic correlation between traits we performed multi-trait and multi-environment linear mixed models (Covarrubias-Pazaran 2016) with treatment, population, and days since inoculation as fixed factors, and cowpea genotype as random effect.