Genetic and Physiological Characterization of Soybean-Nodule-Derived Isolates from Bangladeshi Soils Revealed Diverse Array of Bacteria with Potential Bradyrhizobia for Biofertilizers

Genetic and physiological characterization of bacteria derived from nodules of leguminous plants in the exploration of biofertilizer is of paramount importance from agricultural and environmental perspectives. Phylogenetic analysis of the 16S rRNA gene of 84 isolates derived from Bangladeshi soils revealed an unpredictably diverse array of nodule-forming and endosymbiotic bacteria—mostly belonging to the genus Bradyrhizobium. A sequence analysis of the symbiotic genes (nifH and nodD1) revealed similarities with the 16S rRNA gene tree, with few discrepancies. A phylogenetic analysis of the partial rrn operon (16S-ITS-23S) and multi-locus sequence analysis of atpD, glnII, and gyrB identified that the Bradyrhizobium isolates belonged to Bradyrhizobium diazoefficiens, Bradyrhizobium elkanii, Bradyrhizobium liaoningense and Bradyrhizobium yuanmingense species. In the pot experiment, several isolates showed better activity than B. diazoefficiens USDA110, and the Bho-P2-B2-S1-51 isolate of B. liaoningense showed significantly higher acetylene reduction activity in both Glycine max cv. Enrei and Binasoybean-3 varieties and biomass production increased by 9% in the Binasoybean-3 variety. Tha-P2-B1-S1-68 isolate of B. diazoefficiens significantly enhanced shoot length and induced 10% biomass production in Binasoybean-3. These isolates grew at 1–4% NaCl concentration and pH 4.5–10 and survived at 45 °C, making the isolates potential candidates for eco-friendly soybean biofertilizers in salty and tropical regions.


Introduction
Bacterial diversity analysis is important to understand, maintain, and conserve global genetic resources. Bacteria are not only a vital part of soil and environment, but also the main agents of different nutrient cycling [1]. As new environments and regions are explored, the abundance and fruitfulness of microbial diversity will become increasingly evident [2]. Root nodule bacteria, comprising different types of Alpha-proteobacteria, Beta-proteobacteria, and possibly Gamma-proteobacteria, are of significant importance in biodiversity studies [3].

Collection of Soil Samples
Soil samples were collected from 11 districts of different agro-ecological zones in Bangladesh, focusing mainly on areas where soybeans have been cultivated (Table 1 and Figure 1). For each location, soil sample was collected from a depth of~10-15 cm from different positions of a field and mixed. The soil samples were stored at 4 • C until use.

Collection of Nodules Using Soil Samples in Pot Experiment
The seeds of soybean cultivar Glycine max c.v. Enrei (Japanese variety) were sterilized by immersion in 70% ethanol for 30 s, and then in 3% sodium hypoch

Collection of Nodules Using Soil Samples in Pot Experiment
The seeds of soybean cultivar Glycine max c.v. Enrei (Japanese variety) were surfacesterilized by immersion in 70% ethanol for 30 s, and then in 3% sodium hypochlorite solution for 3 min, and the seeds were washed exhaustively with sterile water [19]. For seed germination, the sterilized seeds were incubated for two days at 28 • C. From each soil sample, 1 g of soil was suspended in 5 mL of sterilized water and soil suspensions were used as inoculants. Each inoculant was applied to 300 mL glass jars containing germinated seeds and sterilized vermiculite, and the jars were placed in growth chamber and kept under controlled conditions with 12 h light/dark photoperiod at temperatures of 25 • C and 18 • C in the light and dark hours, respectively [19]. Plant growth was supported by adding a sterilized N-free nutrient solution [20] to the jars. After four weeks, whole plants were uprooted from the jars, washed in running tap water to remove vermiculite, and the nodules were harvested. Root nodules were surface-sterilized by immersion in 70% ethanol for 30 s and in 3% sodium hypochlorite for 3 min, then washed five times with sterile water [19]. Each nodule was crushed in 100 µL of glycerol solution (15%, v/v) to obtain a turbid suspension. An aliquot (10 µL) of suspension was plated onto 1.5% yeast extract mannitol agar (YEM) [21] plates and incubated for 3-7 days at 28 • C. The remaining suspension was maintained at −30 • C. Single colonies were selected and checked for purity by repeated streaking onto fresh YEM medium. The isolates were maintained for the long-term into 50% glycerol stocks at −80 • C and the short-term into slant stocks at 4 • C.

Temperature Tolerance Test of Isolates
The isolates used in this study were examined for growth under high-temperature conditions. For the temperature tolerance test, isolates were inoculated on YMA plates and incubated at 10 • C, 20 • C, 28 • C, 37 • C, 40 • C, and 45 • C for 3-7 days [22]. After incubation, bacterial growth was examined.

pH Tolerance Test of Isolates
The isolates used in this study were examined for growth at different pHs. For pH tolerance tests, isolates were inoculated on YMA plates at pH 4.5, 6, 7, 8.5, and 10. The pH was adjusted to pH 4.5 and 6 using 0.5 M HCl, while the pH was adjusted to pH 7, 8.5, and 10 using 0.5 M NaOH. After bacterial inoculation, the plates were incubated at 28 • C for 3-7 days [22]. After incubation, bacterial growth was examined.

Salinity Tolerance Test of Isolates
The isolates used in this study were examined for growth under high-salinity conditions. For the salt tolerance test, isolates were inoculated on YMA plates containing 1%, 2%, 3%, and 4% (w/v) NaCl concentration and incubated at 28 • C for 3-7 days [22]. After incubation, bacterial growth was examined.

Analysis of DNA Sequences and Phylogeny Genomic
An analysis and quality check of DNA sequences were performed primarily using Sequence Scanner 2.0 software (Thermo Fisher Scientific) and 4Peaks (Nucleobytes). Sequences derived from forward and reverse primers were individually joined by identifying the overlapping sequence between them, using GENETYX version 11 software (Genetyx Corp., Tokyo, Japan) and the online tool Merger (http://www.bioinformatics.nl/cgi-bin/ emboss/merger; accessed on 17 June 2019). Multiple sequence alignments of nucleotide sequences and bootstrapping to create maximum likelihood and neighbor-joining phylogenetic trees were performed using MEGA X [29]. The evolutionary history of all phylogenetic trees in this study was inferred by the maximum likelihood method using MEGA X software [29]. The percentage of trees in which the associated taxa clustered together is shown next to the branches. The initial tree(s) for the heuristic search were obtained by applying the neighbor-joining method to a matrix of pairwise distances, estimated using the maximum composite likelihood approach. The tree was drawn to scale, with branch lengths measured in terms of the number of substitutions per site [29].

Analysis of Inoculation Effects of Selected Isolates on Plants
Isolates were grown in yeast-extract mannitol (YM) broth for three days at 28 • C, after which the cells were collected by centrifugation. The cells were then washed and diluted with sterile, ultrapure Milli-Q water. The absorbance of the diluted cells was measured using a spectrophotometer (Ultrospec 3300 Pro; Amersham Biosciences, Cambridge, UK) at 600 nm. For each isolate, the amount of bacterial inoculation was kept almost equal by maintaining (absorbance at 600 nm × volume of inoculation = 1) the same scale and measuring the colony-forming unit (CFU) by plate count. The concentration of cells in the solution used for inoculation was maintained at approximately 10 8 cells mL −1 [24,25]. Surface-sterilized seeds of Glycine max cv. Enrei and Binasoybean-3 (Bangladeshi soybean variety, Table S2) were sown in 300 mL of sterile vermiculite (Yoshino Gypsum Co., Ltd., Tokyo, Japan) in glass jars. After sowing, 1 mL of the cell suspension was applied to the seeds. Three replicates were performed for each isolate. Subsequently, the jars were transferred to a growth chamber. Sterilized N-free nutrient solution [20] was added to the glass jars. The moisture content was maintained at~70% throughout the growth period. Plants were grown in a growth room under a 12 h light (25 • C)/dark (18 • C) photoperiod. After five weeks of incubation, whole plants were removed and washed from the vermiculite and different parameters, e.g., shoot length, root length, nodule color, nodule size and number, shoot weight, root weight, nodule weight, and ARA, were measured [23,30]. Experiments were performed using a randomized complete block design with 2 soybean varieties, 14 bacteria and 3 replicates for each treatment [31].

Acetylene Reduction Assay
For the acetylene reduction assay (ARA), whole plant roots with root nodules from each tested plant were introduced into a 300 mL airtight jar. Subsequently, 30 mL of air from the jar was replaced with 30 mL of acetylene, and the whole roots with root nodules were further incubated for 1 h at 30 • C to evaluate nitrogenase activity in the root nodules [19]. The amount of ethylene produced owing to the activity of nitrogenase on acetylene enclosed in the jar was determined using a Shimadzu GC-8A gas chromatograph (Shimadzu, Tokyo, Japan) equipped with a Porapak N column (Chrompack, Middelburg, the Netherlands). The roots of plants without bacterial inoculation were used as controls [19].

Statistical Analysis
The experimental data obtained from plant tests were subjected to statistical analyses, such as multiple comparison using general linear model with Dunnett's Post Hoc test (with Microorganisms 2022, 10, 2282 6 of 21 95% confidence level), a correlation analysis using the bivariate correlations method with Pearson Correlation Coefficient and two-tailed test of significance, and regression analysis using the linear regression method with 95% confidence intervals level using IBM SPSS Statistics, Version 23.0.

Isolation of Bacteria from Root Nodule
Soil samples from 11 districts of Bangladesh were used for root-nodule bacterial isolation. Soybean (Glycine max cv. Enrei) plants were grown in vermiculite using the soil samples. A total of 84 isolates were obtained from the selected root nodules of soybean plants.

Phylogenetic Analysis Based on the 16S rRNA Genes
To determine the phylogenetic position of the isolates, 16S rRNA gene sequencing and phylogenetic analysis with highly similar reference strains were performed using a 1348 nt sequence for each bacterium.
Subsequently, the 84 isolates were classified into eight groups ( Figure 2a) belonging to the Alpha-proteobacteria, Beta-proteobacteria, Gamma-proteobacteria, Firmicutes, and Actinobacteria phylum/class. The majority of isolates (65/84; 77.4%), showed a close relationship with the genus Bradyrhizobium, belonging to the class Alphaproteobacteria and family Rhizobiales. The Bradyrhizobium isolates were classified into three groups: A (16S), B (16S), and C (16S). Forty-five isolates showing a close relationship with Bradyrhizobium diazoefficiens were grouped in A (16S) (Figure 2b). Group B (16S), consisting of 16 isolates, was found to belong to the same clade as Bradyrhizobium liaoningense and Bradyrhizobium yuanmingense (Figure 2c). Four isolates belonging to group C (16S) were observed to have the highest similarity to Bradyrhizobium elkanii. Group D (16S), with six isolates, showed a close relationship with the genus Methylobacterium, which also belongs to the class Alpha-proteobacteria and family Rhizobiales.
Two isolates of group E (16S) and seven isolates of group F (16S) were observed in the same branch as the genus Stenotrophomonas of the Gamma-proteobacteria class and genus Pandoraea of the Beta-proteobacteria class, respectively. Group G (16S), comprising only one isolate, was observed to reside in the same branch as Leifsonia lichenia of the Actinobacteria class. Three isolates of group H (16S) belonged to the genus Bacillus and Firmicutes.

Phylogenetic Analysis Based on the nifH Gene
A total of 65 isolates and seven closely similar type strains were phylogenetically characterized based on 718 nt long DNA fragments from the nifH region. As shown in Figure 3a, the isolates were classified into three groups. Group A (nifH) comprised the majority of the isolates (59/65) positioned in the same branch as the reference strains of B. diazoefficiens (Figure 3b). Isolates of group B (nifH) were spotted in the same clade as the strains of B. yuanmingense. Four isolates from group C (nifH) resembled the B. elkanii strains.
Two isolates of group E (16S) and seven isolates of group F (16S) were observed in the same branch as the genus Stenotrophomonas of the Gamma-proteobacteria class and genus Pandoraea of the Beta-proteobacteria class, respectively. Group G (16S), comprising only one isolate, was observed to reside in the same branch as Leifsonia lichenia of the Actinobacteria class. Three isolates of group H (16S) belonged to the genus Bacillus and Firmicutes.  A total of 65 isolates and seven closely similar type strains were phylogenetically characterized based on 718 nt long DNA fragments from the nifH region. As shown in Figure 3a, the isolates were classified into three groups. Group A (nifH) comprised the majority of the isolates (59/65) positioned in the same branch as the reference strains of B. diazoefficiens (Figure 3b). Isolates of group B (nifH) were spotted in the same clade as the strains of B. yuanmingense. Four isolates from group C (nifH) resembled the B. elkanii strains.

Phylogenetic Analysis Based on the nodD1 Gene
A total of 65 isolates and eight reference type strains were phylogenetically classified based on 353 nt long DNA fragments from the nodD1 region. As shown in Figure 4, four groups were identified. Group A (nodD) consisted of 37 isolates clustered in the same branch as the B. diazoefficiens strains, with B. liaoningense in the neighboring branch, and 21 isolates of group B (nodD) were spotted in the same branch as B. ottawaense. Group C (nodD) and group D (nodD), containing one and four isolates, were formed with type strains of B. yuanmingense and B. elkanii, respectively, which was almost similar to the

Phylogenetic Analysis Based on the nodD1 Gene
A total of 65 isolates and eight reference type strains were phylogenetically classified based on 353 nt long DNA fragments from the nodD1 region. As shown in Figure 4, four groups were identified. Group A (nodD) consisted of 37 isolates clustered in the same branch as the B. diazoefficiens strains, with B. liaoningense in the neighboring branch, and 21 isolates of group B (nodD) were spotted in the same branch as B. ottawaense. Group C (nodD) and group D (nodD), containing one and four isolates, were formed with type strains of B. yuanmingense and B. elkanii, respectively, which was almost similar to the results of previous phylogenetic trees (Figures 2 and 3). One of the isolates, Tan-P1-B2-St1-28, which showed similarity with B. yuanmingense in the 16S rRNA-and nifH-based trees, and Din-P2-M1-M1-25 isolate, a sequence similarity of~83% (100% query coverage) with B. diazoefficiens, were found in a separate branch, with B. diazoefficiens, B. ottawaense, and B. liaoningense strains in the neighboring branches. results of previous phylogenetic trees ( Figure 2 and Figure 3). One of the isolates, Tan-P1-B2-St1-28, which showed similarity with B. yuanmingense in the 16S rRNA-and nifH-based trees, and Din-P2-M1-M1-25 isolate, a sequence similarity of ~83% (100% query coverage) with B. diazoefficiens, were found in a separate branch, with B. diazoefficiens, B. ottawaense, and B. liaoningense strains in the neighboring branches.

Selection of Isolates for Further Genetic and Physiological Characterization
From the 16S rRNA sequence analysis, it was found that, among the 84 isolates, 71 were from Alpha-proteobacteria, seven from Beta-proteobacteria, two from Gamma-proteobacteria, three from Firmicutes, and one from Actinobacteria ( Figure S1). Of the 84 isolates from 11 districts/locations in Bangladesh, 13 representative Bradyrhizobium isolates were selected based on location and the phylogenetic analysis of the 16S rRNA, nifH, and nodD1 genes for pot experiments, as well as further genetic and physiological characterization.

Phylogenetic Analysis of the rrn Operon (16S-ITS-23S)
Phylogenetic analysis of the partial rrn operon (16S-ITS-23S) was performed by concatenating the almost complete 16S rRNA gene, full-length 16S-23S internal transcribed spacer (ITS) region, and partial 23S rRNA gene sequences of the selected isolates ( Figure  5). The length of the sequences of the isolates was 2610-2685 nt. The selected isolates and

Selection of Isolates for Further Genetic and Physiological Characterization
From the 16S rRNA sequence analysis, it was found that, among the 84 isolates, 71 were from Alpha-proteobacteria, seven from Beta-proteobacteria, two from Gamma-proteobacteria, three from Firmicutes, and one from Actinobacteria ( Figure S1). Of the 84 isolates from 11 districts/locations in Bangladesh, 13 representative Bradyrhizobium isolates were selected based on location and the phylogenetic analysis of the 16S rRNA, nifH, and nodD1 genes for pot experiments, as well as further genetic and physiological characterization.

Phylogenetic Analysis of the rrn Operon (16S-ITS-23S)
Phylogenetic analysis of the partial rrn operon (16S-ITS-23S) was performed by concatenating the almost complete 16S rRNA gene, full-length 16S-23S internal transcribed spacer (ITS) region, and partial 23S rRNA gene sequences of the selected isolates ( Figure 5). The length of the sequences of the isolates was 2610-2685 nt. The selected isolates and seven reference type strains were phylogenetically characterized using the maximum likelihood method, and four groups were obtained. Group A (16S-ITS-23S) comprised seven isolates that branched with B. diazoefficiens. Three isolates from group B (16S-ITS-23S) were closely related to B. liaoningense, and one isolate from group C (16S-ITS-23S) was closely related to Bradyrhizobium arachidis and Bradyrhizobium guangxiense. Two isolates clustered with B. elkanii in Group D (16S-ITS-23S). seven reference type strains were phylogenetically characterized using the maximum likelihood method, and four groups were obtained. Group A (16S-ITS-23S) comprised seven isolates that branched with B. diazoefficiens. Three isolates from group B (16S-ITS-23S) were closely related to B. liaoningense, and one isolate from group C (16S-ITS-23S) was closely related to Bradyrhizobium arachidis and Bradyrhizobium guangxiense. Two isolates clustered with B. elkanii in Group D (16S-ITS-23S).

Multi-Locus Sequence Analysis of the Housekeeping Genes
Multi-locus sequence analysis (MLSA) was performed by concatenating the almost complete atpD, glnII, and gyrB genes (total of 2559 nt length of sequence of each isolate) of selected Bradyrhizobium isolates with reference strains using the maximum likelihood method ( Figure 6). The selected isolates and nine reference strains were phylogenetically characterized using the maximum likelihood method, and four groups were obtained.

Multi-Locus Sequence Analysis of the Housekeeping Genes
Multi-locus sequence analysis (MLSA) was performed by concatenating the almost complete atpD, glnII, and gyrB genes (total of 2559 nt length of sequence of each isolate) of selected Bradyrhizobium isolates with reference strains using the maximum likelihood method ( Figure 6). The selected isolates and nine reference strains were phylogenetically characterized using the maximum likelihood method, and four groups were obtained.
Group A (MLSA) comprised of seven isolates that branched with B. diazoefficiens. Three isolates from Group B (MLSA) and one isolate from Group C (MLSA) were closely related to B. liaoningense and B. yuanmingense, respectively. Two isolates clustered with B. elkanii in Group D (MLSA).

Symbiotic Performances
To assess the symbiotic ability of the 13 selected isolates, pot experiments were performed under controlled laboratory conditions using two soybean varieties, Glycine max cv. Enrei and Binasoybean-3. Different parameters, such as shoot length, main root length, nodule color ( Figure S2), nodule size and number, shoot weight, root weight, nodule weight, and ARA, were measured.
The data of nodule number of different sizes, (e.g., medium size (2-5 mm) and small size (<2 mm)) of both varieties inoculated with bacteria are presented in Table 2. Isolates Bho-P2-B2-S1-51 and Mym-P2-M3-S1-45 produced the highest number of medium-and small-sized nodules in Enrei plants, respectively. Similar to the Enrei variety, in the case of Binasoybean-3, Bho-P2-B2-S1-51 produced the highest number of medium-sized nodules, whereas Din-P2-M1-M1-25 produced the highest number of small-sized nodules. The total number of nodules produced was higher in the Binasoybean-3 variety than that in the Enrei variety.

Symbiotic Performances
To assess the symbiotic ability of the 13 selected isolates, pot experiments were performed under controlled laboratory conditions using two soybean varieties, Glycine max cv. Enrei and Binasoybean-3. Different parameters, such as shoot length, main root length, nodule color ( Figure S2), nodule size and number, shoot weight, root weight, nodule weight, and ARA, were measured.
In the Binasoybean-3 variety, most of the isolates (10/13) enhanced shoot growth in inoculated plants compared to B. diazoefficiens USDA110. However, among them, only the Tha-P2-B1-S1-68 isolate showed significantly higher shoot length than B. diazoefficiens USDA110 ( Figure S3b). In the case of root length, no significant increase was observed in comparison with B. diazoefficiens USDA110, and Bog-P3-B1-S1-29 produced the highest root length among all bacteria.
The data on shoot dry weight (DW), root DW, and nodule DW of soybean plants inoculated with isolates and B. diazoefficiens USDA110 are presented in the bar chart in Figure S4. With regard to the soybean Enrei variety, most of the plants inoculated with isolated bacteria accumulated similar or higher amounts of shoot dry weight in comparison to B. diazoefficiens USDA110. However, in the case of root DW and nodule DW, only a few isolates produced higher amounts ( Figure S4a). None of the isolates showed significantly higher shoot, root, or nodule DW gain than B. diazoefficiens USDA110. Din-P2-M1-M1-25, Bog-P3-B1-S1-29, and Nat-P3-M1-S1-79 isolates stimulated the highest shoot DW, root DW, and nodule DW, respectively.
The shoot DW, root DW, and nodule DW of the Binasoybean-3 variety measured from pot experiments are displayed in Figure S4b. Compared to B. diazoefficiens USDA110, a slight increase in shoot DW yield was observed in plants inoculated with some isolates (6/13). Pan-P1-B1-S1-69 isolate stimulated the highest shoot DW. In contrast to the Enrei variety, more isolates induced a higher root and nodule DW in the Binasoybean-3 variety and the Nil-P2-B1-M1-36 isolate produced a significantly higher amount of nodule DW in comparison to B. diazoefficiens USDA110 ( Figure S4b). This isolate also produced the highest root DW.
The results of biomass production (shoot + root + nodule DW) and ARA activity in the Enrei variety and Binasoybean-3 inoculated with isolated bacteria and B. diazoefficiens USDA110 are shown in Table 3. In the case of the Enrei variety, no substantial difference in biomass production was observed between plants inoculated with isolates and B. diazoefficiens USDA110. Some isolates produced higher amounts of plant biomass in comparison with B. diazoefficiens USDA110, and Din-P2-M1-M1-25 isolate produced the highest plant biomass among isolates, which was 7.4% higher than B. diazoefficiens USDA110. Regarding ARA activity, 8 out of 13 isolates demonstrated higher ARA activity than B. diazoefficiens USDA110, and Bho-P2-B2-S1-51 showed significantly higher ARA activity than B. diazoefficiens USDA110. Table 3. Biomass dry weight (DW) and ARA activity in soybean Enrei variety and Binasoybean-3 variety inoculated with isolated bacteria and B. diazoefficiens USDA110. All experiments were performed in triplicate and the data are expressed as the mean ± STDEV.

Discussion
The principal objective of this study was to achieve a fundamental understanding of the diversity, physiology, and symbiotic characteristics of nitrogen-fixing bacteria associated with soybean root nodules in Bangladeshi soil samples. To accomplish these goals, experiments related to bacterial isolation from root nodules, the phylogenetic characterization of bacterial strains with the 16S rRNA, rrn operon (16S-ITS-23S), nifH, nodD1, atpD, glnII, and gyrB genes, stress-tolerance tests, and plant tests in pots were performed.
Soil samples were collected from 11 different districts in Bangladesh, focusing on areas where soybean had been cultivated. Soybean plants were grown in pots using suspensions of the soil samples, and 84 isolates were isolated from the nodules of the soybean plants. Phylogenetic analysis of 84 soybean root-nodule isolates using the 16S rRNA gene revealed that the majority of the isolates (65/84) showed a phylogenetic resemblance with the genus Bradyrhizobium ( Figure 2). Among these isolates, 45 showed a close relationship with B. diazoefficiens and B. japonicum, 16 with B. liaoningense and B. yuanmingense, and four with B. elkanii. The prevalence of Bradyrhizobium strains, especially B. diazoefficiens and B. japonicum, in root nodules of soybean is a common scenario and is scientifically expected.
B. diazoefficiens (and/or B. japonicum) and B. elkanii species have been found in diverse climatic regions across the world, and strains of B. liaoningense are abundantly found in alkaline soils of temperate to subtropical climates in South and Southeast Asia, while the warm tropical climates of India and Nepal support B. yuanmingense [32]. Furthermore, in a research expedition for the genetic analysis of 76 bacteria isolated from root nodules of soybean inoculated with soil samples from Myanmar, India, Nepal, and Vietnam, Vinuesa et al. (2008) showed that 75 of the isolates were phylogenetically similar to B. japonicum (USDA110) or currently, B. diazoefficiens (in Nepal), B. liaoningense (in Myanmar, Vietnam), B. yuanmingense (in India, Myanmar, and Vietnam), and B. elkanii (in Myanmar) [28]. Another study on root-nodule isolates of soybean from five regions in India, performed by Appunu et al. (2008), showed that 36% of isolates were identified as B. yuanmingense, 26% as B. liaoningense, and 38% were different from all described Bradyrhizobium species but showed the same symbiotic genotype as B. liaoningense and B. japonicum [33].
Besides these slow-growing Bradyrhizobium isolates of the Alpha-proteobacteria class, some other bacterial strains, such as Methylobacterium of Alpha-proteobacteria, Pandoraea of Beta-proteobacteria, Stenotrophomonas of Gamma-proteobacteria, Bacillus of Firmicutes, and Leifsonia of Actinobacteria phylum/class, were also found in this study. Methylobacterium of the Methylobacteriaceae family and Rhizobiales order, belonging to the class Alpha-proteobacteria, is a well-established group of bacteria associated with the root nodules of leguminous plants [34][35][36][37]. Bacterial strains of the genus Burkholderia of the Burkholderiaceae family of Beta-proteobacteria have also been isolated by many researchers from root nodules [6,8,38,39]. Therefore, the isolation of bacteria belonging to the Pandoraea genus of the Burkholderiaceae family of the Beta-proteobacteria class is new but not unnatural. Stenotrophomonas bacteria of the Gamma-proteobacteria class have been isolated from the rhizosphere soil of leguminous plants [40]. Therefore, the existence of bacteria of this genus in root-nodules is also reasonable. Different species of Bacillus are commonly found in the root nodules of different leguminous plants and are thus treated as a well-known phenomenon [24,25,41]. In addition, many bacterial strains belonging to the phylum Actinobacteria, such as Frankia, Rhodococcus, Arthrobacter, Brevibacterium, Micromonospora, and Streptomyces, have also been isolated from root nodules and characterized for their nitrogen-fixing or plant-growthpromoting capabilities [3,[42][43][44][45][46].
According to researchers, the analysis of symbiotic genes, such as nifH, nifD, nifK, nodA, nodB, nodC, and nodD, of root nodule bacteria can provide valuable insights into their host range and symbiotic relationships, which may vary from those expected by rRNA-based phylogenies [6,20,38,47]. In the nifH phylogenetic analysis, 59 out of 65 Bradyrhizobium isolates possessed the nifH gene, similar to that of B. diazoefficiens (Figure 3). Except for the four isolates that showed 16S rRNA gene resemblance with B. elkanii (group C (16S) of Figure 2), and two isolates out of 16 of Bradyrhizobium liaoningense and Bradyrhizobium yuanmingense (group B (16S) of Figure 2c), all other isolates showed nifH gene similarity to B. diazoefficiens. In an analysis of a database consisting of 32,954 aligned sequences of the nifH gene, Gaby et al. (2014) revealed that genomes that had >97% similarity in the 16S rRNA gene could possess up to 23% dissimilarity in the nifH sequences [48]. Therefore, B. yuanmingense and B. liaoningense, which have a 16S rRNA gene similarity of approximately 99% with B. diazoefficiens, can logically possess the nifH gene resembling B. diazoefficiens, which has approximately 10% sequence dissimilarity. This transfer of the nifH gene may be the outcome of vertical gene transfer in the symbiotic region of Bradyrhizobium species. In the case of nodD1 gene phylogeny, almost similar results were obtained for the nifH gene tree, with some discrepancies. Most isolates (n = 37) clustered in the same clade as B. diazoefficiens (in group A (nodD) of Figure 4b). The Tan-P1-B2-St1-28 and Din-P2-M1-M1-25 isolates showed a sequence similarity of~83% (100% query coverage) with B. diazoefficiens and were placed in a different branch adjacent to Bradyrhizobium branches. No change was observed in the four isolates that showed similarity with B. elkanii in the nodD1 gene tree or the 16S rRNA and nifH tree, suggesting a strong phylogenetic relation with B. elkanii. Therefore, further analysis should be performed to verify the symbiotic gene transfer process.
Some researchers have shown that the analysis of the rrn operon (16S-ITS-23S) signifies more than the analysis of the 16S rRNA gene for the resolution of taxa at the species level [49,50]. Phylogenetic analysis of the partial rrn operon consisting of the almost complete 16S rRNA gene, full-length 16S-23S internal transcribed spacer (ITS) region, and partial 23S rRNA gene sequences for the selected isolates was congruent with the 16S rRNA gene-based tree ( Figure 5). MLSA was also performed by concatenating the atpD, glnII, and gyrB genes of selected Bradyrhizobium isolates with type strains using the maximum likelihood method (Figure 6). Seven of the 13 isolates branched with B. diazoefficiens, three branched with B. liaoningense, and one branched with B. yuanmingense. Two isolates clustered with B. elkanii. For all these isolates, the MLSA position could be considered the true taxonomic identity. The use of MLSA of the housekeeping genes in bacterial species definition and identification has been strongly endorsed by the scientific community [51]. Unless whole-genome sequencing is performed, MLSA is considered a distinctive identification process for exploring the genetic diversity within a proposed species [51]. Therefore, based on 16S rRNA gene, MLSA and rrn operon phylogenetic analysis, except for the Tan Lak-P1-S1-M1-46, B. liaoningense sp. Bho-P2-B2-S1-51, B. liaoningense sp. Lak-P1-M1-S1-85, B. liaoningense sp. Pan-P1-B1-S1-69 for taxonomic identity. For Tan-P1-B2-S1N-84 isolate further analysis is needed to confirm the taxonomic identity.
In this study, among the different phylogenetic methods, including Bayesian inference, maximum-likelihood, maximum-parsimony, and neighbour-joining, the maximumlikelihood method was applied for phylogenetic reconstruction. Although each of these techniques has its strengths and weaknesses, maximum-likelihood (ML) and Bayesian methods typically provide more realistic and robust phylogenies because they employ explicit models of molecular evolution [51]. The analysis of the highest DNA sequence similarity of the selected isolates with Bradyrhizobium type strains showed no variation in the case of Bradyrhizobium diazoefficiens isolates for all analyzed genes. However, some discrepancies were observed for isolates of other genera (Table 5). Among the housekeeping genes, gyrB was found to be more decisive in identifying the isolates than atpD and glnII, and showed similarities with rrn operon sequence data.
To evaluate the symbiotic ability of the 13 selected isolates, pot experiments were performed using two soybean varieties, Glycine max cv. Enrei and Binasoybean-3, and different parameters were measured. All isolates formed red or pink (cross-sections) nodules, indicating the presence of leghemoglobin in the soybean Enrei variety and Binasoybean-3 variety, except Tan-P1-B2-S1N-84 belonging to B. yuanmingense, which did not form any nodules in the Enrei variety ( Figure S2, Table 2). The variation in nodule production could be due to the variation in the agronomic traits of the two soybean cultivars and differences in the compatibility of the isolates with the cultivars. Wang et al. (2019) described that soybean cultivars can influence nodule formation, and some bacteria can form nodules in some varieties, while others cannot [3]. The number of isolates that formed nodules and the number of nodules formed were higher in the Binasoybean-3 variety than those in the Enrei variety (Table 2). This could be due to the fact that Binasoybean-3 is a local Bangladeshi variety, meaning that isolates of Bangladeshi soil samples have a higher affinity towards it. The performance of the isolates varied in terms of effective symbiosis and biomass production. The nodulation and nitrogen fixation depends on many aspects, especially the symbiotic genes (such as such as nifH, nifD, nifK, nodA, nodB, nodC, and nod, etc.), which can vary from rhizobia to rhizobia; thus, different rhizobial species show different host-specificity and nitrogen fixation performances [3,7]. When inoculated, some of the isolates enhanced biomass production in both the Enrei and Binasoybean-3 varieties in comparison to B. diazoefficiens USDA110, whereas some decreased. In the case of Enrei, the Din-P2-M1-M1-25 isolate produced the highest plant biomass among isolates, which was 7.4% higher than B. diazoefficiens USDA110. For Binasoybean-3, Nil-P2-B1-M1-36 isolate produced the highest plant biomass among all the isolates, which was 14% higher than B. diazoefficiens USDA110. In comparison to B. diazoefficiens USDA110, three isolates produced 10% or higher biomass gain in Binasoybean-3 variety, compared to no isolate in the case of the Enrei variety. These results suggest that the isolates were more compatible with Bangladeshi local varieties. Having a phylogenetic similarity with B. liaoningense, Bho-P2-B2-S1-51 was the only isolate that produced a significant increase in ARA compared to B. diazoefficiens USDA110 in both the Enrei and Binasoybean-3 varieties, increasing plant biomass production by 9% in case of Binasoybean-3 variety. Therefore, this isolate could be a potential biofertilizer candidate. Regarding Binasoybean-3, the Tha-P2-B1-S1-68 isolate belonging to B. diazoefficiens significantly increased shoot length compared to B. diazoefficiens USDA110 ( Figure S3) and enhanced biomass production by 10% (Table 3), making the isolate a suitable candidate as biofertilizer for local varieties. According to Chen et al. (2021), the highly effective symbiosis could only be obtained when the host specificity and habitat specificity cooperated in a rhizobial symbiont; therefore, the rhizobia and legume variety of same locality can demonstrate better symbiotic effects [7].The relationships between different parameters of pot experiments, e.g., nodule number, small and medium size nodule, nodule dry weight, ARA activity, and biomass, were observed to be nonlinear and varied from isolate to isolate (Table 6). Among the different parameters, it was observed that, for both soybean varieties, nodule dry weight had a significant positive correlation with ARA (in µmol/h/plant) and biomass DW, suggesting that with an increase in nodule dry weight, ARA (in µmol/h/plant) and biomass production may increase. It is known that nodule dry weight is a good indicator of symbiotic efficacy, and an important parameter in strain evaluation [52]. Allito et al. also observed similar correlation and mentioned nodule dry weight as one of the key indicators of effective legume-rhizobia symbiosis [53]. Although the medium-size nodule number showed a significantly high positive correlation with biomass DW and ARA (in µmol/h/plant) in Binasoybean-3, small nodule numbers displayed insignificant and/or negative correlations, suggesting that a larger nodule size can significantly contribute to biomass production and ARA (in µmol/h/plant) than immature/smaller ones. No significant positive correlations were observed with other parameters. Similar to these results, some researchers did not find any clear relationship between ARA activity, nodule numbers and plant biomass in their studies [25,54].
In the stress tolerance tests, all the selected isolates grew at 1-4% NaCl concentration. Similarly, most isolates grew at pH 4.5, except four isolates, and almost all isolates grew at pH 4.5-10, which demonstrate that the isolates are tolerant to high salt and pH conditions. Under saline conditions, salt-tolerant rhizobia are beneficial for host plant growth. For instance, the growth of alfalfa is enhanced when inoculated with the salt-tolerant Rhizobium strain Rm1521 rather than the salt-sensitive strain A2 under NaCl treatment [55]. Among the isolates, five isolates, namely, Bho-P2-B2-S1-51, Bog-P3-B1-S1-29, Mym-P2-M3-S1-45, Nat-P3-M1-S1-79, and Tha-P2-B1-S1-68, survived at 45 • C, demonstrating their suitability as biofertilizer candidates for tropical environments. Chen et al. (2002) isolated Rhizobia from soybeans that could grow at 40 • C [56]. Correspondingly, Phaseolus vulgaris nodule bacteria were found to survive at 47 • C [57], and Rhizobium isolated from Sesbania aculeata plants were observed to grow at 50 • C [58,59]. As temperature, salinity, origin of cultivar, etc., are prominent effectors of rhizobial symbiosis [19], the isolates with temperature and salinity tolerance could be potential biofertilizers, especially for local legume varieties. * Denotes significance at 95% confidence level. The correlation analysis and pot experiment results indicated that the isolated bacteria have showed better activity in the Bangladeshi local variety, Binasoybean-3 than the Japanese variety, Enrei. The reason for this could be the higher symbiotic affinity, compatibility and effectiveness between isolates and plant variety of same origin-Bangladesh.

Conclusions
This study represents the first genetic and physiological characterization of soybeannodulating isolates from different regions of Bangladesh, one of the top soybean oil consumer countries of the world. Surprisingly, diverse arrays of both nodule-forming and endosymbiotic bacteria belonging to Bradyrhizobium and Methylobacterium of Alphaproteobacteria, Pandoraea of Beta-proteobacteria, Stenotrophomonas of Gamma-proteobacteria, Bacillus of Firmicutes, and Leifsonia of Actinobacteria phylum/class were found to be present in the root nodules of soybean plants grown in pots using a suspension of soil samples from Bangladesh, with some potential isolates showing a higher symbiotic ability than B. diazoefficiens USDA110 in pot experiments. One of the selected isolates, Bho-P2-B2-S1-51, which showed a phylogenetic similarity to B. liaoningense, produced a significant increase in ARA compared to B. diazoefficiens USDA110 in both the Enrei and Binasoybean-3 varieties, induced plant biomass production by 9% in case of Binasoybean-3 variety, and survived at 45 • C. Therefore, this isolate could be a potential biofertilizer candidate. With regard to Binasoybean-3, the Tha-P2-B1-S1-68 isolate, which showed a phylogenetic similarity to B. diazoefficiens, also significantly increased shoot length compared to B. diazoefficiens USDA110, enhanced biomass production by 10%, and could grow at 45 • C, making this isolate a suitable candidate for use as biofertilizer for local Bangladeshi soybean varieties. Further studies, such as whole-genome sequence analysis and field trials, will need to be performed for prospective candidate isolates as biofertilizers to acquire an in-depth understanding of the genetic make-up and functional characteristics of soybean-rhizobium symbiosis.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms10112282/s1, Table S1: Information of primers used in this study.; Table S2: Basic information of soybean varieties used in this study.; Table S3: List of selected Bradyrhizobium isolates for pot experiments and further characterization.; Figure S1: Doughnut graph (a) and map of Bangladesh (b) showing the distribution of different types of bacteria.; Figure S2. Pictures of dissected nodules of soybean Enrei variety and Binasoybean-3 variety inoculated with selected isolates and B. diazoefficiens USDA110.; Figure S3: Shoot and root length of soybean Enrei variety (a) and Binasoybean-3 variety (b) inoculated with selected isolates and B. diazoefficiens USDA110. Data are presented in average value with error bar denoting the STDEV of three plants each. '*' Denotes significance with "B. diazoefficiens USDA110" at 95% confidence level using Dunnett's test.; Figure S4: Shoot dry weight (DW), root DW and nodule DW of soybean Enrei variety (a) and Binasoybean-3 variety (b) inoculated with isolates and B. diazoefficiens USDA110. Data are presented in average value with error bar denotes STDEV of three plants each. '*' Denotes significance with "B. diazoefficiens USDA110" at 95% confidence level using Dunnett's test.