Genotypic composition and performance of pea-nodulating rhizobia from soils outside the native plant-host range

Cultivated soils need to shelter suitable rhizobia for legume cropping, especially in areas outside of the plant-host native range, where soils may lack efficient symbiotic partners. We analyzed the distribution patterns and traits of native rhizobia associated with Pisum sativum L. in soils of Hebei Province, a region that has recently experienced an expansion of pea production in China. A total of 43 rhizobial isolates were obtained from root-nodules and characterized genetically and symbiotically. The isolates discriminated into 12 genotypes as defined by PCR-RFLP of IGS DNA. Multiple locus sequence analysis (MLSA) based on the 16S rRNA, recA, atpD and gyrB of representative strains placed them into five clusters of four defined species (R. sophorae, R. indicum, R. changzhiense, and R. anhuiense) and a novel Rhizobium genospecies. R. sophorae was the dominant group (58%) followed by R. indicum (23%). The other groups composed of R. changzhiense (14%), R. anhuiense (1 isolate) and the new genospecies (1 isolate), were minor and site-specific. Based on nodC phylogeny, all representatives were intermingled within the symbiovar viciae with R. sophorae and R. changzhiense being a new record. All the tested strains showed efficient symbiotic fixation on pea plants, with half of them exhibiting better plant biomass performance. This suggests that the pea-nodulating rhizobia in Hebei Province form a specific community of efficient symbiotic rhizobia on pea, distinct from those reported in other countries.


Introduction
Pea (Pisum sativum L.) is one of the most significant pulse grains around the world, accounting for 16% of global pulse production (Drew et al., 2012) due to the high protein content (up to around 30%) and quality of its seeds (Hacisalihoglu et al., 2020). The pea is believed to be native to the Middle and Near East, where domestication began about 10,000-9,000 years B.C. It then spread progressively to other countries with settlements (Smýkal et al., 2012). Chinese peas are thought to have originated from domestication in the Middle East, and introduced in China via an ancient trade road through southern Asia (Wu et al., 2017). As an exogenous species, the pea has been widely cultivated in China for over 2,000 years. Currently, it is grown on an area of about 2.4 million ha with a total production of 12.9 million t of both dry and fresh peas per year (FAOSTAT, 2022).
This legume species has been recorded as a nodule-forming plant since the 1880s in North European soils (Nutman, 1987) and its microsymbiont, Rhizobium leguminosarum, has been the earliest reported member of rhizobia studied since that time. Soil-dwelling rhizobia are capable of fixing dinitrogen in symbiosis with legumes by inducing undetermined nodules on pea roots (Oldroyd et al., 2011). This symbiotic relationship is based on a balance of nutrient exchanges between the partners. Plants supply carbohydrates (primarily succinate and malate, products of photosynthesis) and microhabitats to compatible bacteria, while bacteria differentiated into bacteroids supply ammonia [products of biological nitrogen fixation (BNF)] by secretion to plant hosts (Oldroyd et al., 2011;Rogel et al., 2011;Udvardi and Poole, 2013;Schulte et al., 2021).
The oldest described and most studied species R. leguminosarum was further classified in three symbiovars (sv.) depending on the host plants and symbiotic gene phylogeny: sv. viciae infects plants of tribe Vicieae, like pea and vetches (Vicia spp.); sv. trifolii forms nodules with clovers (Trifolium spp.) and sv. phaseoli, with common beans (Phaseolus vulgaris L; Rahi et al., 2020). Subsequently, pea-nodulating rhizobia were found in species Rhizobium indicum (Rahi et al., 2020), Rhizobium ruizarguesonis (Jorrin et al., 2020), Rhizobium pisi (Ramirez-Babena et al., 2008), Rhizobium anhuiense (Zhang et al., 2015), Rhizobium laguerreae (Flores-Félix et al., 2020), and Ensifer meliloti (Ilahi et al., 2021), all belonging to the sv. viciae. Except E. meliloti, R. pisi, R. indicum and R. anhuiense, all other species belonged to the R. leguminosarum complex (Rlc; Young et al., 2021), a group of closely related species that could only be differentiated by phylogeny of multilocus (or concatenated) sequence analysis (MLSA) of housekeeping genes, but not 16S rRNA genes (de Lajudie et al., 2021;Gurkanli, 2021). The presence or relative occurrence of these species varied according to different geographical aereas often at region or country scales depending on soil and crop history. For instance, R. laguerreae was found dominant in both northwest Spain (100% of the isolates) and Tunisia (gsR, 56%;Flores-Félix et al., 2020;Ilahi et al., 2021), while it formed a widespread minor group in Turkey (12.5%; Gurkanli, 2021). In Turkey, the most common pea microsymbiont was the R. leguminosarum complex genospecies B (Rlc gsB, 47.5%) followed by R. ruizarguesonis (Rlc gsC, 20%). In addition, genospecies Rlc gsA and Rlc gsE were found as sitespecific groups nodulating pea in Turkey (Gurkanli, 2021) and were not recovered in Tunisia (Ilahi et al., 2021). These previous studies demonstrated that diversity and uneven distribution of rhizobial communities associated with pea plants showed distinct geographic patterns, shaped by the adaptation of rhizobia to local conditions that may impact symbiosis interaction and functioning .
Application of chemical fertilizers has effectively enhanced crop yields, but their long-term excessive use has led to increased costs, reduced fertilizer efficiency in agricultural production, as well as severe environmental and biodiversity degradation Olivares et al., 2013). Consequently, the employment of rhizobial strains as inoculants to replace commonly used N-fertilizers has become of primary importance and a goal to be achieved for economic and environmental implications (Mia and Shamsuddin, 2010;Bhardwaj et al., 2014). In recent years, the exploitation and improvement of BNF of field pea is a key element for eco-sustainable agriculture (Holt-Gimenez and Altieri, 2013) based upon its worldwide cultivation and high BNF efficiency. It was estimated that more than half of the biologically fixed nitrogen worldwide is yielded by Rhizobium-legume symbioses (Coyne, 2005;Erana, 2021). BNF in pea crops can reach one of the highest levels with more than 80% of the N provided to plants, while an average of only 25-35 kg/ha of N is introduced in the soil, depending on the tillage and cropping system. Therefore, isolation and phenotypic screening of rhizobial strains to select as bio-inoculants is a crucial strategy. Such screening is also important because the efficiency of the symbiosis can be affected by many other factors, like host-specificity and the ability of the selected strain to compete with local rhizobia (Abi-Ghanem et al., 2013;Bourion et al., 2018). The plant growth benefits provided to the same legume species by different rhizobial strains at a given location can vary up to 10-fold (Denison and Kiers, 2004).
Considering all the aforementioned aspects, and the fact that pea-nodulating rhizobia in China have not been systematically studied since the extensive revision of the Rhizobium and Rlc taxonomic system, we decided to conduct the present study. The aim of this work was to evaluate the diversity, relative abundance, and geographic distribution of native rhizobia that nodulate P. sativum in Hebei Province, a non-native area of pea. Traditionally, the pea has been grown sparingly in Hebei as an economic crop for fresh vegetables and grains, and the pea cultivation area has been rapidly expanded there to meet consumer demand. For this purpose, the taxonomic status of the isolated strains from root-nodules was determined through ribosomal intergenic typing, phylogenies analyses of housekeeping genes (recA, atpD, and gyrB), the 16S rRNA gene and a symbiotic gene (nodC). Additionally, the distribution of rhizobia in relation to soil properties and the potential of representative strains to induce effective symbiosis, were also investigated.

Materials and methods
Field and soil sampling and soil characterization Soils were sampled from three fields cultivated with local pea in the Guyuan County, Hebei Province of China. These sites were at the Guyuan Institute of Crop Science, Chinese Academy of Agricultural Sciences (E 115°39′, N 41°40′12″, altitude 1,410 m); Xiaochang town (E 115°46′48″, N 41°24′36″, altitude 1,519 m); and Dashila village (E 115°54′, N41°30′36″, altitude 1,415 m). At each site, soils were sampled to a depth of 10-20 cm from the fields during the flowering stage of Pisum sativum in June 2020. For each site, five randomly taken soil sub-samples of equal volume were thoroughly mixed to constitute a representative soil sample of the field site. Then, soil samples were crushed to a uniform state, and transported to the laboratory in an ice-filled cooler (Zhang et al., 2018). Part of each representative soil sample was chemically analyzed for pH, total salts (TS), organic matter (OM), available phosphorus (AP), available potassium (AK) and total nitrogen (TN) as described previously (Appunu et al., 2009;Han et al., 2009).

Rhizobial isolation and conservation
Surface-sterilized seeds (2.5% w/v NaClO solution for 5 min) of pea (variety of Zhongqin No.1) were germinated and the seeding was Frontiers in Microbiology 03 frontiersin.org sown in surface-sterilized plastic pots (15 cm high × 10 cm diameter) filled with a soil-sterilized vermiculite mixture (1/5 vol/vol) as described previously (Zhang et al., 2019). All plants were grown under greenhouse conditions of 25/20°C (day/night) with a 16-h photoperiod. Sterilized water was added to the pots throughout the experiment as required. After 45 days, all plants in all soils were up-rooted and bacterial strains were isolated from nodules according to the standard protocol (Zhang et al., 2019). Five plants from each soil treatment were selected for the collection of root nodules and isolation of rhizobia. Root nodules were surface sterilized and the individual sterilized root nodule was crushed in sterile water and the bacterial suspension was streaked onto a Yeast extract-Mannitol Agar (YMA) plate. After incubation at 28°C for 3 to 5 days, single colonies representing the dominant bacteria in each plate were picked up and purified by cross-streaking on new YMA plates until pure cultures were obtained. All purified isolates were conserved in tryptone yeast (TY) broth (tryptone 5 g; yeast extract 3 g; CaCl 2 0.6 g; distilled water 1 L, pH 7.0) supplied with glycerol (20%, v/v) at −80°C for long-term storage and on YMA slants at 4°C for temporary storage.

Genomic fingerprinting of Rhizobial isolates
Genomic DNA (gDNA) of each rhizobial isolate was extracted according to Terefework et al. (2001). The gDNA was estimated qualitatively and quantitatively by using a nanodrop (Thermo Fisher) and high quality gDNA was used as template for PCR amplifications of the 16S-23S rRNA intergenic spacer (IGS) region with primers IGS1490' (forward) and IGS132' (reverse; Laguerre et al., 1996). PCR amplification was carried out in a standard 50 μL reaction mixture including 1 μL of DNA template and 5 U of Taq DNA polymerase (Sangon Biotech (Shanghai) Co., Ltd.). Aliquots of amplified PCR products were visualized after electrophoresis in a 1.0% (w/v) agarose gel labeled with GoldView type I. Then PCR products were digested separately with the endonucleases HaeIII, MspI and HhaI (Laguerre et al., 1994) at 37°C for 10 h. The 16S-23S rDNA IGS type of each strain was designated after separation and visualization of restriction fragments by electrophoresis in 2.5% (w/v) agarose gel and UV-illumination.

Molecular and phylogenetic identification of the isolates
Isolates sharing the same RFLP pattern of 16S-23S rDNA IGS in this study were designed as an IGS type. One representative strain for each IGS type in each sample site was selected for amplification of the 16S rRNA gene using the forward primer P1 and reverse primer P6 (Laguerre et al., 1996). The PCR products were verified as mentioned above, and were sent for commercial sequencing (Sangon Biotech (Shanghai) Co., Ltd.). The acquired sequences were compared in the NCBI database using the online BLASTN tool and the sequences for type strains of defined Rhizobium species sharing similarities greater than 97.0% with the new isolates were extracted. The phylogenetic analysis was conducted in the MEGA 7.0 software (Kumar et al., 2016). Sequences were aligned using Clustal W and the best model of sequence evolution was selected. Then the phylogenetic tree was inferred by using the Maximum-likelihood (ML) and the non-parametric bootstrap methods (500 pseudo-replications). DNA fragments of recA (coding for DNA recombination protein), atpD (encoding for ATP synthase beta chain) and gyrB (structural gene for DNA gyrase beta subunit) were amplified separately by PCR using the primer pairs recA41F/recA640R, atpD255F/atpD782R and AMgyrBf/ AMgyrBr, respectively, (Vinuesa et al., 2005;Stepkowski et al., 2012). The nodulation gene nodC was amplified by using the primer pair nodC540F/nodC1160R (Sarita et al., 2005). The PCR conditions were adopted from previous reports (Sarita et al., 2005). PCR products verification, sequencing and tree construction of each sequenced gene were performed as mentioned above. Furthermore, sequences of atpD, recA and gyrB genes were concatenated and aligned using Clustal W (Thompson et al., 1997). Distance calculation, construction of the concatenated housekeeping gene ML tree, bootstrap analysis were performed in MEGA 7.0 as described above. The sequences resulted from the alignment of the related genes have been deposited in the NCBI database (accession numbers are indicated on the trees,

Correlation analysis between soil properties and Rhizobial communities
A distance-based redundancy analysis (RDA) was performed using CANOCO software version 4.54 to investigate the relationships between soil properties (TN, AP, AK, OM, TS, and pH) and the rhizobial community composition based on IGS genotypes.

Nodulation, cross-inoculation and symbiotic efficiency test
The representative strains of the IGS types were tested for their ability to nodulate their original host plant P. sativum and also Vicia faba L., Vicia sativa L., Vigna radiata L., Glycine max L., Arachis hypogaea L. and Cicer arietinum L. Symbiotic efficiency of the strains was tested on pea. Briefly, surface sterilized seedlings were aseptically transferred in pots (1 plant/pot) containing sterile vermiculite as substrate and inoculated with 1 mL of rhizobial suspension (OD 600 = 1.0). Plants were grown under greenhouse conditions and watered with sterile N-free nutrient solution as required . Nodulation was evaluated 40 days after inoculation for all plant species above. Pea growth was estimating by weighting their dry root and shoot biomass and measuring leaf chlorophyll contents (SPAD chlorophyll meter). Uninoculated plants as a negative control, were included and all plant treatments were performed in three replicates. Data were analyzed by one-way ANOVA followed by an LSD post hoc test (p = 0.001).

Physicochemical characteristics of soils
All three sites differed from each other in their levels of organic matter, total nitrogen, available phosphorus, and available potassium (

IGS PCR-RFLP analysis
Totally, 43 rhizobial isolates were obtained from in 3 locations. The 43 isolates were distinguished into 6, 8, and 7 RFLP patterns obtained from the restriction enzymes HaeIII, MspI and HhaI, respectively. By combining all RFLP patterns, the IGS PCR-RFLP fingerprinting allowed the classification of all isolates into 12 IGS types (Table 2). Among the 12 IGS types, IGS type 1 represented the most abundant population (with 20 isolates), followed by IGS type 2 (5 isolates), type 3 (4 isolates), types 4 and 5 (3 isolates each), and type 6 (2 isolates). All the other IGS types have just 1 isolate. Isolates of IGS type 1 distributed over all the three sampling sites, with the H-ES site having the largest number (15 isolates, 34.9%) and the H-CX the smallest (1 isolate, 2.3%). IGS type 2 isolates were detected at sites H-DSL (3 isolates, 7.0%) and H-XC (2 isolates, 4.7%; Table 2). IGS types 3, 5, and 7 were only distributed in site H-XC (8 isolates, 18.6%); IGS types 4 and 8 were only found in H-DSL (4 isolates, 9.3%) and IGS types 6, 9, 10, 11, and 12 were only distributed in H-ES (6 isolates, 14.0%). Site H-ES was inhabited by the most diverse P. Maximum-likelihood phylogenetic tree based on 16S rRNA gene sequences (1,217 base pairs) from studied Pisum sativum L. rhizobia and related rhizobial sequences retrieved from Genbank. The tree was inferred under the best-fit model (T92 + G + I) and rooted with Bradyrhizobium americanum CMVU44 T . Scale bar indicates 0.02 nucleotide substitutions per site. Bootstrap confidence values (%) calculated for 500 replications >70% are shown at the internodes.

Phylogenetic analysis of core and symbiotic genes
Nearly full-length 16S rRNA genes were successfully amplified and sequenced for 13 rhizobial isolates representing all 12 IGS types and sites of origin (  (Figure 1).
The representative isolates divided into five clades (C1-C5) in the phylogenetic tree based on their concatenated recA-atpD-gyrB sequences (Figure 2). Five representative isolates from IGS types 1, 4, 7, and 11 grouped together in C1 cluster which also included the type strain of R. sophorae. They shared 98.5-99.8% similarities with each other and similarities ≤97.0% with all type strains of other closely related Rhizobium species (Table 2). Thus, cluster C1 was identified as R. sophorae. Cluster C2 contained R. indicum JKLM 12A2 T and IGS types 2, 5, and 6 representing 10 isolates (23% in total) which shared 97.3-98.1% similarities with each other and less than 96.7% similarities with the other Rhizobium type strains; they were therefore identified as   Figure S4) and was therefore assigned to a novel Rhizobium genospecies. The phylogenetic analyses of the single gene of recA, atpD, and gyrB were shown in Supplementary Figures S1-S3.
The nodC genes were amplified from all the 13 representative isolates, confirming their genetic basis for rhizobial symbiosis. As shown in the nodC phylogeny, three strongly supported groups (N1-N3) were defined among the representatives (Figure 3) (Figure 3).
The nodC and concatenated recA-atpD-gyrB gene trees were not congruent. Symbiotype N1 covered four representatives belonging to C1 (R. sophorae) and C2 (R. indicum) with 100% nodC phylogeny similarities. Representative isolates of symbiotype N2 that displayed 100% nodC phylogeny similarities, belonged to C1 or C3. The single representative of C4 or C5 both belonged to the symbiotype N3 while other N3 strains were also part of C2 and C3 (Figure 2; Supplementary Figure S3; Table 2). It is worth noting that the three nodC types showed a variable distribution among the sampling sites: all representative strains isolated from H-DSL clustered in symbiotypes N1 to N3, representatives from H-XC in symbiotypes N1 and N3, and representatives from H-ES in symbiotypes N2 and N3. In this study, representatives of R. changzhiense (C3) were found in two nodC gene

Correlation analysis between Rhizobial distribution and soil properties
Redundancy analysis was used to explore the relationships between soil chemical properties and the rhizobial community composition based on IGS genotypes. The analysis was limited to 3 field sites with a narrow pH range (pH 7.1-7.5) and similar total salts across the sites. Maximum-likelihood phylogenetic tree based on symbiotic gene nodC (355 base pairs) showing the relationships of the rhizobia isolated from nodules of Pisum sativum L. in Hebei Province of China. The tree was constructed using the maximum-likelihood method under the best-fit model (T92 + I). Scale bar indicates 0.05 nucleotide substitutions per site. Bootstrap confidence values (%) calculated for 500 replications >70% are indicated at the internodes.
Frontiers in Microbiology 08 frontiersin.org Nevertheless, the sites differed from each other, in particular, organic matter, total nitrogen, available phosphorus and available potassium were lower at Xiaochang town than that at Dashila village and the experimental site of the Guyuan Institute of Crop Science. According to RDA results (Figure 4), all ribosomal IGS types were separated into five distinct groups (Groups I to V) and soil chemical factors had different effects on the distribution of these rhizobial groups and IGS types. Group I includes the IGS types 4 and 8 encompassing 4 isolates. Its distribution was positively correlated with higher OM, AK, and TN values (Figure 4) that coincided with the soil properties of H-DSL (Table 2). Thus, nutriments with higher OM, TN and AK content (Table 1) appeared more suitable for the survival or P. sativum nodulation of Group I rhizobia. Group II (including IGS type 1) gathered the majority of isolates (20) distributed across the three sites (Table 2): it presented a positive association with low TS contents, higher AP, but intermediate OM, AK, and TN values ( Figure 4). Group III was formed by five IGS types (6, 9, 10, 11 and 12) in a balanced population size (6 isolates) and found only in H-ES (Table 2). Meanwhile, strains belonging to Group III were correlated with higher AP and pH, but lower TN, AK, OM and TS (Figure 4). Group IV comprised rhizobia characterized by IGS types 3, 5 and 7, and trapped exclusively from site H-XC (Table 2): it was associated with higher pH and TS contents but lower OM, AK, TN, and AP values (Figure 4). Finally, isolates of Group V (including only IGS type 2 of sites H-XC and H-DSL) were positively distributed with higher TS, but lower AP and pH values. Thus, we found that the different rhizobial populations were influenced differently by edaphic factors, with a potential positive or negative correlation. According to the IGS type biogeography obtained in this study, the distribution of most strains belonging to R. sophorae (20/25 were in Group II) might be more suitable to higher AP, TN, AK and OM, while rhizobia of R. changzhiense (Group III and IV) could be favored by higher levels of soil pH. Besides, no obvious preference was shown for R. indicum (inside Groups III, IV, and V; Figure 4; Tables 1, 2). Thus, rhizobial species seemed to have different responses to soil chemical characteristics. Furthermore, the distribution pattern found at the scale of IGS genotypes may have occurred in the soil independently of that of the rhizobial species.

Discussion
Although pea plants have been cultured widely in different ecoregions in China, the diversity and geographic distribution of pea-nodulating rhizobia in China have rarely been studied (Ramirez-Babena et al., 2008). Previously, R. leguminosarum, two unnamed Rlc genomic species and R. etli were detected using the phylogeny of 16S rRNA genes and PCR-amplified ribosomal IGS restriction patterns (Yang et al., 2008). More recently, R. anhuiense was described based on a polyphasic study (Ramirez-Babena et al., 2008) of pea-nodulating rhizobia in two subtropical regions of China (South East). Unlike the previous studies, we systematically investigated the diversity of P. sativum nodulating microsymbionts across three sampling sites within a temperate region of Hebei Province (North East China).

FIGURE 4
Symbiotic performance of 13 representative strains evaluated on pea grown under greenhouse conditions 40 days after inoculation, CK is the negative control (uninoculated plants). The experiment was conducted in triplicates. In addition, three subsamples were performed for each SPAD replication.  Taking all the results of the 16S rRNA sequence analysis, the phylogeny of concatenated recA-atpD-gyrB gene sequences, IGS PCR-RFLP analysis, and the nodC gene phylogeny, all the isolates obtained in this study were classified as a diverse community consisting of 12 IGS types within 5 genomic species of the genus Rhizobium, all of which belong to the sv. viciae. The dominant species of pea-nodulating rhizobia in the three sampling sites were R. sophorae (58% in relative abundance) and R. indicum (23%), which were found to be distributed universally across all three sites. In contrast, R. changzhiense, R. anhuiense, and a novel Rhizobium genospecies within the Rlc were identified as minor groups and were site-specific. Based on the concatenated housekeeping gene tree, most of the isolates including the new R. changzhiense ones, belong to the complex Rlc (Olivares et al., 2013;de Lajudie et al., 2021;Young et al., 2021;Zhang et al., 2021) with only two isolates falling outside of it (R. anhuiense and the new genospecies C5). This provides new information on the diversity and geographic distribution of pea-nodulating rhizobia. Firstly, a unique community composition of pea-nodulating rhizobia was found in the temperate area studied in China, compared with previous reports in the subtropical zone in China (Ramirez-Babena et al., 2008) and other countries (Bergey et al., 1984;Ramirez-Babena et al., 2008;Jorrin et al., 2020;Rahi et al., 2020;Gurkanli, 2021;Ilahi et al., 2021;Young et al., 2021). Secondly, the identification of R. sophorae and R. changzhiense as novel pea-nodulating rhizobia expands the range of known rhizobial species associated with this plant. R. sophorae was initially isolated from effective nodules of the shrubby Sophora (Sophora flavescens) in Changzhi City (Shanxi Province, China), and was found to be able to effectively nodulated not only Sophora but also Phaseolus vulgaris (Jiao et al., 2015). Later, this rhizobial species was also recovered from root nodules of Vicia faba L. grown in Panxi (southwest China; Chen et al., 2018), southwest China (Zhang et al., 2019), and Hebei Province (northeast China) (Zhang et al., 2022). Our results extend the promiscuous nitrogen-fixing host plant list for R. sophora to include P. sativum. Together with the above-mentioned studies, these findings confirm the widespread existence of sv. viciae in R. sophorae. This emphasizes the wide distribution of this rhizobial species over space and time, as well as its high adaptability to soil and environmental factors in China. R. changzhiense was first isolated from V. sativa cultivated in the Changzhi city of China (Zhang et al., 2021). Together with our results, it provides evidence that rhizobial strains originated from both Vicia sativa (Zhang et al., 2021) and P. sativum belonged to the same species and symbiovar. R. indicum was isolated from root nodules of pea (P. sativum) cultivated in the Indian trans-Himalayas (Rogel et al., 2011). Meanwhile R. anhuiense was first found in rhizobia of pea and V. faba in Anhui and Jiangxi Provinces (Li et al., 2016), and has also been widely recorded in Shandong Province (Li et al., 2016). Our results in the present study might demonstrate that these species have a wide geographic distribution in the world and that pea plants have selected the chromosomal background or communities adapted to the local conditions during the long-term legume cultivation. Thirdly, the detection of highly conserved nodC genes across the 5 genomic species identified in this study offered further evidence that pea plants have stringently selected symbiosis genes in their microsymbionts. This selective pressure may favor symbiosis genes toward the most adapted indigenous rhizobia, as reported in other cases (He et al., 2011). Fourthly, our results also imply the necessity to screen out high-quality rhizobial strains with strong adaptability to local conditions to produce inoculants to increase the yield of pea or other legumes.
The levels of OM, TN, AP, and AK among the three field soils sampled were significantly different in this study, and the distribution of IGS types was related to these soil characteristics. We also found that the same species isolated from different soil types had different responses to soil physicochemical characteristics, as evidenced by a distinct pattern observed at the IGS genotypes within species ( Figure 5). Soil chemical factors had varying effects on the distribution of the rhizobial species and IGS types. For instance, we observed a positive correlation between IGS type 1 (R. sophorae) and a negative correlation of IGS types 3 (R. changzhiense), 5 (R. indicum), and 7 (R. sophorae) with AP, which may influence their separation into Groups II and IV ( Figure 5). Similarly, a positive correlation of IGS type 2 (R. indicum) and a negative correlation of IGS types 6 (R. indicum), 9 and 12 (R. changzhiense), 10 (Rhizobium sp. I), and 11 (R. sophorae) with TS influenced their separation in two distribution patterns. These relationships were consistent with previous reports that soil traits were selective abiotic factors for the biogeography of rhizobial species or IGS types (Han et al., 2009).
In conclusion, our study demonstrated the existence of indigenous rhizobia that form an effective symbiosis with pea cultivated in North China, in which R. sophorae and R. changzhiense were new records as pea-nodulating rhizobia. R. sophorae was highly abundant in the studied soils (58%) and formed a unique species assemblage of pea-nodulating rhizobia along with four less frequent species or genospecies (R. indicum, R. changzhiense, R. anhuiense and a novel genospecies). Moreover R. sophorae and R. indicum occurred in all tested soil types, indicating that strains of these species could be better competitors and adapt to different soil conditions. Finally, all strains belonged to the sv. viciae, Redundancy analysis (RDA) to relate the distribution of the 12 IGS types of isolates to physicochemical factors of soils collected from the different sites. Blue arrows indicate IGS types of rhizobia and red arrows represent soil properties. The longer the arrow was, the greater the influence of the soil property presents on the distribution of the IGS types. The smaller the angle between the arrow and the IGS type was, the stronger the effect of the soil property on distribution of the IGS type.
Frontiers in Microbiology 10 frontiersin.org regardless of their species affiliation. The genetic diversity of pea-associated rhizobia from Hebei Province provided a significant guideline for local pea cultivation and breeding efforts, as well as for the site-specific selection of efficient rhizobial genetic types.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.