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

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. © FEMS 2019.

INTRODUCTION

Legumes are a major component of agricultural systems throughout the world due to their N2-fixinpan>g ability when inpan> symbiosis with soil bacteria called rhizobia. This makes nodulated legumes good candidates for improvinpan>g 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). N2 fixation is usually performed by the bacteria in root nodules, and the NH3 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 inpan>digenous legume crops inpan> Africa, with both wild and cultivated forms growinpan>g across the continpan>ent (Mohammed, Jaiswal and Dakora 2018). West Africa is the centre of originpan> of this crop (Jaiswal and Dakora 2019; Puozaa, Jaiswal and Dakora 2019). Cowpea can meet 66% to 96% of its N nutrition from symbiotic fixation (Jaiswal and Dakora 2019).

The high cost of N fertilizers and the environmental pollution associated with their frequent use to inpan>crease grainpan> 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 N2, solubilize soil P, and produce plant growth-promoting 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 inpan>volved inpan> root formation and root elongation for inpan>creased uptake of water and nutrients (Etesami et al. 2009). Thus, the identification of IAA-producing and P-solubilizing traits among N2-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 N2 fixation (Chidebe, Jaiswal and Dakora 2018). Such studies have inpan>creased our unpan>derstandinpan>g of the genetic diversity of rhizobia and signpan>ificantly improved rhizobial classification through the re-definpan>ition of species. Previous reports from the different parts of the world such as Botswana, Ghana, South Africa, Senegal, Mozambique, Ethiopia, Greece, Spainpan> 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 N2-fixing rhizobia in Southern Africa. Since cowpea is an indigenous legume crop of the vast and agro-ecologically 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.

MATERIALS AND METHODS

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, 26o 43ʹ 44.8ʹʹ E and 15o 09ʹ 12.9ʹʹ S, 30o 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 m2. 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 intra-row spacing of 15 cm.

Figure 1.

Test locations in South Africa and Mozambique, and rhizobial distribution in different clusters based on BOX-PCR analysis.

Table 1:

Colony types, biochemical and molecular characterization of cowpea-nodulating isolates.

			Colony characteristic					Biochemical characteristic		Box-PCR Type
Isolate	Location	Genotype	Appearance day	Size	Texture	Shape	colour	IAA µg ml−1	PSI
South Africa
TUTVUMp5	Marapyane	IT97K-1069–6	8	1.5	Gummy	Irregular	Translucent	26.38±0.05	1.43	1
TUTVUMp6	Marapyane	IT97K-390–2	4	1	Gummy	Conical	Whitish	2.16±0.05	2.72	30
TUTVUMp7	Marapyane	IT97K-390–2	7	3	Gummy	Domed	Milky	15.60±0.05	2	5
TUTVUMp8	Marapyane	IT00K-126–3	6	2	Dry	Irregular	Whitish	19.28±0.5	2.13	33
TUTVUMp11	Marapyane	IT06K-128	5	1	Dry	Domed	Translucent	17.25±3	3.33	5
TUTVUMp12	Marapyane	IT00K-126–3	7	2	Gummy	Irregular	Milky	26.77±0.05	1.47	5
TUTVUMp26	Marapyane	IT00K-126–3	8	1.5	Dry	Conical	Whitish	13.42±1	1.47	5
TUTVUMp28	Marapyane	Sudan 1	7	2	Gummy	Watery	Translucent	28.01±2.5	1.27	5
TUTVUMp37	Marapyane	IT00K-126–3	8	4.5	Gummy	Irregular	Milky	8.89±1.5	1.62	2
TUTVUMp41	Marapyane	IT97K-390–2	6	4	Gummy	Watery	Translucent	31.55±0.005	1.67	5
TUTVUMp48	Marapyane	IT18	9	2	Gummy	Domed	Whitish	35.02±3	0	16
TUTVUMp49	Marapyane	IT97K-390–2	4	1	Gummy	Irregular	Translucent	44.63±3	1.54	22
TUTVUMp50	Marapyane	Sudan 1	7	2.5	Gummy	Irregular	Milky	20.85±2	1.55	5
TUTVUMp53	Marapyane	IT18	8	6	Gummy	Domed	Whitish	39.21±5	2.5	7
TUTVUMp56	Marapyane	Sudan 1	8	2	Gummy	Irregular	Milky	17.89±4.5	1.76	6
TUTVUMp57	Marapyane	IT18	4	6	Gummy	Watery	Translucent	41.54±1.5	2.5	8
TUTVUMp58	Marapyane	Sudan 1	6	3	Gummy	Irregular	Translucent	41.88±0.5	2.16	35
TUTVUMp59	Marapyane	Sudan 1	3	6	Gummy	Irregular	Milky	35.02±3	1.68	3
TUTVUMp60	Marapyane	IT97K-1069–6	5	4	Elastic	Irregular	Milky	8.25±3	0	32
TUTVUMp65	Marapyane	IT97K-1069–6	8	2	Gummy	Domed	Whitish	26.22±1	2	17
TUTVUMp75	Marapyane	IT97K-1069–6	4	2.5	Dry	Irregular	Translucent	35.48±0.5	2.9	25
Mozambique
TUTVUMr80	Ruace	IT00K-126–3	6	3.5	Gummy	Domed	Whitish	1.01±2	3.55	18
TUTVUMr82	Ruace	IT00K-126–3	5	3	Gummy	Irregular	Translucent	37±4	0	20
TUTVUMr83	Ruace	IT97K-1069–6	7	1	Dry	Conical	Whitish	40±2	2	19
TUTVUMr84	Ruace	IT97K-1069–6	6	1	Dry	Conical	Milky	14.32±2	0	29
TUTVUMr85	Ruace	IT97K-1069–6	6	2	Gummy	Domed	Whitish	29.59±5.5	0	23
TUTVUMr89	Ruace	IT18	10	0.5	Dry	Circular	Whitish	2.88±2	2.08	24
TUTVUMr90	Ruace	IT18	3	6	Gummy	Domed	Milky	22.22±3	1.22	11
TUTVUMr93	Ruace	IT97K-1069–6	3	3	Gummy	Domed	Milky	0.64±0.1	1.38	9
TUTVUMr94	Ruace	IT97K-1069–6	6	4.5	Gummy	Irregular	Whitish	4.29±0.05	0	31
TUTVUMr95	Ruace	IT97K-1069–6	3	6	Elastic	Irregular	Milky	6.926±3.5	0	27
TUTVUMr99	Ruace	IT97K-1069–6	6	6	Gummy	Irregular	Milky	5.7±1	0	10
TUTVUMr100	Ruace	IT97K-1069–6	8	1	Dry	Conical	Whitish	4.20±0.05	2.5	26
TUTVUMr101	Ruace	IT18	7	3	Gummy	Flat	Whitish	2.58±0.05	2.33	19
TUTVUMr102	Ruace	IT00K-126–3	3	6	Elastic	Domed	Milky	46.38±1.5	1.6	12
TUTVUMr103	Ruace	IT16	4	6	Gummy	Irregular	Milky	1.94±1.5	0	21
TUTVUMm104	Muriaze	IT97K-1069–6	10	1	Dry	Conical	Whitish	34.98±3	1.81	19
TUTVUMm105	Muriaze	IT97K-1069–6	9	2	Gummy	Irregular	Whitish	40.76±2.5	0	19
TUTVUMm106	Muriaze	IT97K-1069–6	5	2	Elastic	Domed	Whitish	4.49±0.05	0	19
TUTVUMm107	Muriaze	IT97K-1069–6	7	3	Gummy	Domed	Whitish	46.73±0.05	1.83	19
TUTVUMm108	Muriaze	IT97K-1069–6	3	6	Elastic	Irregular	Milky	24.63±0.05	0	6
TUTVUMm109	Muriaze	IT97K-1069–6	6	4	Gummy	Irregular	Milky	21.67±2	2.36	6
TUTVUMm110	Muriaze	IT97K-1069–6	6	2.5	Gummy	Irregular	Milky	40.45±0.05	2.18	22
TUTVUMm111	Muriaze	IT97K-1069–6	7	5	Gummy	Watery	Translucent	35.04±0.05	4	22
TUTVUMm112	Muriaze	IT16	6	3	Elastic	Domed	Whitish	1.07±0.05	1.5	28
TUTVUMm113	Muriaze	IT16	7	2.5	Gummy	Domed	Whitish	46.01±1	0	13
TUTVUMm114	Muriaze	IT16	5	2	Gummy	Domed	Whitish	22.89±1.5	0	14
TUTVUMm115	Muriaze	IT18	6	2	Gummy	Domed	Whitish	42.19±1.5	0	15
TUTVUMm116	Muriaze	IT04K-227–4	6	2.5	Gummy	Domed	Whitish	26.56±2.5	0	19
TUTVUMm117	Muriaze	IT04K-227–4	6	5	Gummy	Irregular	Milky	20.13±1.5	0	4
TUTVUMm118	Muriaze	IT97K-1069–6	5	3	Gummy	Domed	Milky	50.12±0.01	1.57	19
TUTVUMm119	Muriaze	IT16	4	6	Gummy	Domed	Milky	56.46±1.5	1.71	19
TUTVUMm120	Muriaze	IT97K-1069–6	5	4.5	Gummy	Irregular	Milky	52.56±4	1.47	19

Colony types, biochemical and molecular characterization of n class="Species">cowpea-nodulatinpan>g isolates.

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 inpan> runpan>ninpan>g 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 cm3) 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 (∼ 106–107 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.

Biochemical characterization of rhizobia

Two biochemical properties were assessed, namely, the production of the root growth-promoting hormone auxin (indole-3-acetic acid) and P-solubilization. Production of IAA by rhizobial cells was estimated, as described by Raddadi et al (Raddadi et al. 2008). Standard curves were generated from absorbances of 0.5,10,10,20,50 and 100 µg mL−1 IAA measured at 530 nm (Fig. S1, Supporting Information).

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 inpan>cubated at 28°C for 7 days. The formation of a clear halo zone arounpan>d a bacterial spot was taken as inpan>dicative of P solubilisation. The diameter of each halo zone was measured to identify isolates havinpan>g high P-solubilisinpan>g 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 usinpan>g Sigma's Bacterial Genomic DNA Kit according to the manufacturer's instructions (GenEluteTM), 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 DocTM 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, proteinpan>-codinpan>g (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) accordinpan>g to the manufacturer's inpan>struction. 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 cowpea-nodulatinpan>g bradyrhizobia was examinpan>ed usinpan>g canonical correspondence analysis (CCA) with vegan (version 2.4–2) (Oksanen et al. 2010) of R software (R Core Team 2015). Anpan>alysis was done with only the soil parameters which showed signpan>ificant ecological contribution to the distribution of the test cowpea-nodulating bradyrhizobia. The general permutation test was used to assess the statistical significance of the ordination axes

RESULTS

Plant nodulation data

A total of 102 isolates were obtained from the cowpea nodules collected from all three study sites inpan> South Africa and Mozambique. Out of the102 isolates, 61 could elicit nodule formation inpan> 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.

Colony characteristics and biochemical behavior of the rhizobial isolates

Colony texture showed that 72% of the isolates were gummy, 17% dry, and 11% elastic. Colony sizes also differed, with 55% having a diameter of 2 to 4 mm, 28% a diameter of 4.5 to 6 mm, while 17% showed ≤1.5 mm in diameter. About 64% of the rhizobial isolates appeared in 6 to 10 days, while 36% took 3 to 5 days to appear on YMA plates. In terms of colony shape, 42% were irregular, 36% domed, 11% conical, 7.5% watery, 1.8% circular and 1.8% flat. Colony colour also differed, with 41% being whitish, 41% milky, and 18% translucent on YMA plates (Table 1).

IAA production and phosphate-solubilising activities by authenticated isolates

The test rhizobial isolates varied in their ability to produce IAA, with concentrations ranginpan>g 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 3PO4). 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 (inpan>cludinpan>g the commercial isolate TUTVUI122) resulted inpan> distinpan>ct bandinpan>g 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 BLASTn 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–3T and B. subterraneum 58 2–1T were the closest related species with 99.1%–99.5% sequence identity. Cluster II comprised two sub-clusters (IIa and IIb). Isolates TUTVUMp75, TUTVUMr82 and TUTVUMr99 which shared 99.7% to 100% sequence identity with B. elkanii USDA 76T and B. pachyrhizi PAC48T in sub-cluster IIa, and isolate TUTVUMm117 grouped with B. ferriligni in sub-cluster IIb. Isolates TUTVUMr80, TUTVUMr85, TUTVUMr93 and TUTVUMr95 shared 99.7%–100% identical sequences with the reference type strains of B. ganzhouense RITF806T, B. gangdongens CCBAU 51649T and B. manausense BR 3351T in Cluster III. Cluster IV comprised isolates TUTVUMm108 and TUTVUI122 which had a 99.8% sequence similarity with reference type strains B. centrosematis A9T, B. guangxiense CCBAU 23303T, B. huanghuaihaiense CCBAU 23303T, B. ingae BR 10250T, B. stylosanthis BR 446T and B. iriomotense NBRC 102520T. 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–3T, and B. subterraneum 58 2–1T with 75% bootstrap support in Cluster V while isolate TUTVUMp5 showed outgroup of Cluster V (Fig. 2).

Figure 2.

Maximum-likelihood molecular phylogenetic analysis of 16S rRNA nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together in the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site. The analysis involved 76 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

Maximum-likelihood molecular phylogenetic analysis of 16S rRNA nucleotide sequences of rhizobia nodulatinpan>g cowpea collected from South Africa and Mozambique. The associated taxa clustered together in the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site. The analysis involved 76 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

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 05623T and B. centrosematis A9T 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).

Figure 3.

Maximum-likelihood molecular phylogenetic analysis of recA nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together in the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site. The analysis involved 82 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

Maximum-likelihood molecular phylogenetic analysis of recA nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together inpan> the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed usinpan>g the Kimura 2-parameter method and are inpan> the unpan>its of the number of base substitutions per site. The analysis inpan>volved 82 nucleotide sequences. Codon positions inpan>cluded were 1st+2npan>d+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

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 linpan>eage. A maximum likelihood phylogeny was constructed usinpan>g 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 10071T 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 PAC48T. 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–3T and B. arachidis CCBAU 45332T as the closest related species (Fig. 4).

Figure 4.

Maximum-likelihood molecular phylogenetic analysis of concatenated (atpD+glnII+gyrB) nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together in the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site. The analysis involved 67 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7. Clusters lines in the map showed isolates belong to that locations.

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

Figure 5.

Maximum-likelihood molecular phylogenetic analysis of nifH nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together in the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site. The analysis involved 73 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

Figure 6.

Maximum-likelihood molecular phylogenetic analysis of nodC nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together in the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site. The analysis involved 45 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

Maximum-likelihood molecular phylogenetic analysis of nifH nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together inpan> the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed usinpan>g the Kimura 2-parameter method and are inpan> the unpan>its of the number of base substitutions per site. The analysis inpan>volved 73 nucleotide sequences. Codon positions inpan>cluded were 1st+2npan>d+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

Maximum-likelihood molecular phylogenetic analysis of nodC nucleotide sequences of rhizobia nodulating cowpea collected from South Africa and Mozambique. The associated taxa clustered together inpan> the 1000 bootstrap support and percentage are shown next to the branches. The evolutionary distances were computed usinpan>g the Kimura 2-parameter method and are inpan> the unpan>its of the number of base substitutions per site. The analysis inpan>volved 45 nucleotide sequences. Codon positions inpan>cluded were 1st+2npan>d+3rd+Noncoding. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7.

The isolates in Cluster IV were closely related to B. arachidis CCBAU051107T 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–3T with 98.5% sequence identity in the nifH phylogeny. Isolates in Cluster V showed a close relationship with B. subterraneum 58–2-1T 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 A9T 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.

Figure 7.

CCA plot representation of association between soil environmental variables and BOX-PCR banding pattern data of rhizobial isolates.

DISCUSSION

Cowpea is a very important legume inpan> southern Africa and is knpan>own to be promiscuous inpan> 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. Cowpea-nodulating 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-producinpan>g 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 inpan> this study grouped inpan>to eight distinpan>ct clusters inpan> the Box-PCR analysis. Of the 54 isolates, 45 unpan>ique BOX-types were obtainpan>ed, which suggests that the isolates had variable distances inpan> repetitive and conserved sequences, a finpan>dinpan>g consistent with the report by Chidebe et al. (Chidebe, Jaiswal and Dakora 2018). The isolates’ distribution inpan> clusters withinpan> the dendrogram showed that the South African isolates originpan>atinpan>g 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 usinpan>g 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 agro-ecosystem, 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 10071T 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 PAC48T 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 cowpea-nodulating 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 055107T, 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 055107T, 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 solubilisinpan>g P for plant growth. The monophyletic groupinpan>g of the isolates withinpan> phylogenetic trees could suggest their geographic isolation inpan> the African continpan>ent, and hence prove their allopatric speciation. The acidic soils and other unpan>ique edaphic factors were probably responsible for the presence of these diverse inpan>digenous and novel bradyrhizobia nodulatinpan>g 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.

FUNDING

This work was supported with grants from the Bill and Melinda Gates Foundation Project on Capacity Building in Legume Sciences in Africa, the South African Department of Science and Technology, the Tshwane University of Technology, the n class="Chemical">National Research Foundation inpan> Pretoria, and the South African Research Chair inpan> Agrochemurgy and Plant Symbioses.

Click here for additional data file.