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).

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 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₂ 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 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.

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, Figures 1, 2; Supplementary Figures S1–S5).

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.

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₆₀₀ = 1.0). Plants were grown under greenhouse conditions and watered with sterile N-free nutrient solution as required (Zhang et al., 2012). 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).

All three sites differed from each other in their levels of organic matter, total nitrogen, available phosphorus, and available potassium (Table 1). Field soil at H-DSL contained the highest contents of OM (41.4 g/kg soil), TN (2.5 g/kg soil) and AK (195 mg/kg soil) while soil from site H-ES contained the highest AP (57.8 mg / kg soil). At the opposite, site H-XC exhibited the lowest proportions of organic carbon and nitrogen, and available P and K minerals in soil. Two of the sites had weakly alkaline soil pH values of 7.5, while the remaining site had a neutral pH of 7.1. All three sites had comparable and very low total soluble (TS) contents of 0.6–0.7 g/kg of soil.

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. sativum-nodulating rhizobial community with 6 IGS types, while sites H-XC and H-DSL harbored 5 and 4 IGS types, respectively, (Table 2), indicating that the richness and evenness of rhizobial IGS types varied across the sampling area.

Nearly full-length 16S rRNA genes were successfully amplified and sequenced for 13 rhizobial isolates representing all 12 IGS types and sites of origin (Table 2). All representative isolates clustered together with several defined Rhizobium species in a single clade that showed 99.8–100% similarity in their 16S rRNA gene sequences. This clade comprises the type strains of R. hidalgonense FH14^T, R. anhuiense CCBAU 23252^T, R. indicum JKLM 12A2^T, R. laguerreae FB206^T, R. leguminosarum LMG 14904^T, R. ruizarguesonis UMP1133^T, R. sophorae CCBAU 03386^T, Rhizobium acidisoli FH13 ^T and R. changzhiense WYCCWR 11279^T (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 R. indicum. Cluster C3 (6 isolates from 3 IGS types) showed 97.1–99.2% similarities with each other and with R. changzhiense WYCCWR 11279^T, while these isolates presented less than 95.7% similarities with the other Rhizobium species; so, this cluster was identified as R. changzhiense. Clusters C4 and C5 contained only the IGS type 10 and IGS type 8 respectively, each covering a single isolate. The C4 isolate WYCCWR 13265 was identified as belonging to the species R. anhuiense since it displayed 99.1% similarity with type strain CCBAU 23252^T. Finally, the C5 isolate WYCCWR 13271 had similarities of less than 94.9% with all type strains of the closest defined Rhizobium species (Table 2; Supplementary 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). First, symbiotype N1 clustered together with nodC sequences from strains R. binae BLR195^T, R. fabae CCBAU 33202^T, R. anhuiense CCBAU 23252^T, R. laguerreae FB206^T, Rlv 3,841, and Rlv 248 with similarity values higher than 98.3%. Second, symbiotype N2 contained nodC sequences 100% similar to those of R. indicum JKLM12A2^T and Rlv sv. viciae SWD14-4. Finally, symbiotype N3 encompassed nodC sequences 100% similar to those of R. multihospitium CCBAU 83401^T or R. changzhiense WYCCWR 11279^T, and 99.1% similar to the nodC of R. bangladeshense BLR175^T (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 haplotypes (N2 and N3), while those of R. sophorea (C1) fell into nodC haplotypes N1 and N2, those of R. indicum (C2) into nodC haplotypes N1 and N3, and those of R. anhuinense (C4) and Rhizobium sp. I (C5) into nodC haplotype N3.

All representatives from the 12 different IGS types nodulated Pisum sativum, Vicia faba and Vicia sativa. None of them nodulated Vigna radiata, Glycine max, Arachis hypogaea or Cicer arietinum. Compared to the uninoculated control (p < 0.001), representative isolates of the IGS types increased leaf chlorophyll units (+14–23%) and total plant dry weight (+45–86%) of pea indicating that they were effective rhizobial symbionts. Six strains (WYCCWR13245, −57, −61, 65, −70, −74) tended to show higher plant dry biomass (+72–86%), while four of them (WYCCWR13247, −69, −81, −84) exhibited significantly lower pea dry biomass (+45–53%; Figure 3; Table 3).

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. 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.