Phylogenetic evidence of allopatric speciation of bradyrhizobia nodulating cowpea (Vigna unguiculata L. walp) in South African and Mozambican soils

ABSTRACT The legume host and soil environment play a major role in establishing effective symbiosis with diverse rhizobia for plant growth promotion and nodule formation. The aim of this study was to assess the morpho-physiology, distribution and phylogenetic position of rhizobia nodulating cowpea from South Africa and Mozambique. The results showed that the isolates were highly diverse in their appearance on yeast mannitol agar plates. The isolates tested also showed an ability to produce IAA at concentrations ranging from 0.64 to 56.46 μg.ml−1 and to solubilise phosphorus at levels from 0 to 3.55 index. Canonical correspondence analysis showed that soil pH and mineral nutrients significantly influenced bradyrhizobial distribution. Analysis of BOX-PCR placed the isolates in eight major clusters with 0.01 to 1.00 similarity coefficient which resulted in 45 unique BOX-types. Phylogenetic analyses based on 16S rRNA, atpD, glnII, gyrB and recA gene sequences showed distinct novel evolutionary lineages within the genus Bradyrhizobium, with some of them being closely related to Bradyrhizobium kavangense, B. subterraneum and B. pachyrhizi. Furthermore, symbiotic gene phylogenies suggested that the isolates’ sym loci probably relates to the isolates’ geographical origin. The results indicated that geographical origin did affect the isolates’ phylogenetic placement and could be the basis for allopatric speciation


INTRODUCTION
Legumes are a major component of agricultural systems throughout the world due to their N 2 -fixing ability when in symbiosis with soil bacteria called rhizobia. This makes nodulated legumes good candidates for improving soil fertility, especially nitrogen (N) (Mohale, Belane and Dakora 2013). On a global scale, symbiotic bacteria contribute the largest amount of N through BNF to terrestrial ecosystems, and are thus vital for agricultural production (Unkovich and Pate 2000;Echevarria-Zomeño et al. 2009). N 2 fixation is usually performed by the bacteria in root nodules, and the NH 3 produced is exchanged with the host plant for photosynthate (Prell and Poole 2006). The distinctive feature of the legume symbiosis is the formation of root or stem nodules following a series of molecular exchanges between the microsymbiont and the host plant (Chidebe, Jaiswal and Dakora 2018).
Cowpea is one of the most important indigenous legume crops in Africa, with both wild and cultivated forms growing across the continent (Mohammed, Jaiswal and Dakora 2018). West Africa is the centre of origin of this crop Puozaa, Jaiswal and Dakora 2019). Cowpea can meet 66% to 96% of its N nutrition from symbiotic fixation .
The high cost of N fertilizers and the environmental pollution associated with their frequent use to increase grain legume production can be overcome through the use of highly effective rhizobia that can efficiently meet the crop's N demand in a sustainable and environmentally friendly manner. In Africa, the most common factor limiting legume productivity is low availability of mineral nutrients, especially N and phosphorus (P) (Dakora and Keya 1997). Plants tend to perform better in the presence of sufficient numbers of effective soil rhizobia that can fix atmospheric N 2 , solubilize soil P, and produce plant growthpromoting compounds such as indole acetic acid (IAA) which promote root growth and water nutrient uptake.
Plant growth-promoting bacteria can therefore influence plant development directly by synthesizing root-growth hormones thus facilitating the supply and uptake of mineral nutrients from the soil. Solubilization of P compounds in soil by rhizobia is an important strategy for enhancing plant growth and productivity (Marra et al. 2012). Although P is an important macronutrient for plant growth, about 95% to 99% of soil P occurs in the insoluble form that cannot be utilized by plants (Vassilev et al. 2001). Some rhizobial bacteria are however capable of solubilizing unavailable soil P complexes for plant uptake and development, just as some rhizobia are also capable of producing IAA, a hormone that is involved in root formation and root elongation for increased uptake of water and nutrients (Etesami et al. 2009). Thus, the identification of IAA-producing and P-solubilizing traits among N 2 -fixing rhizobia is important for selecting rhizobia for inoculant production.
A robust rhizobial taxonomy and phylogeny is key to identifying microsymbionts for improving legume productivity using biological N 2 fixation (Chidebe, Jaiswal and Dakora 2018). Such studies have increased our understanding of the genetic diversity of rhizobia and significantly improved rhizobial classification through the re-definition of species. Previous reports from the different parts of the world such as Botswana, Ghana, South Africa, Senegal, Mozambique, Ethiopia, Greece, Spain and Brazil have identified Bradyrhizobium as the main microsymbiont nodulating cowpea (Pule-meulenberg et al. 2010;Bejarano et al. 2014;Wade et al. 2014;Degefu et al. 2017;Marinho et al. 2017;Tampakaki et al. 2017;Chidebe, Jaiswal and Dakora 2018;Mohammed, Jaiswal and Dakora 2018;Jaiswal and Dakora 2019;Puozaa, Jaiswal and Dakora 2019). However, knowledge gaps still exist in the genetic diversity and biogeographic distribution of indigenous N 2 -fixing rhizobia in Southern Africa. Since cowpea is an indigenous legume crop of the vast and agroecologically divergent African continent, identification and characterization of native rhizobia nodulating cowpea cultivated in different geographic regions of Africa will help to understand and unravel the factors shaping rhizobial diversity and distribution in African soils. Therefore, the aim of this study was to describe the morpho-physiological and molecular characteristics of rhizobial isolates nodulating cowpea from different ecological environments as well as to assess these bacteria for P solubilisation and IAA production.

Field experimentation, nodule collection and isolation
Nodules were collected from cowpea genotypes grown in the field at Marapyane in Nkangala District, west of the Mpumalanga Province of South Africa and at Ruace and Muriaze of Mozambique. Marapyane is located at the geographical coordinates of 24 • 57 60 S, 28 • 45 76 E and elevation 1022 meters above sea level, while Ruace and Muriaze are located at 15 o 14 17.5 S, 26 o 43 44.8 E and 15 o 09 12.9 S, 30 o 19 20 E geographical coordinates, respectively (Fig. 1). The test sites had been under cultivation of various cereals and grain legumes for the past five years. Before planting cowpea, soil samples were collected at a depth of 0-30 cm for physico-chemical analysis (Table S1, Supporting Information). The experiment was laid out as a randomised complete block design with three replications. Six and five cowpea genotypes were used to trap microsymbionts in the fields of South Africa and Mozambique, respectively (Table 1). There were three common cowpea genotypes (IT97K-1069-6, IT00K-126-3, and IT18) grown in the test locations of both South Africa and Mozambique (Table 1). Each experimental plot measured 1.2 × 3 m 2 . There were four rows in each plot, with 1.2 m row length. Seeds were planted with an inter-row spacing of 75 cm and intrarow spacing of 15 cm.
At 50 days after planting, plants were carefully dug out from the Marapyane site (South Africa) and from two farmers' fields at Ruace and Muriaze (Mozambique) (Fig. 1), and transported to the laboratory. The roots with attached nodules were separated and washed in running tap water to remove adhered soil. Healthy, unbroken and firm root nodules were selected and detached from the root using forceps and scissors and stored in vials containing silica gel at 4 • C until bacterial isolation. During isolation, the nodules were hydrated in water, and surface-sterilized by immersion in 75% ethanol for 1 minute, then washed in 3.5% NaOCl for 3 minutes, followed by rinsing five times with sterile distilled water. Each surface-sterilized nodule was then crushed in a loop of sterile distilled water in a sterile petri dish, and the nodule macerate streaked on yeast mannitol agar (YMA) plates and incubated at 28 • C. Colony appearance was observed from two to twelve days after incubation, and colony characteristics such as diameter, shape, colour and texture were recorded on YMA plates after 15 days of incubation for each isolate.

Nodulation test with isolates
To fulfil Koch's postulates, the single-colony bacterial isolates were examined for their ability to form root nodules on cowpea, their homologous host plant. Healthy cowpea seeds were surface-sterilized by immersion in 75% alcohol for 10 seconds. The alcohol was drained, and the seeds placed in 3% sodium hypochlorite for 1 to 2 minutes. The seeds were then rinsed six times with sterile distilled water and two seeds planted in sterile autoclaved sand contained in sterile pots (5006.05 cm 3 ) using a randomised complete block design with three replications under glasshouse conditions. After germination, the seedlings were thinned to one plant per pot. Six-day-old seedlings were inoculated with 1 mL exponential phase-grown bacterial suspension (∼ 10 6 -10 7 cells/ml) using sterile micropipettes. Uninoculated, plants were included as control. The test seedlings were supplied with sterile N-free nutrient solution (Broughton and Dilworth 1971) twice a week and sterile distilled water when necessary. At 8 weeks after planting, the plants were harvested, and nodulation data recorded. Deep green leaf colour and pink nodules were considered as an indication of effective nodulation.
The ability of rhizobial isolates to solubilize P was assessed by the Double-layered plate method (Yanni et al. 2001). At the exponential growth stage, 10 μL rhizobial culture was spotted onto double-layered agar medium plate and incubated at 28 • C for 7 days. The formation of a clear halo zone around a bacterial spot was taken as indicative of P solubilisation. The diameter of each halo zone was measured to identify isolates having high P-solubilising activity. Phosphate-solubilizing capacity was calculated in terms of phosphate solubilization index (PSI) as PSI = A/B, where A is the total diameter of the halo zone, and B is the colony diameter. The isolates showing PSI ≥ 2 were considered as high phosphate-solubilizing bacteria (Rahman et al. 2014).

DNA extraction and BOX-PCR fingerprinting of bacterial isolates
Bacterial genomic DNA was extracted using Sigma's Bacterial Genomic DNA Kit according to the manufacturer's instructions (GenElute TM ), and used for rep-PCR (BOX-PCR) with a Box-A1R primer (Table S2, Supporting Information). The final volume of the PCR reaction mixture was 25 μL and contained 1 μL (50-70 ng μL −1 ) genomic DNA, 12.5 μL MyTaq PCR master-mix (2x), 1μL (10 μM) BOX-A1R primer, and 10.5 mL sterile ultrapure water. PCR amplification was done in a Thermal cycler (T100 BIORAD, USA) following standard temperature profiles (Table S2, Supporting Information). The PCR-amplified products were electrophoresed in 1.2% agarose gel for 6 hours, visualized and photographed using gel documentation system (BIO-RAD Gel Doc TM XR+). The sizes of bands were determined using the Image Lab software (Bio-Rad version 4.1). The band pattern was recorded in a binary form (1, 0) and cluster analysis was carried out with the UPGMA (Unweighted Pair Group Method with Arithmetic mean) algorithm using the NtSyspc 2.1 software (Rohlf FJ, 2009). PCR amplification of the 16S rRNA, housekeeping (atpD, glnII, gyrB, recA) and symbiotic (nifH and nodC) genes PCR amplification of 16S rRNA, protein-coding (atpD, glnII, gyrB and recA) and symbiotic (nifH and nodC) genes were each done in 25 μL reaction volume. The reaction mixture contained 1 μL DNA (50-70 ng μL −1 ), 3 μL of myTaq buffer (5x), 1 μL (10 μM) each of forward and reverse primers of the gene of interest, 0.1 μL Taq polymerase (5U) (Bioline, USA) and 18.9 μL double distilled ultrapure water. The PCR process was carried out with Bio-Rad T100 thermal cycler using standard temperature profiles ( Table  S2, Supporting Information). The amplified gene products were confirmed by gel electrophoresis in 1.2% agarose gel stained with ethidium bromide in 1X TAE buffer at 85V for 1 hour.

Sequencing and phylogenetic analysis
For sequencing, the amplified PCR products were purified using PCR Cleanup kit (NEB, USA) according to the manufacturer's instruction. The purified amplified DNA was sent to Macrogen (Netherlands) for sequencing. Thereafter, the quality of sequences was verified using the software BioEdit 7.0.9.0 (Hall 1999). The BLASTn program was used to identify closely related species in the NCBI database. Pairwise and multiple sequence alignments were done with CLUSTALW, and phylogenetic trees constructed by means of the maximum likelihood statistical method using MEGA 7 software (Kumar, Stecher and Tamura 2016). The robustness of branching was estimated using 1000 bootstrap replicates (Felsenstein 1985). The obtained sequences were deposited in the NCBI GenBank to get the accession numbers which were indicated in phylogenetic tress. The test gene sequences of the isolates used for phylogenetic tree constructions are shown in Figure S2 (Supporting Information).

Statistical analysis
The effect of soil factors on the distribution of cowpeanodulating bradyrhizobia was examined using canonical correspondence analysis (CCA) with vegan (version 2.4-2) (Oksanen et al. 2010) of R software (R Core Team 2015). Analysis was done with only the soil parameters which showed significant ecological contribution to the distribution of the test cowpeanodulating bradyrhizobia. The general permutation test was used to assess the statistical significance of the ordination axes

Plant nodulation data
A total of 102 isolates were obtained from the cowpea nodules collected from all three study sites in South Africa and Mozambique. Out of the102 isolates, 61 could elicit nodule formation in cowpea (their homologous host) under glasshouse conditions. Isolates that produced green leaves and pink nodules were considered as effective rhizobia. Only effective nodule-forming isolates were used for further studies.

IAA production and phosphate-solubilising activities by authenticated isolates
The test rhizobial isolates varied in their ability to produce IAA, with concentrations ranging from 0.64 μg ml −1 to 56.46 μg ml −1 ( Table 1). The highest IAA producer was isolate TUTVUMm119 (56.46 μg ml −1 ) and the least isolate TUTVUMr93 (0.64 μg ml −1 ) ( Table 1). Most of the isolates tested showed their ability to solubilize phosphate (Ca 3 PO 4 ). The halo zone produced by P solubilization appeared within 3 days after spotting culture on the plate. Out of 53 isolates tested, 36 were able to solubilize inorganic phosphate. The phosphate solubilisation index (PSI) differed significantly (P < 0.05) among the isolates, with an index range of 0 to 4.0 ( Table 1). The highest solubilization index was produced by the Marapyane isolates, followed by Muriaze, and then Ruace. Of the isolates selected from Marapyane for this study, 90% were found to solubilize P, with 10% lacking solubilisation ability (no halo zone formation). The selected isolates from Mozambique showed 53% solubilization, with 47% not being able to solubilize P ( Table 1). The total solubilisation by the test isolates stood at 66%.

Correlation analysis
There was a significant correlation between P solubilization and IAA production. Isolates that elicited high IAA production also showed greater P solubilization index (r = 0.187 * ), suggesting that greater P accumulation supported IAA production (Fig. S3, Supporting Information).

Box-PCR fingerprinting analysis of genomic DNA
PCR amplification of the BOX region of the rhizobial genomic DNA from 54 representative isolates (including the commercial isolate TUTVUI122) resulted in distinct banding patterns that ranged from 1 to 14 fragments per DNA profile, with band sizes of 10 037 to 440 bp. Based on banding pattern results, the bacterial isolates were grouped into eight clusters (Clusters I-VIII) with 0.01 to 1.00 similarity coefficient, and yielded 45 unique BOX-fingerprints in the test rhizobial populations (Table 1; Fig. S4, Supporting Information). The rhizobial strains obtained from the three study sites were distributed across almost all the identified clusters, except for Cluster IV which consisted only of Mozambican isolates. The majority of the rhizobial isolates from the three experimental sites grouped in Cluster V. Isolates obtained from Marapyane were highly diverse and occupied seven clusters (Clusters I, II, III, V, VI, VII and VIII) while those from Muriaze and Ruace in Mozambique were distributed in five clusters within the dendrogram (Fig.  S4, Supporting Information). Based on the test cowpea genotypes, variety IT97K-1069-6 attracted highly diverse microsymbionts which occupied six clusters (Clusters I, II, III, V, VI, and VII), followed by genotypes IT00K-126-3 and IT18 which attracted diverse isolates contained in five (I, II, IV, V, and VIII) and four (III, IV, V and VI) clusters, respectively (Table 1; Fig.  S4, Supporting Information). Based on BOX-PCR fingerprinting, there was a 100% genomic similarity for isolates TUTVUMp7 and TUTVUMp12; isolates TUTUVMp26 and TUTVUMp28; isolates TUTVUMm108 and TUTVUMm109; isolates TUTVUMr101, TUTVUMm104, TUTVUMm105 and TUTVUMm107; as well as isolates TUTVUMm110 and TUTVUMm111, with all those isolates originating from the same site, except for isolate TUTVUMr101 (Fig. S4, Supporting Information).

Phylogenetic positons of the test isolates
Isolates were selected from each cluster of BOX-PCR dendrogram for further phylogenetic analysis. The PCR amplification of 16S rRNA gene yielded ∼1.5 kb amplicons of full length for the test isolates. The BLAST n 16S rRNA sequence analysis of these isolates showed that all the strains had high sequence similarities with Bradyrhizobium with 73.8% conserved and 24.4% variable sequence sites (Table S2, Supporting Information). The maximum likelihood phylogeny of the 16S rRNA gene placed the test isolates into five distinct clusters (Clusters I to V) (Fig. 2). Cluster I included isolates TUTVUMm113, TUTVUMm115, TUTVUMm116, TUTVUMp6 and TUTVUMr94 which did not group with any of the reference type strains, but B. kavangense 14-3 T and B. subterraneum 58 2-1 T were the closest related species with 99.1%-99.5% sequence identity. Cluster II comprised two sub-clusters (IIa and IIb  stylosanthis BR 446 T and B. iriomotense NBRC 102520 T . Isolates TUTVUMp8, TUTVUMp48, TUTVUMp53, TUTVUMp56, TUTVUMp65, TUTVUMr89, TUTVUMr103, TUTVUMm112 and TUTVUMm114 recorded a 100% sequence similarity with the type reference strains of B. kavangense 14-3 T , and B. subterraneum 58 2-1 T with 75% bootstrap support in Cluster V while isolate TUTVUMp5 showed outgroup of Cluster V (Fig. 2).
Sequence and phylogenetic analysis of atpD, glnII, gyrB and recA genes PCR-amplification of atpD, glnII, gyrB and recA genes, which are responsible for ATP synthase subunit beta, glutamine synthetase II, DNA gyrase subunit B and recombinase A (in that order), yielded 600, 600, 750 and 700 bp band lengths, respectively. The results of sequence analysis showed that atpD gene had the highest conserved sequences (63.2%) and gyrB gene the lowest (48.7%). The highest (33.9%) parsimony-informative sites were found in gyrB gene (Table S3, Supporting Information). Except for recA phylogeny, the topography of the individual (atpD, glnII and gyrB) housekeeping gene phylogenies were very similar and congruent with the 16S rRNA phylogeny with little inconsistency ( Fig. S5-S7, Supporting Information). For example, isolates TUTVUMm110 and TUTVUMr84 grouped with B. yuanmingense in Cluster I of atpD phylogeny, with B. arachidis in Cluster IV of glnII and gyrB phylogenies ( Fig. S5-S7, Supporting Information).
The recA phylogeny showed discrepancy compared to the other gene phylogenies. For example, isolates TUTVUMr99 and TUTVUMr82 grouped with B. yuanmingense CCBAU 05623 T and B. centrosematis A9 T in Cluster I and Cluster III, respectively (Fig. 3), while they showed proximal relation with the B. elkanii group in Cluster II of the atpD, glnII and gyrB phylogenies (Figs. 3, S5-S7, Supporting Information). Similarly, isolate TUTVUMm105 in Cluster IV aligned with isolates TUTVUMm108 and TUTVUI122 clustered with B. arachidis in recA phylogeny, as an outgroup while it grouped with isolate TUTVUMr94 in Cluster I in atpD and glnII phylogenies which had B. yuanmingense as the closest related species. In the same way, isolate TUTVUMr93 grouped with isolates of Cluster I in recA phylogeny, while it grouped with isolates of Cluster III in the other phylogenies (Figs. 3; S5-S7, Supporting Information).

Concatenated gene phylogeny
If each separate gene tree represents a partial view of the information contained in the genome, all four housekeeping genes together could lead to a better understanding of the relationships between the rhizobial strains. A concatenated gene sequence analysis was therefore performed to refine the phylogenetic positions of the test isolates within the Bradyrhizobium lineage. A maximum likelihood phylogeny was constructed using atpD+glnII+gyrB concatenated aligned sequences of the test isolates and type reference strains common across the three housekeeping genes. Due to inconsistency in recA sequences, we did not include it in the concatenated sequence and phylogenetic analysis. The concatenated consensus sequence (1421 bp) comprised 818 conserved, 603 variables, 422 parsimonies informative and 181 singleton sites (Table S3, Supporting Information). The phylograms constructed from concatenated sequences of housekeeping genes was congruent to individual gene phylograms. Like the individual gene phylogenies, the isolates occupied five clusters in the concatenated gene tree. Isolates TUTVUMr103, TUTVUMm112, TUTVUMm113, TUTVUMm114, TUTVUMm115, TUTVUMm116 and TUTVUMp6 grouped together in Cluster I with type reference strain B. yuanmingense CCBAU 10071 T with 98.8%-99% sequence similarity and 100% bootstrap support (Fig. 4), while in the same Cluster I, isolates TUTVUMm105 and TUTVUMr94 clustered together as an outgroup and had B. yuanmingense as the closest related species with a low 95.8% sequence similarity. The isolates in Cluster II belonged to the B. elkanii group formed two sub-clusters (IIa and IIb) but with a low (<97%) sequence similarity; only isolate TUTVUMp60 showed >97% sequence identity with B. pachyrhizi PAC48 T . The isolates in Cluster III, IV and V did not tightly group with any reference type strains, but they had 94%-96.8% sequence similarity with B. kavangense 14-3 T and B. arachidis CCBAU 45332 T as the closest related species (Fig. 4).

Phylogenetic analysis based on symbiotic nifH and nodC gene sequences
The PCR amplification of full-length amplicons of the symbiotic nifH and nodC genes resulted in clear bands of 650 and 300 bp, respectively, as the amplified products. We did not find nodC amplification of all selected isolates. Therefore, we could not include all selected isolates with nodC sequences in the phylogenetic analysis. The non-amplification of nodC genomic region could be due to incompatibility of used primer pairs with the isolates' genomes. In the maximum likelihood phylogeny, the isolates clustered into five groups (Figs 5-6) which were congruent with the phylogenies of the individual housekeeping genes, as well as with the phylogeny based on concatenated gene sequences (Figures 2 and 4; Figures S5-S7, Supporting Information).
The isolates in Cluster IV were closely related to B. arachidis CCBAU051107 T with a sequence identity of 97% in both nifH and nodC phylogenies, while isolate TUTVUMr84 in the same cluster was highly related to B. kavangense 14-3 T with 98.5% sequence identity in the nifH phylogeny. Isolates in Cluster V showed a close relationship with B. subterraneum 58-2-1 T with 99.5%-100% sequence similarity, but this group stood alone in the nodC phylogeny due to the absence of B. subterraneum in the tree. Isolates TUTVUMr93, TUTVUMr95, TUTVUMr80 and TUTVUMr85 grouped with B. centrosematis A9 T in Cluster III in nifH phylogeny, but TUTVUMr93 and TUTVUMr95 stood separately and showed proximal relation with sv. glycinearum in nodC phylogeny. Isolates TUTVUMp60 and TUTVUMp75 in Cluster II formed sv. tropici group in both nifH and nodC phylogenies. The isolates in Clusters I and III stood alone, and did not show close relationship with any reference type strains (Figs. 5 and 6).

Influence of the soil environment on bradyrhizobial distribution
CCA analysis was used to correlate bradyrhizobial distribution with soil pH and soil mineral nutrients (namely, Ca, Mg, Na, N, P and K) (Fig. 7). A biplot ordination graph was constructed with only those variables which showed significant influence on bradyrhizobial distribution (p ≤ 0.005 at permutation 999). In the CCA plot, the total mean square contingency coefficient (inertia) was 7.29, of which 0.42 was explainable and 6.8 was unexplainable. Axes CCA1 and CCA2 were able to explain the combined 0.41 of the total explainable variation. The results of the CCA plot indicated that soil pH and Na concentration showed stronger affinity with the first canonical axis (CCA1), although  the effect of Na was in the opposite direction (Fig. 7). According to the length of the arrows and angles with CCA 1 axis, Na showed stronger correlation with the axis. Soil nutrients such as K, Mg, N and P showed significant correlation (P ≤ 0.001) with both canonical axes (CCA1 and CCA2), as well as intermediate influence. However, soil K and N contents were orthogonal to each other (Fig. 7). From the CCA axis, the CCA biplot seemed to suggest that soil pH, as well as Na and Mg levels strongly influenced the genomic variation of isolates collected from Muriaze and Marapyane.

DISCUSSION
Cowpea is a very important legume in southern Africa and is known to be promiscuous in its nodulation. Increased cowpea production for global food security depends on increasing and maintaining higher grain yield, which is currently not possible without major nutrient inputs, especially N and P. Cowpeanodulating rhizobia can be preferred by farmers include those producing growth-promoting hormones such as IAA and those solubilizing phosphate complexes to increase rhizosphere P availability.
This study focused on the analysis of genetic diversity, plant growth-promoting traits and the phylogeny of rhizobia needed to define the criteria for bacterial use in inoculant production for increased cowpea production. Inorganic phosphate solubilisation and IAA production were thus explored in cowpea isolates obtained from South Africa and Mozambique. Although there are a few reports on IAA-producing microsymbionts responsible for Vigna radiata nodulation (Anjum et al. 2011;Ahmad et al. 2013), there is little information on IAA-producing microsymbionts that nodulate Vigna unguiculata L. Walp. Yet Hunter (Hunter 1989) has suggested that IAA in nodules might be involved in several stages of the legume symbiosis. The genes induced by IAA are probably responsible for the many cellular functions of plants (Abel et al. 1996). It has also been suggested that IAA released during rhizobial infection triggers the division of cortical cells in the roots, leading to nodule formation (Datta and Basu 2000).  In this study, most of the IAA-producing rhizobial isolates also exhibited phosphate-solubilizing activity. However, the relationship between the two traits remain unknown. But Chaiharn et al. (Chaiharn and Lumyong 2011) found that the isolates that had the ability to solubilize inorganic phosphate could also produce IAA. It was therefore not surprising that the two traits were highly correlated in this study (see Fig. S3, Supporting Information). The phosphate solubilization results showed a strong variation in the radial distance of the clear zone surrounding the colony periphery, probably indicating that isolate diversity closely mirrored isolate ability to solubilize phosphate.
In this study, the BOX-PCR data revealed the presence of high genetically-diverse rhizobial populations in soils from the agro-ecological regions of South Africa and Mozambique. The 54 authenticated cowpea nodule bacteria used in this study grouped into eight distinct clusters in the Box-PCR analysis. Of the 54 isolates, 45 unique BOX-types were obtained, which suggests that the isolates had variable distances in repetitive and conserved sequences, a finding consistent with the report by Chidebe et al. (Chidebe, Jaiswal and Dakora 2018). The isolates' distribution in clusters within the dendrogram showed that the South African isolates originating from Marapyane exhibited greater diversity when compared to the Mozambican isolates, which showed similar distribution to the results reported by Pule-Meulenberg et al. (Pule-meulenberg et al. 2010). Many isolates were grouped closely together in clusters according to their location of origin in this study, a finding similar to the reports from Senegal, Greece, South Africa and Mozambique (Kouyaté et al. 2014;Tampakaki et al. 2017;Chidebe, Jaiswal and Dakora 2018;Mohammed, Jaiswal and Dakora 2018). The genome of the Biofix (TUTVUI122) commercial inoculant strain did not match with any of the isolates tested, indicating poor nodule occupancy possibly from weak competitiveness.
Environmental factors, especially soil pH and mineral nutrients were found to influence rhizobial community structure and function (Mohammed, Jaiswal and Dakora 2018). The genomic DNA analysis of bradyrhizobial isolates using BOX-PCR and the CCA analysis seem to suggest the existence of a relationship between bradyrhizobial distribution and the physicochemical properties of soil. For example, bradyrhizobial isolates collected from Marapyane, South Africa, were highly influenced by Na content of soil, and by soil pH of isolates from Muriaze in Mozambique. Ndungu et al. (Ndungu et al. 2018) reported that the occurrence and abundance of diverse cowpea rhizobial populations in Kenya was influenced by soil pH, which often affected the bioavailability of mineral nutrients in soils. Recently, Mohammed et al. (Mohammed, Jaiswal and Dakora 2018) and Puozaa et al. (Puozaa, Jaiswal and Dakora 2019) also reported that soil pH and mineral nutrients can influence the diversity of bradyrhizobial populations in South African and Ghanaian soils, and found that N, P and Na concentrations could specifically influence South African isolates. The results of this study also suggested that bradyrhizobial distribution in soil was strongly controlled by soil physico-chemical properties such as pH, and mineral nutrients. Therefore, in the African agroecosystem, edaphic effects can alter the bardyrhizobial distribution in soils (Mohammed, Jaiswal and Dakora 2019).
The evolutionary relationship of the bacterial isolates from cowpea, based on 16S rRNA, housekeeping (atpD, glnII, gyrB and recA) and symbiotic nifH and nodC gene sequences, revealed the presence of diverse bradyrhizobial species in both South African and Mozambican soils. The isolates were grouped in diverse monophyletic clusters within individual gene trees constructed from 16S rRNA, housekeeping and symbiotic gene sequences.
Except for recA phylogeny, all other single gene phylogenies were very congruent to each other. The congruency of isolates from soils of South Africa and Mozambique in all five clusters of the phylogenetic trees could suggest the acquisition of symbiotic genes through vertical gene transfer, except for isolate TUTVUMm110 in atpD phylogeny. In the recA phylogeny, isolates TUTVUMr93, TUTVUMr85, TUTVUMm105 and TUTVUMr94 showed inconsistencies which could be attributed to horizontal gene transfer, subsequent recombination events, and/or differences in the evolutionary history of the gene (Menna and Hungria 2011;Naamala, Jaiswal and Dakora 2016;Jaiswal, Msimbira and Dakora 2017;Chidebe, Jaiswal and Dakora 2018). An important consideration is that, in general, all groups in the phylogenies of the housekeeping genes atpD, glnII and gyrB, as well as the symbiotic nifH and nodC genes were similar to the 16S rRNA gene tree, thus confirming the usefulness of these genes for phylogenetic reconstruction of the genus Bradyrhizobium. In this study, geographic origin did affect the phylogenetic placement of isolates in the housekeeping gene trees, which diminished the hypothesis of a common origin for Bradyrhizobium and supported the evolution of these isolates by allopatric speciation. For example, isolates from Marapyane, Muriaze and Ruace did not cluster together, and/or aligned with any known Bradyrhizobium reference type strains.
To get a clearer resolution in the relationship between test isolates and type reference strains, phylogenetic information was deduced from combined sequence analysis of more than one gene (Chidebe, Jaiswal and Dakora 2018). In this study, concatenated atpD+glnII+gyrB gene sequence analysis yielded the same results as the individual phylogenetic trees, with the isolates grouping in five clusters. In Cluster I, isolate TUTVUMp6 obtained from South Africa, and isolates TUTVUMr103, TUTVUMm112, TUTVUMm113 and TUTVUMm114 from Mozambique, all showed very close relationship with B. yuanmingense CCBAU 10071 T which was originally isolated from the genus Lespedeza in China (Yao et al. 2002). Isolates TUTVUMp60, TUTVUMp75 from South Africa and TUTVUMm117 from Mozambique shown in Cluster II, was closely related to B. pachyrhizi PAC48 T originally obtained from root nodules of Pachyrhizus erosus in Costa Rica (Ramírez-Bahena et al. 2009). Recently, B. pachyrhizi was also identified as a cowpeanodulating microsymbiont in Angola (Grönemeyer et al. 2014), Spain (Bejarano et al. 2014), South Africa (Mohammed, Jaiswal and Dakora 2018) and Mozambique (Chidebe, Jaiswal and Dakora 2018). Isolate TUTVUMm110 in Cluster IV obtained from Mozambican soil and the commercial Biofix strain TUTVUI122 were very close to B. arachidis CCBAU 055107 T , a finding consistent with the results of Chidebe et al. (Chidebe, Jaiswal and Dakora 2018). Due to the absence of gyrB sequence of B. arachidis CCBAU 055107 T , it was excluded in the concatenated phylogeny. The isolates in Clusters IIa, III and V did not show any close relationship with known reference type strains of Bradyrhizobium, and could therefore be novel Bradyrhizobium species in South African and Mozambican soils, this indicating allopatric speciation.
Taken together, the bradyrhizobial isolates from this study exhibited varying morpho-physiological and molecular characteristics, which clearly suggested that uncultivated soils of the South African site and farmers' fields in Mozambique harboured diverse and novel bradyrhizobia of African origin capable of releasing IAA and solubilising P for plant growth. The monophyletic grouping of the isolates within phylogenetic trees could suggest their geographic isolation in the African continent, and hence prove their allopatric speciation. The acidic soils and other unique edaphic factors were probably responsible for the presence of these diverse indigenous and novel bradyrhizobia nodulating cowpea in Africa. These isolates should be further studied systematically to ascertain their identity and symbiotic performance for subsequent testing under field conditions to identify strains for inoculant production.

SUPPLEMENTARY DATA
Supplementary data are available at FEMSEC online.