The Rhizobial Microbiome from the Tropical Savannah Zones in Northern Côte d’Ivoire

Over the past decade, many projects have been initiated worldwide to decipher the composition and function of the soil microbiome, including the African Soil Microbiome (AfSM) project that aims at providing new insights into the presence and distribution of key groups of soil bacteria from across the African continent. In this national study, carried out under the auspices of the AfSM project, we assessed the taxonomy, diversity and distribution of rhizobial genera in soils from the tropical savannah zones in Northern Côte d’Ivoire. Genomic DNA extracted from seven sampled soils was analyzed by sequencing the V4-V5 variable region of the 16S rDNA using Illumina’s MiSeq platform. Subsequent bioinformatic and phylogenetic analyses showed that these soils harbored 12 out of 18 genera of Proteobacteria harboring rhizobia species validly published to date and revealed for the first time that the Bradyrhizobium genus dominates in tropical savannah soils, together with Microvirga and Paraburkholderia. In silico comparisons of different 16S rRNA gene variable regions suggested that the V5-V7 region could be suitable for differentiating rhizobia at the genus level, possibly replacing the use of the V4-V5 region. These data could serve as indicators for future rhizobial microbiome explorations and for land-use decision-making.


Introduction
Since the advent of sequencing technologies, the determination of microbial diversity has become a major topic of interest [1]. Over the last decade, for example, many small-or broad-scales initiatives have been launched around the world to decipher the composition and function of soil microbiome [2][3][4][5][6][7], including the African Soil Microbiome (AfSM) Project implemented in Sub-Saharan Africa (SSA) [8,9]. This unique multi-national project, implemented in a dozen SSA countries, is the first such study to ever be undertaken in Africa at this scale [9,10]. It was launched in 2016 to provide new insights into the presence and the distribution of key groups of soil bacteria, including the rhizobia, by using the high-throughput amplicon sequencing (HTAS) and phylogeny of the 16S rRNA gene [9,10].
Rhizobia are Gram-negative saprophytic Alphaand Beta-proteobacteria that play a key role in nitrogen biochemical cycling [11,12]. They form a polyphyletic group of bacteria among the lineages of prokaryotes capable of reducing atmospheric dinitrogen (N 2 ) Sudan savannah in the north, both of which constitute more than two-thirds of the entire savannah region, and the Guinean savannah. The Guinean savannah, which is located in the southern part of the savannah biome, is sometimes referred to as the transition zone, even though the entire savannah region is transitional between the narrow belt of forest paralleling the coastline and the Sahara [47]. While the Guinean savannah zone has been well studied, including the phylogenetic relationships, ecological niches and functional roles of N-cycling bacteria [48,49], the biodiversity and ecology of microorganisms from the two other savannah zones remain relatively unexplored.
The aim of this study was to assess the taxonomy, diversity and distribution of rhizobial genera in soils from the Sudan and the Sub-Sudan tropical savannah zones in Northern Côte d'Ivoire using high-throughput sequencing of the 16S rDNA variable V4-V5 region. In addition, the different 16S rRNA gene variable regions were compared in silico to assess their effectiveness for differentiating all the genera of rhizobia validly published to date.

The Study Area
The study was carried out in the context of the African Soil Microbiome (AfSM) project, and the studied area was located in the savannah zones in Northern Côte d'Ivoire. Côte d'Ivoire is divided roughly into two large agro-ecological regions, of which the northern savannah region, where food crops, cotton and livestock predominate, and the fertile forest zone of the south, where most of the country's cash crops are produced [50]. The boundary that marks the transition from forest to savannah is remarkably irregular (Figure 1a). It is characterized by the presence of an inverted triangular-like structure known as the « V-baoulé » (see [51,52]). Historically, the country was also divided into five zones according to the vegetation types, including (from the far north to the south) the Sudan savannah (I), the Sub-Sudan savannah (II), the Guinean savannah (III), the Semi-deciduous moist forest (IV) and the Evergreen moist forest (V) [46,53,54] (Figure 1b).
The study area covered the Sudanian savannah (I) and the Sub-Sudanian savannah (II) zones ( Figure 1c). The annual rainfall is among the lowest in the country [53,55], ranging approximately from 1000 to 1750 mm per year [52,53,55]. The Sub-Sudanian savannah and the Sudanian savannah zones are also characterized by an average annual humidity of 60-70%, annual mean temperature of 24-27 • C and ferralitic and ferruginous soils [55,56]. The vegetation consists of grasslands, wooded grasslands and gallery forests [46,56]. Narrow gallery forests extend along watercourses and drainage basins, where very tall trees, such as Ceiba pentandra, Sterculia tragacantha and Triplochiton scleroxylon, dominate. The dominant tree species of the wooded grasslands are Acacia albida, Khaya senegalensis, Parkia biglobosa, and Tamarindus indica, and herbaceous plants include Andropogon tectorum and Pennisetum purpureum. Dominant trees in savannahs consist also of Butyrospermum parkii, Daniellia oliveri and Lophira lanceolata, as well as Andropogon ivorensis, Loudetia simplex and Panicum phragmitoides for herbaceous species [53,56]. As for the cultivated fields, they consist mainly of cashew trees and cereals (maize and rice). The sampled soils belong to the rhizosphere of all the vegetation types we described (Table 1).

Soil Sampling
Soil samples were collected in August-September 2017 from seven sites located in five administrative regions (Table 1) and alongside the national roads. Each sampled soil belongs to the rhizosphere of a natural herbaceous or wooded vegetation and/or a cultivated plant species (cashew, maize, rice etc.) ( Table 1). The distance between sampling sites spanned 50-300 km. Each sampling site was represented by an area of approximately 100 m × 50 m with four independent sample locations (a virtual 1 m 2 quadrat) at the corners of the oblong ( Figure S1). At each of the four independent sample locations, four topsoil cores (2 cm in diameter and 5 cm in depth) (pseudo-replicate samples) were collected, pooled together, and homogenized into a composite sample of approximately 25 g (replicate sample). Four independent replicate samples (4 × 25 g) obtained from four sample locations at each sampling site were kept in a labelled sterile plastic bag and formed an independent soil sample. This process was repeated for all seven sampling sites. In total, seven independent soil samples were obtained. Each soil sample taken from the two savannah zones from northern Côte d'Ivoire (CI) is referred to by the soil numbers 11, 13, 14, 17, 18, 20 or 44 (Table 1). After sampling, the soil samples were transported to the laboratory, where they were stored at 4 • C prior to a shipment to South Africa for further analysis.

Analysis of Soil Physicochemical Properties
The analysis of soil physico-chemical characteristics was carried out by Bemlab (Strand, Cape Province, South Africa) using standard methods. Briefly, prior to analyses, the samples were air-dried at room temperature for four days, separated from roots and debris, and passed through a 2 mm sieve. The sieved replicate samples of each sampling site were subsequently pooled together to obtain a composite soil sample. Physical characteristics (fractions of clay, sand and silt) were analyzed using the Bouyoucos sedimentation method (hydrometer method) [57]. The classification of soils according to texture was based on the standard USDA particle-size classification using the Soil Texture online Calculator (https: //www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/survey/?cid=nrcs142p2_054167; accessed on 8 April 2021). The pH (aqueous) was measured as described [58], while the oxidizable carbon content was determined using the dichromate oxidation method (the Walkley-Black method) [59]. Soil chemical parameters (exchangeable and soluble Na, K, Ca, Mg, Al, Fe, Mn and P) were analyzed using the Mehlich No. 3 soil test extractant with Inductively Coupled Atomic Emission Spectrometry (ICP-AES) procedures [60].

DNA Extraction, PCR, MiSeq Sequencing, and Sequence Data Analysis
Genomic DNA (gDNA) extraction, amplification and high-throughput amplicon sequencing were carried out as in Nkuekam et al. [61], with few modifications. DNA extraction was conducted at the Centre for Microbial Ecology and Genomics (University of Pretoria, South Africa). Briefly, soil samples were first ground with Powerlyser (Mo Bio Laboratories Inc.), and the genomic DNA was extracted from 0.25 g of soil using the PowerSoil DNA isolation kit (Mo Bio Laboratories Inc., Carlsbad, CA, USA). The success of the extraction was verified by 1% agarose gel electrophoresis visualizing under UV light. DNA amplification was conducted at the MRDNA sequencing facility (www. mrdnalab.com, accessed on 25 June 2021, Shallowater, TX, USA) in a 30-cycle PCR using the HotStarTaq Plus Master Mix Kit (Qiagen, Germantown, MD, USA) [61]. The V4-V5 variable region of the 16S rRNA gene was amplified and sequenced using the alternative forward primer 515F-Y (5 -GTGYCAGCMGCCGCGGTAA-3 ; [62]) and the universal reverse primer 909-928 (5 -CCCCGYCAATTCMTTTRAGT-3 ; [63]), with 12 nucleotides unique barcode at 5-end of 515F-Y for each soil sample. High-throughput amplicon sequencing was performed using an Illumina MiSeq platform at the MRDNA sequencing facility.
For the processing of the sequencing data, raw sequences were first checked for reads quality using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/; accessed on 1 May 2021). Reads were then sorted based on unique soil sample tags using Sabre 1.0 program (https://github.com/najoshi/sabre; accessed on 1 May 2021) with default parameters and trimmed for primer and barcode removal using cutadapt 2.10 [64]. The trimmed sequences were subsequently denoised using the DADA2 algorithm [65] that resolves Illumina-sequenced amplicon errors to generate amplicon sequence variants (ASVs). ASVs were classified using the RDP classifier [66] with default parameters. The assignment of ASVs to rhizobia taxa was processed through a search for similar sequences conducted with BLAST v. 2.9.0+ [67] against the SILVA 138.1 database [68].
The taxonomic affiliations obtained with SILVA was manually validated at the genus level using several approaches, including (i) BlastN with NCBI/GenBank online standard databases (nucleotide collection (nr/nt) and Whole-genome shotgun contigs (wgs)) to select closely related reference sequences, (ii) phylogenetic reconstructions and (iii) similarity level calculations between the ASVs and selected rhizobial reference sequences. Briefly, V4-V5 edited sequences were aligned with MUSCLE as implemented in MEGA software v. 7 and phylogenies were inferred subsequently with evolutionary trees reconstructed using the maximum likelihood (ML) and neighbor-joining (NJ) methods [69,70]. Bestfit nucleotide substitution models were selected according to the Bayesian information criterion (BIC) [71] and the uncorrected genetic distances calculated as in Rashid et al. [72]. Phylogenetic analyses also included 16S rRNA gene sequences derived from archived genome data of the type species of all the 18 alphaproteobacterial and betaproteobacterial genera harboring rhizobia isolates [15], except for Paraburkholderia. P. graminis PHS1 16S rRNA gene data was used in the analysis as a surrogate of P. graminis type strain C4D1M, for which no full rRNA gene sequence was accessible at the time of writing (June 2021). Details of the type species of all the 18 bacterial genera harboring rhizobia are listed in Table 2.

Statistical and Diversity Index Analyses
The statistical and diversity index analyses were performed using R v. 4.0.3 [95], including the R packages vegan [96], phyloseq [97] and ggplot2 [98]. The rarefaction curves were computed using the vegan function rarefy, which is based on Hurlbert's formula [99] to evaluate the sequencing efforts provided. As a normalization step to reduce bias associated with different sequencing depths, all samples were subsampled down to the size of the smallest sample. Each sample was rarefied to 1384 reads. Indices of richness (Chao1) and alpha diversity (Shannon, Simpson and Fisher) were calculated by savannah zone, and a non-parametric Wilcoxon test was used to compare the mean values at the significance level of 5%. The degree of community differentiation (beta-diversity) was evaluated to calculate Jaccard's similarity coefficient and the Bray-Curtis index of (dis)similarity for each ASV. The relationship between the ASVs of rhizobia and the environmental parameters that characterize the soils of the savannah zones of Northern Côte d'Ivoire was assessed by canonical correspondence analysis (CCA) using ten physico-chemical parameters. Prior to drawing the relative abundance of rhizobia taxon per soil sample, the sample counts of ASVs were used to calculate relative abundance by computing the ratio of the count of each sample by the sum of the counts of all samples. The obtained relative abundance of counts was used to draw the bar plot of relative abundance of ASVs by genera and family between samples. The heatmap was created using the ecologically organized plot_heatmap function of the phyloseq package, which is a variant of the heatmap provided by the NeatMap package [100]. To draw the heatmap, we used the NMDS ordination method and the Bray-Curtis distance. The 16S rRNA gene sequences used in this study are available in the NCBI SRA database under accession number SRR13623326 (CI11), SRR13623324 (CI13), SRR13623323 (CI14), SRR13623320 (CI17), SRR13623319 (CI18), SRR13623317 (CI20) and SRR13623335 (CI44).

In Silico Evaluation of the 16S rRNA Gene V-Regions Discriminatory Power for Rhizobia
The aim of this analysis was to compare in silico the discriminatory power at the genus level of nine commonly used 16S rRNA gene V-regions for rhizobia. Prior to this analysis, we estimated the current number of species of rhizobia. We counted the number of species of rhizobia validly published within the 18 alphaproteobacterial and betaproteobacterial genera harboring rhizobial species and provided on the List of Prokaryotic names with Standing in Nomenclature (LPSN), also known as bacterio.net [101]. We also documented the prevalence of species with nodulation (Nod) and/or N 2 -fixation (Fix) capacities using the original publications describing novel taxa of rhizobia and accessible on the LPSN website, as well as the more recent publications that reviewed the symbiotic features of rhizobial taxa [102][103][104]. All the data are reported in Table 2. As for the evaluation of the discriminatory power of V-regions, nine V-regions (V1 to V9) spanning the entire 16S rRNA gene and commonly targeted in microbial metagenomic analyses were selected [63,[105][106][107] ( Table 3). The corresponding universal primers targeting the selected V-regions and their relevant characteristics are reported in Table 3. The V-regions were compared to the fulllength size of the 16S rRNA gene sequences with a method used in VanInsberghe et al. [108]. Briefly, the full-length 16S rRNA gene sequences were aligned for all the 18 genera with MAFFT version 7 using the Q-INS-I method [109] and followed by a maximum likelihood phylogenetic tree reconstruction as well as by a pairwise similarity distances calculation as in Rashid et al. [72]. Similarity values were used to identify the uniquely distinguishable taxa at 97%, 99% or 100% cut-offs. Subsequently, the full-length 16S rRNA gene sequence alignment was edited to the total number of positions that corresponds to those of each V-region, in addition to that of the V1-V9 region, which is the near-full-length size of the 16S rRNA gene ( Table 3). The total number of positions in each edited dataset was used to calculate the similarity values that served to identify the uniquely distinguishable taxa for the given V-region.

Results
We carried out a high-throughput amplicon sequencing (HTAS) analysis of the 16S rRNA gene V4-V5 region to assess the taxonomy, diversity and distribution of rhizobial taxa in seven soils in the Sudan savannah (I) and the Sub-Sudan savannah (II) zones in Northern Côte d'Ivoire ( Figure 1). These two zones have been largely neglected in terms of fundamental research in microbial ecology, and this study provides their first comprehensive rhizobial microbiome analysis. The sampled soils from the two zones were analyzed for their physico-chemical properties prior to the HTAS analysis.

Physico-Chemical Properties of Soil Samples
The soils' physico-chemical data are reported in Supplementary Table S1. The seven studied localities have similar soil textures characterized by a high proportion of sand (>70%) but can be divided into two subgroups: soils CI13, CI14, CI17, CI18 and CI20 were sandy loams, while CI11 and CI44 were loamy sand. The pH of the seven soils ranged from 5 to 7, being consistent with the soil pH range expected in tropical humid regions (https://www.qld.gov.au/environment/land/management/soil/soil-properties/ ph-levels; accessed on 8 April 2021). Soils CI11 and CI13, both of which were from the locality of Bouna in the north-east (Figure 1), were neutral (pH 6.6), while CI14, CI17, CI18, CI20 and CI44 soils were acidic (pH < 6.5). CI14 (pH = 6.4) was the least acidic soil (nearly neutral). The distribution of the soil samples according to the chemical properties was more heterogeneous. CI11 was among the soil samples having the highest values of calcium (Ca 2+ ), magnesium (Mg 2+ ), potassium (K + ) and sodium (Na + ), while CI44 and CI-17 had the lowest values for the same mineral elements.

Sequence Data and Taxonomic Affiliation
The amplification of the total DNA extracted from the seven soil samples using V4-V5 primers yielded ca. 400-500 nucleotide length products, as expected for bacteria ( Table 3). The rarefaction curves reach the asymptote with less than 1000 sequences, suggesting that the sequencing effort of each amplicon was sufficient ( Figure S2). From a total of 900,760 sequences obtained through Illumina's high throughput sequencing platform, a total of 786,283 sequences were considered for the clustering after sequence trimming. When clustered and quality-controlled, the 786,283 sequences yielded 5997 ASVs in total, of which 80 (less than 2%) matched to rhizobia in the SILVA database. This assignment of the ASVs to rhizobia taxa was further refined using a multi-step approach that includes phylogenetic analyses of the ASVs (Figures 2 and S3) as well as a genetic distance comparison (Table S2) and online blastN analyses. Phylogenetic assignments of the ASVs performed with a subset of 86 closely related sequences (99 to 100% similar) and/or 18 type species of Proteobacteria harboring rhizobia species validly published to date yielded similar taxonomic affiliations (see Figures 2 and S3, respectively). Together, these different analyses improved the taxonomic identification, with 77 ASVs (equivalent to 15,886 sequences) being confirmed as rhizobia (Figure 2; Tables S3 and S4). The 77 ASVs belonged to 12 genera of rhizobia (Bradyrhizobium, Cupriavidus, Devosia, Ensifer, Mesorhizobium, Methylobacterium, Microvirga, Neorhizobium, Paraburkholderia, Rhizobium, Shinella and Trinickia) ( Figure 2) and were present in six families (Burkholderiaceae, Devosiaceae, Methylobacteriaceae, Nitrobacteraceae, Phyllobacteriaceae and Rhizobiaceae) of the classes Alphaproteobacteria (09 genera) and Betaproteobacteria (03 genera) (Table S3). In silico taxonomic assignments of the 77 ASVs revealed that many families in the class Alphaproteobacteria such as Bradyrhizobiaceae (renamed Nitrobacteraceae), Methylobacteriaceae or Phyllobacteriaceae are not accurately assigned in SILVA 138 database, as reported elsewhere [114]. Indeed, these three families were misidentified, including Xanthobacteraceae, Beijerinckiaceae and Rhizobiaceae, respectively (Table S2). These weaknesses were compensated using the validly published names reported on the LPSN website [101]. Beta-proteobacteria genera harboring described rhizobia species. Evolutionary relationships were inferred using the Maximum Likelihood method based on the Tamura 3-parameter using a discrete Gamma distribution with invariant sites (T92+G+I). Bootstrap values ≥ 70% based on 1000 replicates are indicated, and the scale bar represents the number of substitutions per site. Type species of the 18 genera are displayed with strain ID followed by the GenBank 16S rRNA gene accession number. DNA sequences for ASVs used in this tree are provided in Table S3, and they are taken from the complete sequencing data archived in the NCBI SRA database. All ASVs that could not be accurately identified in the tree are enclosed in quotation marks.

Richness and Diversity Indices
The seven soils had a comparable number of ASVs (ASVs richness) which ranged from 24 (soil # CI11 and # CI17) to 35 (# CI20) ( Table S5). The indices of richness (Chao1) and alpha diversity (Shannon, Simpson and Fisher) analyzed per savannah biome are similar among the Sudanian savannah and the Sub-Sudanian savannah (Table S6). The community diversity indices showed that the sites CI18 and CI44 shared the lowest value of Bray-Curtis dissimilarity (calculated value of 0.29), meaning that these two sites shared the highest number of ASVs together when the composition of all the seven sites was compared. In contrast, CI14 and CI17 had the lowest number of shared ASVs (Table 4).  The Jaccard distance from the community diversity analysis revealed that the sites CI11 and CI17 were the most dissimilar (Jaccard distance of 0.82) while CI44 and CI20 were the least dissimilar among all the seven sites (Jaccard distance of 0.52) ( Table 5). The canonical correspondence analysis showed that the pH, C, Ca 2+ , K + , Mg 2+ and Na + were the soil properties that most strongly influenced the distribution of rhizobial taxa from the savannah soils in Northern Côte d'Ivoire ( Figure 6).

In Silico Evaluation of 16S rDNA V-Regions Discriminatory Power for Rhizobia
In an attempt to assign the 77 rhizobial ASVs detected from the savannah soils in Northern Côte d'Ivoire to the 18 genera of rhizobia validly published to date, we found that some genera type species, including Aminobacter aminovorans DSM 7048 T and Mesorhizobium loti DSM 2626 T , had identical 16S rRNA gene V4-V5 region ( Figure 2). Thus, we struggled to cluster at the genus level all the ASVs that belong to the Aminobacter -Mesorhizobium clade, including ASV_185, ASV_252 and ASV_1356. These weaknesses were compensated using several approaches, including BlastN with NCBI/GenBanK online nr/nt and wgs databases. However, the taxonomic affiliation of nine ASVs that belonged to the clades of Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium and Burkholderia-Caballeronia-Parabukholderia was partially solved (Figure 2, Table S2). In silico comparisons of different 16S rRNA gene variable regions to assess their effectiveness for differentiating rhizobia confirmed that the V4-V5 region has an insufficient resolution for separation of all genera of rhizobia at the genus level. As expected [115], the V4 region alone did also not perform well, regardless of the threshold used for the delineation. Both V4 and V4-V5 primer pairs were the only sets of primers that could not discriminate all 18 genera at the ASV level (one nucleotide polymorphism level). In contrast, the analysis showed that the V5-V7 region was the best target for genus discrimination (Figure 7): the V5-V7 region discriminated all the 18 genera of rhizobia at 100% and 99% thresholds, and 16 genera at the 97% threshold (Figure 7). It is noted that the two genera which were not discriminated by the V5-V7 region at the 97% threshold were Aminobacter and Mesorhizobium, suggesting that a threshold higher than the classical 97% should be used when targeting these two genera in HTAS analyses.

Discussion
A few studies reporting a characterization of soil rhizobial communities using the HTAS of 16S rDNA variable regions were carried out in temperate arable soils in East Europe [45] and coniferous forest soils in North America [41]. Up to now, little is known of the soil microbiome of semiarid areas commonly known as savannahs [116,117], although there are considered important biodiversity hotspots, including for microorganisms [118,119]. Two examples of savannah biomes have been neglected in terms of research in microbial ecology for decades: the Brazilian Cerrado savannah and the African savannah [116,120,121]. The microbiome of the Cerrado savannah is relatively more explored, including for Archaea [122], Bacteria [118,123], Fungi [124] and Protists [125], unlike that of the African savannah, which has not been studied in a systematic manner [35,48,126]. The current study is among the pioneer studies on African savannahs microbiome [35,127], and it provides new insights into the presence and distribution of taxa of rhizobia across the Sudanian and the Sub-Sudanian savannah zones. It revealed that the rhizobial diversity in the savannah zones in Northern Côte d'Ivoire is considerable in terms of richness and relative abundance of genera and families detected. These findings are similar to those observed in the Brazilian Cerrado savannah [118,123] and the African Miombo Woodlands in Mozambique [127], where rhizobacteria, including rhizobia, were found genetically diversified and abundant [118,127]. However, these results contrast with a similar study carried out in the Mopane woodlands, another important savannah ecosystem in southern Africa [35].
Of the 18 Alphaand Beta-proteobacterial genera harboring the described rhizobial species, only Allorhizobium, Aminobacter, Azorhizobium, Ochrobactrum, Pararhizobium and Phyllobacterium were not detected in soils from the savannah zones in Northern Côte d'Ivoire. As these six genera have also a low relative number of validly published species to date, all these data suggested that they are probably less abundant and diversified in soils and/or are associated with a limited set of legumes species. From all the 12 genera detected, Bradyrhizobium was found more abundant and ubiquitous, together with Microvirga and Paraburkholderia. Contrasting findings have been reported on the prevalence, genetic diversity and the ubiquity of these three genera of rhizobia. An HTAS study of the potential nitrogen-fixing bacteria in Polish soils detected Devosia, Mesorhizobium, Methylobacterium, Microvirga, Phyllobacterium, and Rhizobium (alpha-rhizobia), as well as Burkholderia sensu lato (s.l.) and Cupriavidus (beta-rhizobia), but noted the absence of Bradyrhizobium [45]. In contrast to Wolińska et al. [45], a recent atlas established for dominant soil bacteria classified Bradyrhizobium and Devosia among the most abundant and ubiquitous bacteria worldwide, with an apparent paucity of Burkholderia s.l in soils [6]. Nevertheless, a survey of the top 20 most abundant genera found in soil samples revealed that Bradyrhizobium and Burkholderia s.l. are, respectively, the first and the second most prevalent genera of soil bacteria [128]. Although all these data indicated that Bradyrhizobium and/or Burkholderia s.l. and/or Microvirga were frequently detected among the most dominant bacteria genera in soil samples, to our knowledge, the current study is the first showing that the Bradyrhizobium genus dominates in tropical savannah soils, together with Microvirga and Paraburkholderia. The predominance and the ubiquity of rhizobia genera, including Bradyrhizobium and Burkholderia s.l., is thought to be due to their genetic diversity, and their catabolic versatility that enables them to degrade recalcitrant compounds and survive in oligotrophic environments [41,129,130]. Since Moulin et al. [131] described two Burkholderia nodule-forming strains isolated in French Guiana and in South Africa, beta-rhizobia have been routinely identified from soils, mainly in South Africa, South America and southeast Asia [40]. Some of these studies even reported the dominance of Paraburkholderia when compared to cosmopolitan Bradyrhizobium in several soils, depending on the biome (e.g., the Cerrado, Caatinga and Forest Atlantic biomes in Brazil; the Fynbos biome in South Africa), the legumes species (Mimosa spp.; Lebeckia spp.) and the soil types [39,40,104,132]. Several studies have demonstrated that the beta-rhizobia are well adapted to poor and acidic soils [37,133,134]. Our study suggests that Paraburkholderia and Trinickia are more abundant in the mildly acidic soils (pH 5.7 < pH< 6.0), all of which harbored anthropogenic activities (fields of cashew and cereals etc.). Despite this observation in the cultivated soils, the impact of the savannah types on the dynamics of rhizobia diversity and abundance was not established in this study.
Although the variable regions of the 16S rRNA gene (e.g., individual V-regions, adjacent V-regions, pairs of non-contiguous V-regions) are well known [106,135,136], the selection of the most efficient variable region (s) for microbiome analysis is still debated [63,106,107,110,[137][138][139]. Many studies indicated that the efficiency of the variable regions for HTAS analysis depends on multiple parameters, including the microorganisms of interest and the extent to which their 16S rRNA genes have evolved [105,140,141]. For rhizobia, our study suggested that the V5-V7 region could be suitable for differentiating strains at the genus level, possibly replacing the use of the V4-V5 region. In a previous study, Eardly et al. [142] identified the V7 region alone as highly polymorphic in the Rhizobiales. Taken together, we suggest that the V5-V7 region contains sufficiently polymorphic DNA sequences to resolve the genetic complexity of the full 16S rRNA gene in rhizobia.
Many studies had reported the use of single-copy housekeeping genes in microbiome analyses to improve resolution at species and subspecies levels [143][144][145][146][147]. A multigenic approach that includes at least one housekeeping gene (e.g., rpoB) and one variable region of the 16S rRNA gene [148] is also considered a promising methodology. Taking into account these recommendations, we further propose the use of the V5-V7 region to analyze the rhizobial microbiome in combination with one of the four housekeeping genes (atpD-gyrB-recA-rpoB) that have been used for resolving ambiguous cases of identification among Rhizobium strains [149].

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/microorganisms9091842/s1, Figure S1: Covering area of each sampling site and its corresponding features, Figure S2: Rarefaction curve of the seven samples, indicated by the number of ASVs depending on the size of the sequence sample, Figure S3: 16S rDNA V4-V5 phylogenetic tree showing the relationship between all selected reference alpha and beta-rhizobia strains and the 77 ASVs detected in savannah soils of Northern Côte d'Ivoire, Table S1: Physico-chemical properties of samples soils, Table S2: Levels of similarity between the V4-V5 sequences of the 77 ASVs of rhizobia and all 18 alphaproteobacterial and betaproteobacterial genera harboring rhizobia strains, Table S3: Major characteristics of the 77 ASVs of rhizobia detected in savannah soils of Northern Côte d'Ivoire, Table S4: ASVs and their corresponding md5_hash identities, Table S5: ASVs richness, cumulative relative abundance of sequences and prevalence of ASVs per soil, Table S6: Measure of the richness and the alpha diversity per savannah zone.

Data Availability Statement:
The high-throughput amplicon sequencing reads of the 16S rRNA gene used in this study are available in the NCBI SRA database under accession numbers SRR13623317, SRR13623319, SRR13623320, SRR13623323, SRR13623324, SRR13623326 and SRR13623335.