Association Study of Symbiotic Genes in Pea (Pisum sativum L.) Cultivars Grown in Symbiotic Conditions

Abstract

In garden pea (Pisum sativum L.), several symbiotic genes are known to control the development of mutualistic symbioses with nodule bacteria (NB) and arbuscular mycorrhizal fungi (AMF). Here, we studied whether the allelic state of the symbiotic genes was associated with the growth parameters of pea plants under single inoculation with NB and under double inoculation with NB + AMF. Using different statistical methods, we analyzed the dataset obtained from a pot experiment that involved 99 pea cultivars, 10 of which were characterized as having shortened internodes due to the presence of the natural mutation p.A229T in the developmental gene Le. The plant’s habitus strongly influenced most of the studied growth and yield parameters and the effectiveness of the symbiotic interactions under NB and NB + AMF inoculation. Double inoculation had different effects on Le⁺ (normal) and le⁻ (dwarf) plants with regard to nitrogen and phosphorus content in seeds. Regardless of the Le-status of plants, allelic states of the symbiotic gene LykX encoding the putative receptor of Nod factors (bacterial signal molecules) were shown to be associated with seed number, thousand-seed weight, and pod number at the level of FDR < 0.001, whereas associations of allelic states of the other studied symbiotic genes were less significant.

1. Introduction

The legumes (family Fabaceae) form a peculiar group among plants because they are able to actively interact with a wide range of beneficial soil microorganisms (BSMs), developing three mutually beneficial symbioses: (i) arbuscular mycorrhiza (AM) with fungi of phylum Glomeromycota, (ii) nitrogen-fixing root nodules with rhizobia, so-called legume–rhizobial (LR) symbiosis, and (iii) symbiosis with plant growth-promoting bacteria (PGPB). These symbioses are beneficial for both the host plant and its environment; therefore, legumes are well-suited for cultivation according to the modern concept of sustainable agriculture [1,2]. For successful implementation of this idea, new varieties of legumes that are adapted for effective interaction with microorganisms need to be bred [3,4]. One approach for achieving this goal is the mobilization of plant genetic resources, i.e., searching for valuable alleles of genes that can improve the symbiotic properties and plant growth parameters when introduced in the genotype.

In several species of legumes, numerous symbiotic (Sym) genes were described as necessary for symbioses’ formation and functioning. Using approaches of forward and reverse genetics, over 50 genes were characterized with respect to their role in LR symbiosis and AM development in model legumes Medicago truncatula Gaertn. and Lotus japonicus (Regel.) K. Larsen [5] and crop cultures such as pea (Pisum sativum L.) [6], soybean (Glycine max (L.) Merr.), and common bean (Phaseolus vulgaris L.) [5]. Sym genes are possible candidates for determining the agriculturally important traits in plants grown in symbiotic conditions; indeed, the genomic positions of some Sym genes are similar to the locations of quantitative trait loci (QTLs) related to fruit, root, or nodule traits in M. truncatula [7], Ph. vulgaris L. [8], and L. japonicus [9]. Colocalization with QTLs of seed yield and seed protein content was also demonstrated for some developmental genes of pea, which featured natural mutations that influenced the plant’s habitus (le⁻, causing shortened internodes, and af, causing conversion of leaves into tendrils) [10,11,12,13]. However, no direct role of the genes Le or Af in the development and functioning of either LR symbiosis or AM has been recorded to date.

At present, legume plant breeding is mainly focused on improving the performance of LR symbiosis [3,14,15,16]. The beneficial effect on AM colonization on legumes has also been documented, although it was hardly distinguished from the effect of LR symbiosis [17]; therefore, breeding for improving the effects of AM itself is problematic. For this reason, and because in field conditions plants interact with the indigenous soil microbial community, it was proposed to breed legume crops considering the integral trait called “effectiveness of interactions with beneficial soil microorganisms” (EIBSM) [3,4]. EIBSM is a measure of increase in seed biomass (and/or other agronomically important parameters such as plant biomass, seed protein content, etc.) of plants inoculated with beneficial microorganisms (in pot or field experiments) as compared to uninoculated control [4]. The EIBSM trait was analyzed in 3-year field trials in garden pea, and a total of 26 genotypes were rated according to their ability to increase biomass and yield under complex inoculation with rhizobia (Rh) and AM fungi; thus, high- and low-EIBSM genotypes were identified [3,18]. The molecular mechanisms underlying the manifestation of the EIBSM trait in pea, as well as in other legumes, remain generally undiscovered; however, in our recent studies, we found that high-EIBSM genotypes (also referred to as highly responsive to inoculation with BSM), in contrast to low-EIBSM ones (low responsive), tended to prolong the seed-filling stage under double inoculation with Rh and AM and demonstrated stricter control over microorganisms in roots [19,20].

Over 20 years ago, a pot experiment was conducted at ARRIAM (St. Petersburg, Russia) [21] in which several growth parameters of 99 pea cultivars were evaluated in conditions of (i) single inoculation with rhizobia and (ii) double inoculation with Rh + AM. However, a thorough analysis of the obtained dataset was not performed; thus, we aimed to conduct a deeper analysis of the earlier data with the use of appropriate statistical methods. Also, in our new experiment, we tested whether the allelic state of known Sym genes influenced the growth parameters of pea plants under single and double inoculation and whether the particular allelic variants were associated with increases in any of the investigated parameters under double inoculation as compared to a single inoculation. We found an impact of plant habitus (related to the allelic state of the Le gene) on values of most parameters, including AM-treatment-caused increments. Further, we identified nonsynonymous SNVs (single nucleotide variants) is some Sym genes with statistically significant association with several parameters and with the AM-treatment-caused increment in the thousand-seed weight (TSW).

2. Materials and Methods

2.1. Plant Material

Pea (Pisum sativum L.) cultivars were obtained from the Pea World Germplasm Collection of the N.I. Vavilov All-Russian Institute of Plant Genetic Resources (VIR, St. Petersburg, Russia). The list of accessions with their morphological characteristics is presented in Supplementary Table S1. The low (stem length about 30 cm) determinate pea (P. sativum) cv. Finale (Cebeco, Rotterdam, The Netherlands), having stable yields and wide adaptation [22], was used as a reference genotype for AM development.

2.2. Microorganisms

The fungal isolate Rhizophagus irregularis BEG144 (Syn. Glomus intraradices strain 8) used in the study was previously characterized as forming highly effective AM symbioses with many agricultural crops [23]. The isolate was initially provided by the International Bank for the Glomeromycota (Dijon, France) and was maintained in alternating Plecthrantus australis and Sorghum spp. pot cultures.

The nodule bacteria strain Rhizobium leguminosarum bv. viciae RCAM 1026 [24], characterized by high efficiency and competitiveness, was used for inoculation in the vegetation experiment.

2.3. Vegetation Experiment for Assesment of the Symbiotic Effectiveness of 99 Pea Cultivars

The experiment conducted by Yakobi et al. [21] was carried out in a greenhouse at ARRIAM (Pushkin, St. Petersburg, Russia) in a randomized complete block design in three replicates (five plants per pot for each treatment variant (Rh or Rh + AM)). Humus horizon of sod podzolic light-loamy fallow soil as growth substrate: pH (KCl) = 4.9, organic matter 3%, available phosphorus 1.26 mg/100 g of soil (extraction with 0.2 N HCl). Air-dried soil was sterilized by autoclaving (1 h, 2 atm), supplemented with CaCO₃ (pur. ≥ 98.5%, Lenreactiv, St. Petersburg, Russia) (1.65 g per 1 kg of soil), and mixed with sterile quarry quartz sand (Fraction 0–0.63 mm, SiO₂ ≥ 80%, Russia) (1:1 v/v). No fertilizers were added to the substrate.

The isolate BEG144 of Rhizophagus irregularis was propagated on the roots of Sorghum vulgare Pers. Root–soil mixture after vegetation of S. vulgare with or without AM (15 g per each pea plant) was used for inoculation as AM+ and AM− variants. All plants were also inoculated with the R. leguminosarum strain RCAM 1026 in amounts of 10⁷ CFU (colony-forming units) per plant.

Plants were grown in pots with 4.2 kg of the substrate in three replicates for each treatment. The pots were arranged in three blocks such that each treatment was replicated once in each block and the cultivars were randomly located within a block. Pots in which some plants prematurely died were excluded from the analysis, so that no fewer than 10 plants per variant of treatment were analyzed.

Plants were harvested for the analysis of growth parameters at the stage of complete seed maturation (approximately 3 months post inoculation). Nitrogen content in seeds was determined using the Kjeldahl method, and phosphorus content in seeds was determined photocolorimetrically [21].

2.4. SNVs Detection and Cultivar Screening

For SNP analysis in our new experiment, the genomic DNA was extracted from one-week old pea seedlings (3–5 seeds of each cultivar). Fresh plant tissue was disrupted between glass beads using FastPrep-24 grinder (MP Biomedicals, Irvine, CA, USA). This was followed by lysis and DNA extraction with in-house made CTAB buffer: CTAB 2% (w/v) (pur., Sisco Research Laboratories, Mumbai, India), NaCl 1.4 M (puriss., Lenreactiv), Tris Cl (pH 8.0) 0.1 M (pur., Acros Organics, Geel, Belgium), EDTA 20 mM (puriss. spec., Lenreactiv). The supernatant was twice purified with chloroform (pur., Lenreactiv), and DNA was precipitated with fresh, distilled, ice-cold ethanol (MedChimProm, Balashikha, Russia) and dissolved in Milli-Q water.

In order to test whether the allelic state of symbiotic genes affected the properties of pea plants grown in symbiotic conditions, the set of eight genes that act during the development of root nodules and AM was selected (see Table 1).

Using the extracted total DNA, for cultivars k-925, k-1693, k-3064, k-3358, k-7128, and k-8274, the fragments of the genes of interest were PCR amplified using a T100 PCR Thermal Cycler (Bio-Rad Laboratories, San Francisco, CA, USA) for subsequent sequencing. The primers for amplification were designed based on transcript sequences presented in available pea transcriptome assemblies and sequence databases in accordance with the putative exon–intron structures of the genes constructed based on the structures of homologous genes of M. truncatula (see Supplementary Materials Gene structure schemes). For amplification, we used a Tersus polymerase PCR kit (Evrogen, Moscow, Russia). The optimal annealing temperature for each primer pair was determined empirically in test amplifications. The sequences of primers and the used primer pairs and gene structures are presented in Supplementary Materials Gene structure schemes. The obtained amplicons were sequenced on an ABI Prism 3500XL (Thermo Fisher Scientific, Waltham, MA, USA) with the usage of the original sequencing kit supplied by the manufacturer of the sequencer.

The genomic sequences were assembled into contigs and truncated to CDSs (see Supplementary Sequences). From CDSs translated into peptides, the nonsynonymous SNVs were identified. For the detection of the missense variants in cultivars, we developed CAPS and dCAPS (cleaved amplified polymorphic sequence/derived CAPS [25]) markers suitable for DNA genotyping (Supplementary Table S2).

The target DNA fragments were PCR amplified on the matrix of total DNA using Taq-polymerase-based ScreenMix (Evrogen) individually for each cultivar and digested by corresponding restriction endonucleases in the conditions recommended by the manufacturer. The analysis of the restriction patterns was performed either by agarose gel electrophoresis (the mode was selected individually for each marker) or with use of the MultiNA microchip electrophoresis system (Shimadzu, Kyoto, Japan) using the original DNA-500 or DNA-1000 kits (Shimadzu).

In total, 93 cultivars were screened using the created markers, and the allelic state of each marker was determined. Along with Sym-genes, the gene Le, which encodes gibberellin 3-beta-dioxygenase 1 and defines the stem length in pea [26], was included in the study. Earlier sequencing data of the first exons of LysM-RLK (LysM-containing receptor-like kinase) genes Sym37, K1, and LykX on the same set of 93 pea genotypes [27] were also added to the dataset.

Table 1. The Sym genes selected for association analysis.

Gene Identifier	Accession	Description
IGN1 [28]	KR047192	Ankyrin-repeat protein with transmembrane regions
	GICP01056763 1
SEN1 [29]	KY888171	Control of rhizobial differentiation
	GDTM01047803 2
SST1 [28]	KC008603	Symbiotic sulphate transporter
	GDTM01043951 2
Sym29	AJ495759	CLAVATA1-like receptor
	GICP01072631 1
STR1	GICP01133083 1	ABC transporter G family member STR, required for arbuscule development
STR2	GICP01179536 1	ABC transporter G family member STR, required for arbuscule development
AMT2;3	PsCam053936 3	Ammonium transporter
AMT2;5	GICP01050721 1	Ammonium transporter

¹—accession from NCBI Nucleotide Database from pea transcriptome shotgun assembly GICP00000000 [30]; ²—accession from NCBI Nucleotide Database from pea transcriptome shotgun assembly GDTM00000000 [31]; ³—accession from “Pisum sativum Gene Atlas” [32].

Table 1. The Sym genes selected for association analysis.

Gene Identifier	Accession	Description
IGN1 [28]	KR047192	Ankyrin-repeat protein with transmembrane regions
	GICP01056763 1
SEN1 [29]	KY888171	Control of rhizobial differentiation
	GDTM01047803 2
SST1 [28]	KC008603	Symbiotic sulphate transporter
	GDTM01043951 2
Sym29	AJ495759	CLAVATA1-like receptor
	GICP01072631 1
STR1	GICP01133083 1	ABC transporter G family member STR, required for arbuscule development
STR2	GICP01179536 1	ABC transporter G family member STR, required for arbuscule development
AMT2;3	PsCam053936 3	Ammonium transporter
AMT2;5	GICP01050721 1	Ammonium transporter

¹—accession from NCBI Nucleotide Database from pea transcriptome shotgun assembly GICP00000000 [30]; ²—accession from NCBI Nucleotide Database from pea transcriptome shotgun assembly GDTM00000000 [31]; ³—accession from “Pisum sativum Gene Atlas” [32].

2.5. Statistical Analysis

Statistical analysis included correlation analysis, t-tests, F-tests, ANOVA, and regression analysis. For multiple comparisons, FDR correction was applied. In regression analysis, the LRT (likelihood ratio test) was used for variable significance testing. All statistical analyses were performed using R Statistical Software (version 3.5.3).

2.6. Analysis of AM Development

In order to analyze AM development, a new experiment was carried out in nurse pots. Nurse pots of chives (Allium schoenoprasum L.) with established AM were prepared (see [33] for details), which provided an efficient inoculation procedure for P. sativum plants, including synchronization of root colonization by R. irregularis.

Pea seeds were surface disinfected as follows: 1 min in 96% ethanol (MedChimProm), a rinse with sterile water, 8 min in a 5% NaClO aqueous solution (“Belizna”, Kaustik, Volgograd, Russia), and a thorough rinse with sterile water. Disinfected pea seeds were germinated on sterile humid vermiculite in Petri dishes for 3 days at 27 °C in the dark. Three seedlings of equal size were planted into each nurse pot around a 12-week-old mycorrhizal chive plant (3 pots per genotype). Three hundred milliliter ceramic flowerpots, with silica-rich marl mineral substrate (“Barsik”, NPO Novye Tehnologii Ltd., Krasnodar, Russia) and supplemented with 1 g/L calcium orthophosphate (pur., Lenreactiv) were used. Pots with substrate were sterilized by autoclaving for 60 min at 134 °C and 0.22 MPa. Plants were grown in a growth chamber (model VB 1514, Vötsch Industrietechnik GmbH, Balingen-Frommern, Germany) under the following conditions: day/night, 16/8 h; temperature, 24/22 °C; relative humidity, 75%. Plants were supplemented once a week with ½-strength Hoagland’s solution [34] without phosphate (50 mL per pot) and watered as needed. The experiment was carried out in two replicates.

Pea plants were harvested at 21 days of cocultivation with A. schoenoprasum in the nurse pots. Root systems were washed in cold tap water, and several randomly selected lateral roots with total lengths of approximately 50–100 cm were collected individually from each system. Root samples were cleared with 10% KOH and stained with black ink (Sheaffer Manufacturing Co., LLC, Ft. Madison, IA, USA) according to Vierheilig et al. [35]. After staining, the roots were washed once with distilled water and covered in glycerol (pur., Lenreactiv); 30 cm of root fragments for each plant was laid out on a glass slide, covered with a second slide, and squashed. AM development was examined using a light microscope (Axiovert 35, Zeiss, Jena, Germany). AM development was quantitatively assessed according to Trouvelot et al. [36]: F% = frequency of mycorrhizal colonization (which reflects the proportion of mycorrhized root fragments among all those analyzed), M% = intensity of mycorrhizal colonization (which reflects the proportion of the root length colonized by the fungus), a% = arbuscule abundance in mycorrhizal root fragments (which characterizes the functional state of the fungus), and A% = saturation of the root system with arbuscules (A% = M × a/100). For statistical analysis, the parameters were subjected to arcsine transformation to normalize data.

3. Results

3.1. General Description of the Dataset

The dataset contains the mean values of eight traits measured for 99 cultivars (see Supplementary Table S1) grown in pots with presterilized soil under two conditions: inoculation with rhizobia (hereinafter Rh) and inoculation with rhizobia and arbuscular mycorrhizal fungi (hereinafter Rh + AM).

The effect of mycorrhization on the analyzed parameters was well pronounced (Figure 1;

Table 2. The traits analyzed for pea plants grown in Rh and Rh + AM conditions.

Measured Trait	Mean Value, Rh	Mean Value, Rh + AM	Mean Increase (abs) 1	Mean Increase (rel), %
Shoot weight, g	1.96, SD = 0.668	3.62, SD = 0.725	1.67 ± 0.1	98.90%
Shoot length, cm	65.5, SD = 21.3	79.3, SD = 22.3	13.8 ± 1.9	24.00%
Pod number	1.56, SD = 0.43	2.61, SD = 0.622	1.05 ± 0.097	73.10%
Seed number per plant	4.87, SD = 1.97	10.6, SD = 3.13	5.75 ± 0.43	133.40%
Seed weight per plant, g	0.74, SD = 0.267	1.89, SD = 0.364	1.15 ± 0.061	183.30%
Thousand-seed weight, g	160, SD = 49.0	190, SD = 55.0	30.2 ± 3.7	20.30%
Seed nitrogen content, mg/g	38, SD = 5.25	30.5, SD = 3.17	−7.54 ± 1.2	−18.40%
Seed phosphorus content, mg/g	0.816, SD = 0.116	1.27, SD = 0.283	0.453 ± 0.06	58.30%

¹ value after ± denotes spread of 95% confidence interval.

); the paired samples t-test persuasively (p-value < 1.17 × 10⁻²¹) showed differences for all parameters (see

Table 3. Comparisons between the traits values in Rh and Rh + AM conditions.

Measured Trait	Paired t-Test 1 t-Statistic, p-Value	F-Test 2 F-Statistic, p-Value	Pearson Correlation Coefficient (R)	Coefficient of Determination, Adjusted (R2adj)
PW	31.64, p-value = 3.07 × 10−53	1.18, p-value = 0.425	0.720, p-value = 4.55 × 10−17	0.513
PL	14.71, p-value = 1.53 × 10−26	1.09, p-value = 0.659	0.910, p-value = 8.37 × 10−39	0.826
PN	21.49, p-value = 7.21 × 10−39	2.09, p-value = 3.29 × 10−4	0.625, p-value = 4.93 × 10−12	0.384
SN	26.47, p-value = 1.95 × 10−46	2.53, p-value = 6.53 × 10−6	0.731, p-value = 8.73 × 10−18	0.53
SW	37.32, p-value = 9.76 × 10−60	1.86, p-value = 2.45 × 10−3	0.560, p-value = 1.62 × 10−9	0.307
TSW	16.03, p-value = 3.89 × 10−29	1.26, p--value = 0.253	0.941, p-value = 1.40 × 10−47	0.885
NtS	−12.34, p-value = 1.17 × 10−21	0.37, p-value = 1.13 × 10−6	0.019, p-value = 0.85	−0.01
PhS	15.05, p-value = 3.23 × 10−27	5.96, p-value < 2.2 × 10−16	0.052, p-value = 0.611	−0.008

¹ Two-tailed paired samples t-test with unequal variances, df = 98; ² two-tailed F-test of equality of variances, df = 98/98. PW—plant shoot weight; PL—plant shoot length; PN—pod number; SN—seed number per plant; SW—seed weight per plant; TSW—thousand-seed weight; NtS—seed nitrogen content; PhS—seed phosphorus content.

). The increments were positive for all parameters except NtS. The relative increase due to mycorrhization varied from 20% for TSW to 183% for SW, and for NtS, a relative decrease of −18% was observed (Table 2).

3.2. Correlation Analysis

As was previously seen [21], a significant correlation between values in Rh and Rh + AM treatments was found for most traits (see Table 3) except for NtS (R = 0.019, p-value = 0.85) and PhS (R = 0.052, p-value = 0.611). According to the F-test, for the plant growth traits PW and PL and for the TSW, the sample variance was similar between Rh and Rh + AM treatments (Table 3), which indicates that the addition of AM evenly influenced these traits in all cultivars (i.e., mycorrhization shifted only the mean values while the variance remained the same) (Figure 1). In the case of the traits PN, SN, and SW, mycorrhization significantly affected the variance (Table 3), indicating that the variability of the tested pea cultivars by these parameters was high and could be a subject of breeding work.

Similarly, the coefficient of determination (R²), which reflects the proportion of the variance that can be explained by the individual genetic features of the cultivars (regardless of the treatment), was shown to be high for the parameters PL and TSW (82.6% and 88.5%, respectively). Indeed, AM treatment increased the values of the parameters by approximately 20% only (Table 2). For the traits PN, SN, and SW, the coefficient of determination was much lower (Table 3), indicating that the AM treatment influenced these parameters differently in different cultivars.

To evaluate the correlations between traits within experimental treatments, we performed (1) Pearson correlation analysis under both experimental conditions (Rh and Rh + AM inoculation) and (2) regression analysis with a model combining both conditions as a factor. The scatter plots, correlation coefficients (R), and coefficients of determination (R2) are presented in Figure 2. The values of PW, PL, PN, SN, and SW were significantly correlated with one another under both experimental conditions. The SW trait demonstrated the strongest correlation with other mentioned traits, PW, PL, PN and SN, with a level of determination reaching 96% for PW. Further, this group of parameters showed an interesting pattern of relationship with nitrogen content (NtS) and phosphorus content (PhS) in seeds: a negative correlation with NtS and no correlation with PhS under Rh inoculation (Figure 2), and, oppositely, no significant correlation with NtS (except weak negative correlation for PL) and negative correlation with PhS under double inoculation with Rh + AM. The TSW trait had the least relation with other parameters, which is in line with the observation that seed size is highly genetically determined and is faintly affected by mycorrhizal fungi. However, the values of TSW positively correlated with PL and PW and negatively correlated with PN and SN (Figure 2).

Alongside correlation between traits, we estimated correlations between absolute AM-caused trait increments (ΔT = TRh.AM − TRh) for all traits (see Figure 3). The increments were significantly positively correlated in the group of ΔPW–ΔPL–ΔPN–ΔSN–ΔSW, indicating that the values of these parameters increased together upon AM inoculation. Also, weak positive correlation was found for ΔSW–ΔTSW, and negative correlations were obtained for ΔSW–ΔPhS, ΔNtS–ΔPhS, and ΔNtS–ΔPL.

Some AM-caused trait increments (Δ) also correlated with values of corresponding traits in Rh conditions (Figure 4), and this relationship became clearer when the relative increments (α = ΔT/TRh) were considered (Supplementary Figure S1). The latter relations were not linear and fit the exponential approximation (except NtS, for which the best approximation was linear). Generally, the plants with lesser values of traits in the absence of AM were characterized by higher relative increments of these traits due to mycorrhization.

In the αPL–PL-Rh plot, we noticed a separation of dots into two clusters corresponding to short and tall plants, and the shorter plants were characterized by wider spread of αPL. In pea, one of the genes that determine the stem length is Le, which encodes gibberellin 3-beta-dioxygenase [26]. The natural mutation le (p.A229T) causes shortening of the internodes, while an intact copy of the gene determines the normal length of internodes and a whole plant. We genotyped cultivars by Le allele (Supplementary Table S3, see below) and used its allelic state as a factor in regression models connecting traits’ absolute and relative increments with trait values. The results showed that the Le/le allele status affected absolute increments of SW, PL, TSW, NtS, and PhS (p-value < 0.05) and relative increments of PW, PL, PN, TSW, NtS, and PhS (p-value < 0.05) (Supplementary Figure S1), which implies that plant length might influence pea responsivity to AM inoculation. However, it is worth noting that le⁻ cultivars accounted for only approximately 10% of the analyzed pea genotypes. Therefore, further analysis of a set of pea cultivars with equal ratio of Le⁺–le⁻ plants is needed to support the present result.

3.3. Principal Component (PC) Decomposition

Because of essential relations and correlations between traits, to study symbiotic responsivity, it seems rational to consider not a trait but a characteristic encompassing various measurable parameters and reflecting the interrelationship among these parameters. Here, we used linear principal component (PC) decomposition. The PCs were calculated based on centered, SD-normalized values of all eight traits from both Rh and Rh + AM conditions. The trait vectors in PC coordinates are presented in Figure 5. The first four PCs explained 94.1% of total variance. As expected, the correlation analysis among the PCs showed no significant relationship (Supplementary Figure S2), but correlation between Rh and Rh + AM condition persisted for all PCs (Supplementary Figure S3). The increments in PC coordinates depended on the corresponding PC values of the Rh condition (Supplementary Figure S4) for all PCs except PC3.

3.4. Genotyping of Cultivars

In order to test whether the allelic state of symbiotic genes affected the properties of pea plants grown in symbiotic conditions, we sequenced a set of eight genes in six pea cultivars (k-925, k-1693, k-3064, k-3358, k-7128 and k-8274) and identified the missense SNVs (

Table 4. The selected Sym genes and detected missense variations in putative protein sequences.

Gene Identifier	Protein Substitutions
IGN1	D108E
SEN1	P22A
SST1	S545Y
Sym29	F16S
	P43H
	A53G
	I171F
	S618P
	T623I
STR1	E109K
	P397S
	A401P
	H475Y
STR2	T676I
AMT2;3	A21P *
	N215I
	Y385N
AMT2;5	—

*—not included in analysis as a rare variant.

). We screened the 93 pea cultivars for these allelic variants using created CAPS/dCAPS markers (Supplementary Table S2); the complete genotyping results are presented in Supplementary Table S3.

3.5. Association Analysis between the Allelic State of Sym-Genes and Traits

After filtering out the rare variants (present in less than ten cultivars), 52 nonsynonymic SNVs (single nucleotide variants) were analyzed (Supplementary Table S3).

To estimate the possible influence of the allelic states of selected genes on plant traits, we performed a series of two-way ANOVA tests with inoculation (Rh or Rh + AM) as the first factor and the allelic state of an SNV as the second factor. The full results are presented in Supplementary Table S4. The tests confirmed the effect of AM inoculation on all parameters; potential association of some allelic variants with some traits was also found (

Table 5. Associations between SNVs and traits.

SNV	Trait	p-Value	FDR	SNV	Trait	p-Value	FDR
Le.A229T	PL	4.00 × 10−38	3.20 × 10−37	SEN1.P22A	PL	1.24 × 10−3	9.92 × 10−3
Le.A229T	PW	9.54 × 10−9	3.81 × 10−8	LykX.L42K	TSW	1.53 × 10−3	8.65 × 10−3
Le.A229T	NtS	7.82 × 10−7	2.08 × 10−6	K1.E90D	TSW	1.86 × 10−3	0.0149
LykX.L13V	SN	4.43 × 10−6	3.55 × 10−5	LykX.L42K	SN	2.16 × 10−3	8.60 × 10−3
Le.A229T	PhS	8.36 × 10−6	1.67 × 10−5	Sym37.L56V	SN	2.21 × 10−3	8.80 × 10−3
LykX.L9F	TSW	1.20 × 10−5	9.57 × 10−5	IGN1.D108E	PN	2.36 × 10−3	9.20 × 10−3
LykX.A142I	TSW	2.35 × 10−5	1.88 × 10−4	IGN1.D108E	PL	3.47 × 10−3	9.20 × 10−3
Le.A229T	SW	3.37 × 10−5	5.39 × 10−5	LykX.A142I	SN	3.58 × 10−3	0.0143
LykX.L9F	SN	4.96 × 10−5	1.98 × 10−4	LykX.S134G	SW	5.82 × 10−3	0.0466
LykX.L13V	PN	1.30 × 10−4	5.19 × 10−4	AMT2;3.N215I	SW	6.46 × 10−3	0.0204
LykX.V86I	SN	1.97 × 10−4	1.02 × 10−3	AMT2;3.N215I	PL	7.63 × 10−3	0.0204
LykX.F184S	SN	1.97 × 10−4	1.02 × 10−3	IGN1.D108E	TSW	7.68 × 10−3	0.0154
IGN1.D108E	SN	2.09 × 10−4	1.67 × 10−3	Sym37.L56V	PN	8.71 × 10−3	0.0232
Sym37.S45L	TSW	2.30 × 10−4	1.84 × 10−3	K1.E90D	SN	0.0103	0.0411
LykX.V86I	TSW	2.54 × 10−4	1.02 × 10−3	AMT2;3.N215I	TSW	0.0109	0.0218
LykX.F184S	TSW	2.54 × 10−4	1.02 × 10−3	LykX.L9F	PN	0.0119	0.0316
Le.A229T	SN	4.29 × 10−4	5.72 × 10−3	LykX.L13V	PL	0.0136	0.0273
Sym37.L56V	TSW	5.40 × 10−4	4.32 × 10−3	Sym37.S45L	PN	0.0168	0.0480
Le.A229T	PN	6.08 × 10−4	6.95 × 10−4	Sym37.S45L	SN	0.0180	0.0480
LykX.L13V	TSW	7.74 × 10−4	2.06 × 10−3	LykX.L13V	SW	0.0296	0.0461
AMT2;3.N215I	PW	1.17 × 10−3	9.36 × 10−3	LykX.L13V	NtS	0.0346	0.0461

). All traits except TSW showed a relationship with Le.A229T SNV, which defines the mutant short internode phenotypes. For the traits NtS and PhS, the interaction between the factors LE and inoculation was noticed, the result implying that inoculation may have had different effects on element content in seeds of plants differing in Le:A229T. In addition to Le.A229T, allelic states of LykX, IGN1, K1, Sym37, AMT2;3, and SEN1 were shown to be associated with different SN, TSW, PL, PN, and PW traits (FDR < 0.05, Table 5).

The effect of Le.A229T was also clearly visible in PC coordinates, where it appeared to be associated with PC1, PC3, and PC4 (

Table 6. Associations between SNVs and PCs.

SNV	PC	p-Value	FDR	SNV	PC	p-Value	FDR
Le.A229T	PC3	1.23 × 10−25	4.93 × 10−25	LykX.L42K	PC2	2.40 × 10−3	9.61 × 10−3
Le.A229T	PC1	2.04 × 10−13	4.08 × 10−13	IGN1.D108E	PC2	3.66 × 10−3	0.0130
Le.A229T	PC4	1.66 × 10−8	2.22 × 10−8	K1.E90D	PC2	4.98 × 10−3	0.0199
LykX.L9F	PC2	1.13 × 10−5	4.51 × 10−5	AMT2;3.N215I	PC1	5.13 × 10−3	0.0205
Sym37.S45L	PC2	9.71 × 10−5	3.88 × 10−4	IGN1.D108E	PC1	6.51 × 10−3	0.0130
LykX.A142I	PC2	2.02 × 10−4	8.06 × 10−4	STR2.T676I	PC2	7.70 × 10−3	0.0308
LykX.L13V	PC2	2.21 × 10−4	8.85 × 10−4	LykX.S134G	PC1	0.0112	0.0450
Sym37.L56V	PC2	2.48 × 10−4	9.90 × 10−4	SEN1.P22A	PC4	0.0145	0.0484
LykX.F184S	PC2	2.85 × 10−4	1.14 × 10−3	K1.E90D	PC4	0.0219	0.0437
LykX.V86I	PC2	2.85 × 10−4	1.14 × 10−3	SEN1.P22A	PC3	0.0242	0.0484
LykX.L13V	PC1	4.43 × 10−4	8.87 × 10−4	LykX.L13V	PC3	0.0353	0.0471

and Supplementary Table S5). SNVs in symbiotic genes LykX, Sym37, K1, and IGN1 were associated with PC2, and LykX.13LV was associated with with PC1 (FDR < 0.05). Obviously, the PC2 axis described the variation in the TSW trait (see Figure 5a), and in accordance with this, most of these SNVs were also found to be associated with the TSW trait.

3.6. Regression Analysis of the Association of the AM-Caused Increments with SNVs

As was shown above, the AM-caused increments were heavily dependent on the absolute values of the corresponding parameters in Rh conditions, which was true both for the traits (see Figure 4 and Supplementary Figure S1) and for the PCs calculated based on these traits (Supplementary Figure S4). The AM-caused increments of PL, PN, TSW, NtS, and PhS were also found to be dependent on the Le/le allelic state (Figure 4) (however, this effect was not visible in PC coordinates). Based on these observations, for the analysis of the possible impact of SNVs on AM-caused increments, we used generalized linear regression models, which optionally included Rh-value and the Le allele as covariates (

Table 7. The models chosen for regression analysis of AM-caused increments in traits.

Trait	The Formula Used in the Model
PW	dPW~PW.Rh + LE
PL	dPL~PL.Rh + LE
PN	dPN~LE
SN	dSN~1
SW	dSW~SW.Rh
TSW	dTSW~LE
NtS	dNtS~NtS.Rh + LE
PhS	dPhS~PhS.Rh + LE

). For trait increments, the covariates included in the models were defined based on likelihood ratio tests (LRTs). The final chosen models are presented in Table 7. For PCs, only Rh value was included in the model as a covariate (dPC_i~PC_i^Rh). These models were used to test the significance of SNV impact on traits and PCs.

No strong associations were found for any of the SNVs with the studied traits; however, weak associations of SNVs in genes AMT2;3, K1, and LykX with AM-caused increments of TSW were identified (

Table 8. Selected associations with trait increments.

SNV	Trait	p-Value	FDR
AMT2;3.215NI	dTSW	0.0050	0.0402
K1.31LI	dTSW	0.0047	0.0379
LykX.16VF	dTSW	0.0038	0.0301
LykX.184FS	dTSW	0.0056	0.0447

). The SNVs AMT2;3.215NI, LykX.16VF, and SST1.545SY were associated with dPC2 and/or dPC1 (

Table 9. Selected associations with PC increments.

SNV	PC	p-Value	FDR
AMT2;3.215NI	dPC2	0.0061	0.0243
LykX.16VF	dPC2	0.0093	0.0371
SST1.545SY	dPC1	0.0215	0.0430
SST1.545SY	dPC2	0.0164	0.0430

). Thus, in addition to Le, the genes AMT2;3, K1, LykX, and SST1 can be considered candidate genes that contribute to the determination of the responsivity of pea plants to AM + Rh inoculation as compared to single Rh inoculation.

3.7. Association of the Root Mycorrhization Level and the Stem Length in Pea Cultivars

Since we found that the stem length of pea plants (determined by the Le allele) was strongly associated with the mycorrhiza-derived increments of most of the growth parameters, and because previous studies did not focus on the development of AM [21], we tested whether the stem length was associated with the mycorrhization parameters. To do so, we determined the level of root colonization in six cultivars from the analyzed set (for which the sequencing of the symbiotic genes was performed), along with the cv. Finale, for which the AM development was earlier characterized in detail [33,37,38].

Four parameters of AM development were analyzed: the frequency of mycorrhizal colonization (F%), the intensity of mycorrhizal colonization (M%), the abundance of arbuscules in the mycorrhizal part of the root (a%), and the saturation of the root system with arbuscules (A%). The F% values did not differ in the studied genotypes and were at a very high level (about 95–100%, data not shown), which indicates that the infection load was quite high. The data for M%, a%, and A% are presented in Figure 6. The data on mycorrhization of cv. Finale are consistent with the results of previous studies conducted in similar conditions [33]. The M% values in cultivars k-1693, k-8274, and Finale were significantly lower than those in k-925, k-3064, and k-3358. In the k-8274 cultivar, M% was also significantly lower than in the k-7128 cultivar. Moreover, in k-8274 (earlier characterized as highly symbiotically effective, or high-EIBSM), the value of this parameter was more than three times lower than that in the low-EIBSM cultivars k-3358 and k-3064. Thus, both analyzed low-EIBSM cultivars had a high level of mycorrhizal colonization M%, while the high-EIBSM ones differed in this parameter. The values of the parameters characterizing the development of arbuscules (a% and A%), in general, did not depend on the symbiotic efficiency of the studied cultivars. Thus, the cultivars previously described as high-EIBSM, k-8274, k-1693, and k-925, had a% at the level of that of the low-effective k-3358. Moreover, in the k-3358 cultivar, a% was lower than in the k-7128 cultivar. In cv. k-8274, a% was significantly lower than in cultivars k-7128 and k-3064. The value of the A% parameter in k-8274 was significantly lower than in all other cultivars except for Finale and k-1693. In the cv. Finale, A% was lower than in the k-7128 and k-3064 cultivars, while in cv. k-1693, it was lower than in cv. k-3064. Thus, both analyzed low-EIBSM cultivars had a high level of mycorrhizal colonization, M%, while the high-EIBSM ones differed in this parameter. At the same time, the values of the parameters characterizing arbuscule development (a% and A%), in general, were not related to the symbiotic efficiency of the studied cultivars.

The pea cultivars analyzed for AM development significantly differed in stem length and had different allelic states of the Le gene: le for k-8274 and cv. Finale and Le for the others. As already noted, the Le/le allele status affected relative increments of the values of a number of traits (Figure 4), which implies that plant length might influence pea responsivity to AM-inoculation. We therefore determined the level of Pearson’s correlation between the values of the parameter M% and the values of the stem length (PL) in the studied plants. The level of root mycorrhization had a strong positive correlation with the stem length of the studied pea plants: r = 0.9; p-value < 0.05 (Figure 6d).

4. Discussion

In the present study, we aimed at better understanding the genetic bases of symbiotic responsivity in garden pea (Pisum sativum L.). It is known that plant species and especially varieties within a species differ in their capacities to respond to beneficial soil microorganisms (BSM); this implies the existence of genetic factors determining this responsivity [39]. However, the nature of these factors, as well as mechanisms that underlie the successful interaction of plants with BSM, are currently unclear, mainly because of the complexity of plant–microbial symbiotic systems and difficulties in the assessment of their performance (because the benefits may appear in different forms depending on circumstances and plant properties). Some plant traits can be considered important for plants’ symbiotic effectiveness.

First, host–symbiont specificity plays a key role at least at the very early stage of symbiosis formation-partner recognition. In nitrogen-fixing legume–rhizobia symbiosis, plants can limit the set of bacteria able to penetrate their roots, which may become a substantial strategy in some environments. This is the case for some specimens of pea originating from several Middle Eastern countries, which are able to form nodules only with bacteria producing their signal molecule, which is called nodulation (Nod) factor, with a specific modification [40]. We recently showed that amino acid changes in LysM receptor-like kinase LykX corresponded well with the specificity of interaction with rhizobia strains [41]. The results of mathematical modeling of ligand–receptor interactions of Nod factor and the heterodimeric receptor complex LykX–Sym10 also bespoke the importance of LykX sequence in determination the specificity of interactions of pea with rhizobia [42]. Genes encoding other LysM receptor-like kinases (LysM-RLKs), Sym37 and K1, are located in a cluster together with LykX [27] and also play a role in recognition of Nod factor structure [43,44,45]. The structural variations in LysM domains of the proteins Sym37 and LykX are associated with the specificity of Nod factor recognition and, therefore, have a direct effect on LR symbiosis development [27,41,46].

Arbuscular mycorrhiza generally has no plant host specificity, although some works have reported that specificity of AM fungi might be a key factor in garnering the benefits thereof. Therefore, fungal communities in the roots of medic (M. truncatula) and leek (Allium porrum L.) were usually dominated by one AM species when inoculated by multispecies mixtures, although the composition of the communities depended on both the plant and the time of harvest [47]. AM fungi excrete signal molecules, mycorrhization (Myc) factors, that have a structure similar to that of Nod factors [48] and are presumably recognized in plants by homologs of Nod factor receptors [49]. However, it is unknown whether the structures of any of plant receptors influence the effectiveness of AM symbiosis.

In our work, we found that single nucleotide variants (SNVs) of the genes of Nod factor receptors LykX and K1 possibly influenced symbiotic responsivity to combined inoculation with the rhizobial strain RCAM 1026 and the AM isolate BEG144. The cultivars carrying specific alleles of these genes demonstrated statistically significant increases in the thousand-seed weight (TSW) parameter under Rh + AM inoculation as compared to single Rh inoculation. This result was surprising for us, since this increase was obviously due to mycorrhization, not nodulation. However, phosphorus deficiency in soil (and absence of mycorrhiza in Rh treatment in our experiment) suppresses nitrogenase activity in nodules and, consequently, the nitrogen supply of plants [50], whereas normalization of phosphorus supply (by addition of AM fungi) allows the manifestation of different facets of efficiency of nitrogen-fixing symbiosis. This can explain the observed association of alleles of genes controlling features of legume–rhizobial symbiosis with values of growth parameters under double inoculation (Rh+AM).

In line with our observations, in two independent QTL studies performed in field conditions (i.e., where symbioses with BSM could naturally form), the locus determining the TSW in pea was positioned at the upper part of linkage group (LG) I, which exactly matched the location of a cluster of LysM-RLK genes [13,51]. It seems feasible to analyze the sequences of parental genotypes used in the mentioned QTL studies in the future in order to check whether they differ in the same SNVs that were identified in our work.

Second, the systems of symbiosis regulation by the host must be critical for the symbiosis’s efficiency. Any symbiotic relationship is a risky and complicated trading process in which for each partner it is more expedient to get more than to give. A plant might invest in symbioses and not benefit from them because of a partner’s cheating or objective resource limits [52]. Excessive nodulation has been found to have a negative impact on plant growth, since nitrogen fixation is energy-intensive [53,54]. Legume plants have evolved mechanisms involving long-distance root–shoot signaling that can correct the number of forming nodules, so-called autoregulation of nodulation (AON). The same regulatory elements were shown to be involved in the autoregulation of mycorrhizae (AOM) as well [55,56]. In our study, we did not identify any association between the sequence variations of an important AON gene, Sym29 (CLAVATA1), which controls both nodulation and mycorrhization intensity, and symbiotic responsivity of pea genotypes. Probably, the study of genes encoding CLE peptides, acting as ligands of CLAVATA1 receptor [57,58], could shed light on the involvement of the AON system in the manifestation of symbiotic responsivity in pea.

A plant’s control over the symbiotic process is closely linked with its shoot and root architecture (habitus), which also appears to influence the plant’s responsivity to BSM. As the AON/AOM system’s control assumes the exchange of signal molecules between shoots and roots [55,58], short plants can be expected to react to environmental challenges more efficiently because of a shorter “signal line” and quicker communication between shoot and root. Thus, proper control over AM colonization could be the key to symbiotic efficiency and responsivity. Our data favored this assumption: the mycorrhization level of taller plants was higher than that of shorter plants.

Another example of plant habitus’s impact on symbiotic responsivity is related to the semideterminate stem growth of le⁻ plants. Such plants usually cannot form more than five flowering nodes, unlike Le⁺ plants, the virtually indeterminate growth of which results in the development of many flowering nodes (usually more than ten nodes under optimal mineral nutrition). Recently we suggested that one of the features of the successful interactions of plants with BSM might be the prolongation of the seed-filling period in the so-called highly responsive genotypes [19].

Indeed, previously we demonstrated that double inoculation with Rh + AM increased the TSW of the responsive genotype k-8274, as opposed to that of the nonresponsive genotype k-3358 [59]. The analysis of the dataset obtained from 99 genotypes allowed us to show that shorter plants (le⁻) (this group includes k-8274) demonstrated significantly higher AM-derived increments (both absolute and relative) of TSW and PhS as compared to Le⁺ (this group includes k-3358); these increments were also higher than the average increment calculated for the whole dataset.

Thus, the limited pod number of le⁻ plants determines the strategy of increasing the filling of the already formed limited number of seeds under optimal conditions of mineral supply (in our case, because of the formation of symbioses). Therefore, le⁻ plants are basically predisposed to manifest high responsivity to inoculation in terms of prolongation of the seed filling period and increase in TSW. It is important to note, however, that the le⁻ mutation negatively impacts other characteristics of plants, resulting in lesser increases in total biomass (PW), pod number (PN), and shoot length (PL) under Rh + AM inoculation in le⁻ plants than in Le⁺ plants, and to the average mean calculated for the whole dataset. Also, the fact that the Le gene colocalized with many QTLs in studies involving field trials [10,13] shows that in natural conditions, the symbiotic responsivity mediated by better nitrogen and phosphorus supply of plants may also be manifested.

The last factor we presume to influence symbiotic responsivity is the performance of proteins involved in building the symbiotic compartments and in nutrient exchange. In our work, we discovered that the allele AMT2;3.215N, which encodes a nitrogen transporter, presumably led to more effective acquisition of nutrients when plants were inoculated with AM. Because of complex connectivity between all traits, this enhancement significantly manifested only for the TSW trait. Allelic states of other genes encoding SST1 (symbiotic sulphate transporter 1), IGN1 (ineffective greenish nodules 1, a structural protein that probably anchors other proteins on membranes), STR1 and STR2 (half-size ABC transporters indispensable for arbuscule formation), and SEN1 (stationary endosymbiont nodule 1, an integral membrane protein that possibly transports iron) did not show significant association with the traits’ increments. In L. japonicus, SEN1 allelic states were reported to be responsible for the manifestation of QTLs for nitrogen fixation activity (acetylene reduction activity, ARA) and seed weight [60]. The line Gifu B129 of L. japonicus has a deletion of three nucleotides in comparison with that of L. japonicus Miyakojima MG20, which leads to an absence of an amino acid and, apparently, to lower nitrogen fixation activity and poor seed weight accumulation. In our tested pea cultivars, no amino acid deletions or insertions were found in the putative proteins encoded by the analyzed genes, which probably points to the importance of these genes for proper functioning of symbioses. Notably, in the sequence of the AMT2;5 gene, we were unable to detect any variations (data not shown), which is a possible signal of recent selection pressure that displaced other alleles as not adaptive. Deeper analysis of wild pea varieties may result in the discovery of other alleles of nonvariable symbiotic genes, which should be tested in the future for their influence on the properties of beneficial symbioses formed by pea.

Finally, we suppose that the influence of alleles of symbiotic genes on plant growth parameters was hardly detectable because different varieties in our dataset may have employed different means to increase productivity in symbiotic conditions. Therefore, the particular mechanisms and genes involved in these plants can be elucidated in genetic analysis of progenies of crosses between different genotypes, contrasting in their symbiotic responsivity. We can recommend the crosses of genotypes with similar habitus but contrasting in AM-derived increments by traits of interest, for example, seed biomass (see Figure 7). Also, it is possible that the influence of allelic states of particular symbiotic genes on growth parameters manifested differently at the Le⁺ and le⁻ backgrounds. As our set of lines included only 10 le⁻ genotypes, which seems to be too few for correct testing of this idea, we plan to extend the number of le⁻ genotypes for further studies of symbiotic responsivity in pea. Lastly, since in our recent work we described transcriptome signatures characteristic for symbiotically effective pea cultivars and identified a set of expression markers [20], it seems promising to analyze the promoter regions of the differentially expressed genes on the same set of 99 pea cultivars to identify DNA markers of symbiotic responsivity in pea.

5. Conclusions

The usage of microbial inoculants instead of, or together with, mineral fertilizers is now becoming more and more common. Consequently, the need for new cultivars responsive to inoculation with microbial preparation is increasing, and the responsivity of plants to inoculation with beneficial microorganisms has become an important trait for modern breeding [4]. The employment of state-of-the-art strategies of molecular breeding aimed at creation of highly symbiotically responsive cultivars will be possible after the molecular genetic bases of this trait have been elucidated. In this study, we showed that in pea, the natural mutation in the gibberellin 3-beta-dioxygenase (Le) gene increases the responsivity to inoculation with nodule bacteria and AM fungi so that the increase in individual seed weight under inoculation with Rh + AM is more pronounced in le⁻ than in Le⁺ plants. This fact favors the usage of complex microbial inoculants for modern cultivars of pea, since most such cultivars have short internodes due to a natural Le.A229T mutation and, therefore, should benefit from symbioses with rhizobia and AM fungi. At the same time, the identification of other genes that positively influence symbiotic responsivity and introduction of their “effective” alleles into genomes of new cultivars will be the task of future studies.