Rhizosphere Diazotrophs and Other Bacteria Associated with Native and Encroaching Legumes in the Succulent Karoo Biome in South Africa

Total and diazotrophic bacteria were assessed in the rhizosphere soils of native and encroaching legumes growing in the Succulent Karoo Biome (SKB), South Africa. These were Calobota sericea, Lessertia diffusa, Vachellia karroo, and Wiborgia monoptera, of Fabaceae family near Springbok (Northern Cape Province) and neighboring refugia of the Fynbos biome for C. sericea for comparison purposes. Metabarcoding approach using 16S rRNA gene revealed Actinobacteria (26.7%), Proteobacteria (23.6%), Planctomycetes, and Acidobacteria (10%), while the nifH gene revealed Proteobacteria (70.3%) and Cyanobacteria (29.5%) of the total sequences recovered as the dominant phyla. Some of the diazotrophs measured were assigned to families; Phyllobacteriaceae (39%) and Nostocaceae (24.4%) (all legumes), Rhodospirillaceae (7.9%), Bradyrhizobiaceae (4.6%) and Methylobacteriaceae (3%) (C. sericea, V. karroo, W. monoptera), Rhizobiaceae (4.2%; C. sericea, L. diffusa, V. Karroo), Microchaetaceae (4%; W. monoptera, V. karroo), Scytonemataceae (3.1%; L. diffusa, W. monoptera), and Pseudomonadaceae (2.7%; V. karroo) of the total sequences recovered. These families have the potential to fix the atmospheric nitrogen. While some diazotrophs were specific or shared across several legumes, a member of Mesorhizobium species was common in all rhizosphere soils considered. V. karroo had statistically significantly higher Alpha and distinct Beta-diversity values, than other legumes, supporting its influence on soil microbes. Overall, this work showed diverse bacteria that support plant life in harsh environments such as the SKB, and shows how they are influenced by legumes.


Introduction
The Succulent Karoo Biome (SKB) in South Africa is one of the world's most species rich biomes, with over 5000 plant taxa of which 40% are endemic to the region [1][2][3]. The main vegetation cover includes leaf-succulent, stem-succulent and non-succulent plant species [2,4]. The SKB is semi-arid, experiences winter rainfall, and the high levels of plant diversity have been ascribed to the region's oscillating wet and dry climatic conditions, moderate climatic history and soils that are sandy, acidic and nutrient limited [4,5].
The need for conservation of this ecosystem has directed some research on the influence of rainfall and droughts on vegetation cover [3,6]. However, the fact that the success and survival of plant species in this relatively harsh environment are due to their interactions with microbial communities in the soil, is often underestimated [7,8].
Members of the Fabaceae family form a significant part of the non-succulent plant species in the SKB, particularly due to their ability to survive in extreme environments [9][10][11][12]. These plants generally possess the ability to fix atmospheric nitrogen in association with symbiotic and free-living soil bacteria [13]. For example, diverse species of Ensifer were reported to be the nodulating symbionts of Vachellia jacquemontii, thus promoting the The SKB of South Africa is situated in the western part of the Northern Cape Province. It lies below 800 m but may reach up to 1500 m above the sea level (a.s.l) on the eastern side. It is characterized by semi-arid to arid conditions, receiving 20 to 290 mm of rainfall per year during the winter season. The temperature may rise beyond 40 • C during the summer periods. SKB soils are generally acidic, weakly developed and mostly found on rocks with limited nutrients [2,5,29]. Based on the occurrence of legumes, two sites within this biome were identified, which were Kamieskroon (30 • 38.67 E) were also selected ( Figure 1). The Kamiesberg Center is an outlying area of the Fynbos, mostly above 1200 m a.s.l., with climatic conditions and soil properties that are generally similar to those of the SKB [30,31].
Center is an outlying area of the Fynbos, mostly above 1200 m a.s.l., with climatic conditions and soil properties that are generally similar to those of the SKB [30,31].

Soil Sampling and Physico-Chemical Properties
In total, 24 composite soil samples were obtained from the four sites under the four different legumes. Five soil sub-samples were randomly scooped at a depth of 0-10 cm using a soil auger in the root zone area under each legume. The five sub-samples were bulked to form one composite sample per plant. Such samples were collected for C. sericea from Brakputs (n = 3), Kamiesberg (n = 4), Kamieskroon (n = 2) and Leliefontein (n = 1), while those for V. karroo were collected from Brakputs (n = 5), and Kamieskroon (n = 1). For W. monoptera from Brakputs (n = 3) and L. diffusa from Brakputs (n = 5). The composite soil samples were stored at 4 °C. Prior to analysis, each composite soil sample was split into two proportions. The first proportion was sieved (2 mm mesh) and used for soil properties analysis of pH using the KCl method, total carbon (TC), and total nitrogen (TN) using the dry combustion method, ammonium (NH4 + ), nitrates (NO3 -) measured calorimetrically and soil texture determined using the Hydrometer method at Bemlab, Cape Town, South Africa. The second proportion was kept at −70 °C for subsequent analyses.

Soil Sampling and Physico-Chemical Properties
In total, 24 composite soil samples were obtained from the four sites under the four different legumes. Five soil sub-samples were randomly scooped at a depth of 0-10 cm using a soil auger in the root zone area under each legume. The five sub-samples were bulked to form one composite sample per plant. Such samples were collected for C. sericea from Brakputs (n = 3), Kamiesberg (n = 4), Kamieskroon (n = 2) and Leliefontein (n = 1), while those for V. karroo were collected from Brakputs (n = 5), and Kamieskroon (n = 1). For W. monoptera from Brakputs (n = 3) and L. diffusa from Brakputs (n = 5). The composite soil samples were stored at 4 • C. Prior to analysis, each composite soil sample was split into two proportions. The first proportion was sieved (2 mm mesh) and used for soil properties analysis of pH using the KCl method, total carbon (TC), and total nitrogen (TN) using the dry combustion method, ammonium (NH 4 + ), nitrates (NO 3 − ) measured calorimetrically and soil texture determined using the Hydrometer method at Bemlab, Cape Town, South Africa. The second proportion was kept at −70 • C for subsequent analyses.

Soil Metabarcoding Analysis
A culture-independent (metabarcoding) approach based on next generation sequencing of specific marker genes was used to investigate the diversity and abundance of total and diazotrophic bacteria in the collected soils [32]. Total soil genomic DNA was isolated from 0.5 g of each composite soil sample using the FastDNA TM SPIN Kit for Soil (MP Biomedicals, Solon, OH, USA). The quality and quantity of extracted DNA were evaluated using 1.0% (w/v) agarose (Whitehead Scientific, Cape Town, South Africa) gel electrophore-sis and a NanoDrop™ 2000/2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). For each soil DNA sample, two consistent NanoDrop readings, whose variation was less than 5%, were considered for subsequent determination of the DNA concentration.
The two metabarcoding markers used were the universal house-keeping 16S rRNA gene and the nifH gene specific to diazotrophs [33]. The extracted DNA was split into two proportions. The first proportion was sent to the Biomedical Core Research Facilities at the University of Michigan, Ann Arbor, USA for sequencing of the V4 hypervariable region of the 16S rRNA gene. For library preparation, forward and reverse 515F/806R primers contained Illumina-sequencing P5 and P7 dual-indexing adapters with a unique barcode for every sample for multiplexing and demultiplexing [34]. PCR conditions were as follows: 94 • C for 3 min, followed by 33 cycles of 94 • C for 30 s, 55 • C for 30 s and 72 • C for 30 s with a final elongation step at 72 • C for 5 min.
The second DNA proportion was sent to MR DNA (www.mrdnalab.com, accessed on 2 October 2021, Shallowater, TX, USA) for sequencing of the nifH gene using PolF/PolR primers [35]. For these reactions, the forward primer contained a barcode. The amplifications were carried out using the HotStarTaq Plus Master Mix Kit (Qiagen, Germantown, MD, USA) with cycling conditions as follows: 94 • C for 3 min, followed by 28 cycles of 94 • C for 30 s, 53 • C for 40 s, and 72 • C for 1 min and a final elongation step at 72 • C for 5 min.
For both genes, amplicons from different samples were pooled together in equal concentrations and purified using calibrated Agencourt Ampure XP beads (Agencourt Bioscience Corporation, Beverly, MA, USA). Paired-end sequencing was performed using the Illumina MiSeq platform to produce paired reads that were 250 nucleotides in length. All raw sequence data for the 16S rRNA and nifH genes have been deposited with links to BioProject accession numbers (PRJNA767582, PRJNA767584), respectively.
The 16S rRNA sequence data were processed and analysed using Mothur version 1.35.1 [36]. After joining the paired forward and reverse sequence reads, the sequence quality was checked by specifying a maximum length of 275 base pair (bp) and a minimum of 250 bp, as well as removal of all sequences with ambiguities and with more than 8 bp homopolymers. The sequences were then sorted into their respective samples by allowing of maximum of 2 bp mismatches with the primer and no mismatches in the multiplexing barcodes. The quality-filtered and sorted sequences were then aligned to the SILVA ribosomal RNA reference database [37] provided through Mothur and filtered (filter.seqs option in mothur) to ensure that all the sequences aligned to a similar region of the V4 region of the 16S rRNA gene. These sequences were then processed using a pseudo-single linkage algorithm with the pre.cluster command in Mothur using a 2 bp similarity cut-off [36,38]. In addition, chimeras were detected and removed [39] using the chimera.uchime command in Mothur. The remaining sequences were then associated with particular taxonomic lineages using the Greengenes database [40] with a cut-off confidence level of 80%. After removing sequences of unknown taxonomic lineages, the remaining data were classified into operational taxonomic units (OTUs) using average neighbour algorithm with a similarity cut-off of 97%. OTUs were assigned by obtaining a representative sequence from each OTU (get.oturep command in Mothur) and classifying it using the Greengenes database at a confidence level of 80% [40].
The nifH sequence data were processed using the analysis pipeline of MR DNA (www.mrdnalab.com, accessed on 2 October 2021, Shallowater, TX, USA) according to qiime pipeline using the divisive amplicon denoising algorithm (dada2) [41]. This entailed the joining of paired end sequences and quality filtering. For the latter, sequences shorter than 200 bp, containing more than two ambiguous characters, that had homopolymer runs longer than 6 bp, or that contained more than one mismatch to the sample-specific barcode or to the primers sequences were removed. Following removal of the barcode and primer sequences, the raw nifH reads were denoised to remove sequencing errors, separated into OTUs defined by clustering at 97% sequence similarity [42] using a curated and publicly available database of nifH sequences (https://wwwzehr.pmc.ucsc.edu/, accessed on 15 September 2021) [43], after which chimeras were removed [44][45][46]. Final OTUs were taxonomically classified using BLASTn against this database.

Statistical Analyses of Metabarcoding Data and Soil Physico-Chemical Properties
Prior to statistical analyses, singletons were removed using the filter.shared command in Mothur (version 1.35.1) [36]. Alpha and beta diversity analyses were conducted to compare samples based on OTU assignment. Alpha diversity indices (richness, non-parametric Shannon diversity and Shannon evenness) [47] were calculated using the summary.single command in Mothur [36]. Beta diversity OTU-based (Bray-Curtis and Jaccard) matrices to compare the structure (Bray-Curtis) and membership (Jaccard) of the samples were calculated using the dist.shared command in Mothur [36]. The least number of sequences across all the samples after removal of singletons were set as the sampling depths (i.e., n = 1465 for 16S rRNA and n = 2691 for nifH) and subsampled 1000 times (iters = 1000) to standardize both the alpha and beta diversity calculations. Good coverage estimate for the alpha diversity matrices was used to determine the % coverage within samples after sub-sampling.
Prior to analysis of variance (ANOVA) or permutations analysis of variance (PER-MANOVA), a legume rhizosphere single soil sample from Leliefontein was eliminated. Therefore, four legume species (C. sericea, L. diffusa, V. karroo, and W. monoptera from three sampling sites (Brakputs, Kamiesberg, and Kamieskroon) totalling to 23 samples were considered. Alpha diversity and soil properties data were subjected to normality tests using histograms and Q-Q plots to check the homogeneity of variance in R (v.4.0.0) software (https://www.rstudio.com, accessed on 5 September 2021). One-way ANOVA was conducted to determine whether the Alpha diversity indices and soil properties were significantly different under legumes and sites where the samples were collected. The Akaike information criterion (AIC) was used to select the best-fit model to explain variation in the dependent variables. Table of p values for all possible pairwise comparisons of least square means (lsmeans), and the letters (letter display) were conducted using the Tukey's honestly significant difference (Tukey's HSD) post-hoc test in R (v.4.0.0). Pearson correlation and regression analysis with the soil properties were also conducted using R (v.4.0.0) to test if these data were correlated and by how much the properties explained the changes in the alpha diversity indices. PERMANOVA was conducted to test if legume species and sites influenced the Beta diversity matrices.
For graphical visualization of the structure of communities as shaped by the legumes in the respective sampling sites, non-metric multidimensional scaling (NMDS) and principal co-ordinate (PCoA) were performed in R (v.4.0.0) on OTU-based matrices (Bray-Curtis and Jaccard) [48]. Canonical analysis of principal coordinates (CAP) was also conducted in R (v.4.0.0) to emphasize the influence of soil properties on the structure of bacteria [49]. Heatmaps were generated using relative abundance data of OTUs using the R (v.4.0.0) package heatmap.

Physico-Chemical Properties of the Studied Rhizosphere Soils from the Succulent Karoo Biome
The results of the soil properties analysis are presented in Table SI. Overall, soil properties were not significantly different among the sampling sites (Brakputs, Kamiesberg and Kamieskroon) (p > 0.05). Legume species, however, revealed significant differences for pH and TC (p < 0.05). Overall, pH values ranged from 4.00 (W. monoptera) to 6.47 (V. karroo). The pH values associated with V. karroo were significantly higher than those from W. monoptera and C. sericea (p < 0.05). The pH therefore ranged from extremely acidic for W. monoptera, to strongly acidic for C. sericea, L. diffusa and slightly acidic for V. karroo. TC values ranged from 0.61 to 1.93 (%) for L. diffusa and V. karroo, respectively. TC values under V. karroo were significantly higher than those from C. sericea and L. diffusa (p < 0.05). TC values were generally low for all the species, except V. karroo where it was high. TN values ranged from 0.05 to 0.16 (%) for rhizosphere soils from L. diffusa and V. karroo, respectively.
Generally, TN values ranged from low to very low for all the legume rhizosphere soils analyzed. The NH 4 + and NO 3 − content ranged from 9.46 to 17.1 mg-N kg −1 and 0.57 to 10.9 mg-N kg −1 , respectively. There were no significant differences observed for TN, NH 4 + and NO 3 − values between the different legumes, although V. karroo consistently maintained the highest values (p > 0.05).
Based on the nifH data, several diazotrophic families were measured ( Figure 3B; Table S3). These consisted of OTUs from the families Nostoceae, Microchaetaceae and Scytonemataceae (24.4%, 4.0%, and 3.1%), respectively, (Cyanobacteria) of the total sequences ( Figure 3B; Table S3). Others OTUs belonged to the families Phyllobacteriaceae (39.0%), and Enterobacteriaceae (0.06%) (Proteobacteria) of the total sequences. There was also one OTU belonging to the family Lachnospiraceae (0.2%) (Firmicutes) of the total sequences (Table S3). The high dominance of Phyllobacteriaceae was also supported by the recovery of one OTU in all 24 soil samples examined. This OTU represented a member of the genus Mesorhizobium and formed 27.3% of the total nifH sequences. Two OTUs, members of Proteobacteria representing approximately 27.3% of the total sequences were shared among the SKB samples.

OTU-Based Diversity of Diazotrophs as Influenced by Legume Species
Legumes showed common and distinct influences on the diazotrophs at family level in the arid soils of the Succulent Karoo or Fynbos biomes based on nifH metabarcoding analysis ( Figure 3B; Table S3). For example, the rhizosphere soils of C. sericea and V. karroo shared most members of the identified families. These included OTUs from Bradyrhizobiaceae, Acetobacteraceae, Phyllobacteriaceae, Rhodospirillaceae, Comamonadaceae, Rhizobiaceae, Rhodocyclaceae, Nostocaceae, Methylobacteriaceae, and Sphingomonadaceae. On the other hand, L. diffusa and W. monoptera showed few families whose OTUs were mainly shared either with C. sericea, V. karroo or both. These were OTUs from Bradyrhizobiaceae, Nostocaceae, Microchaetaceae, Rhodospirillaceae, Methylobacteriaceae, Acetobacteraceae, and Scytonemataceae in the rhizosphere soils of W. monoptera and OTUs from Phyllobacteriaceae, Rhizobiaceae, Nostocaceae, and Scytonemataceae in the rhizosphere soils of L. diffusa (Table S3).
In terms of the association of diazotrophs with specific legumes, the rhizosphere soils of V. karroo specifically showed a high diversity and abundance of OTUs representing Alcaligenaceae, Desulfovibrionaceae, Geobacteraceae, Lachnospiraceae, Ectothiorhodospiraceae and Pseudomonadaceae among others (Table S3). Enterobacteriaceae was specifically found in the rhizosphere of W. monoptera (Table S3).

Alpha Diversity
The richness (observed taxa) and structure-based alpha diversity metrics (Shannon Index and Shannoneven Index) were strongly influenced by legumes species (p < 0.01) but not the sampling sites or the interactions between legume species and sites (p > 0.05) for both 16S rRNA and nifH genes (Table S4). Therefore, only the results of the influence of legume species on alpha diversity indices are presented.

Beta-Diversity
Based on 16S rRNA gene, the average dissimilarity in community structure between legumes and sites was 76 ± 9% (Bray-Curtis), and 85 ± 5% (Jaccard) (Figure 4A,B). These high indices suggested a low level of similarity between the samples examined. Prior to PERMANOVA analysis, ordination analysis (NMDS and PCoA) based on Bray-Curtis and Jaccard matrices revealed legumes specific clustering of bacteria particularly for V. karroo (results not shown). PERMANOVA analysis based on both matrices revealed that variability among samples was best explained by the legumes (p < 0.001) (Table S4). CAP analysis showed pH, TN, and TC as important factors influencing the 16S rRNA Betadiversity ( Figure 5A,B). For example, communities clustered together in the direction of high pH, TN and TC particularly those associated with V. karroo. Abbreviations: TN, total nitrogen; TC, total carbon; NH4 + , ammonium; NO3 -, nitrates. Significance levels: ns, not significant at p > 0.05, significant at p < 0.05 *; p < 0.01 **, and p < 0.001 *** levels.

Beta-Diversity
Based on 16S rRNA gene, the average dissimilarity in community structure between legumes and sites was 76 ± 9% (Bray-Curtis), and 85 ± 5% (Jaccard) (Figure 4A,B). These high indices suggested a low level of similarity between the samples examined. Prior to PERMANOVA analysis, ordination analysis (NMDS and PCoA) based on Bray-Curtis and Jaccard matrices revealed legumes specific clustering of bacteria particularly for V. karroo (results not shown). PERMANOVA analysis based on both matrices revealed that variability among samples was best explained by the legumes (p < 0.001) (Table S4). CAP analysis showed pH, TN, and TC as important factors influencing the 16S rRNA Beta-diversity ( Figure 5A,B). For example, communities clustered together in the direction of high pH, TN and TC particularly those associated with V. karroo. The average dissimilarity for the nifH gene was 91 ± 13% (Bray-Curtis) and 89 ± 6% (Jaccard) (Figure 4C,D). These indicated a low level of similarity between the samples. Similar to the 16S rRNA, nifH gene-based ordination analysis on diazotrophic communities revealed legumes specific clustering particularly for V. karroo ( Figure 5C,D). PER-MANOVA analysis based on both matrices revealed that legumes and sites influenced structure and membership matrices (p < 0.01) (Table S4). For example, CAP analysis showed clear clustering patterns among legumes (e.g., L. diffusa, W. monoptera and V. karroo). In addition, nitrates and TN were found as the environmental properties that mainly influenced the nifH Beta diversity, as communities clustered in the direction of high ni-

Bacteria Composition in Rhizosphere Soils Used Is Associated with Environmental Conditions
Based on the metabarcoding of the 16S rRNA, several of the taxa identified (e.g., Proteobacteria, Actinobacteria, Bacteroidetes, and Cyanobacteria) have previously been shown to dominate rhizosphere and desert soils [7,50]. Most of them are reported to have inherent properties making them well suited to the SKB soils, such as adaptation to high temperatures and soil acidity as is the case for the SKB. For instance, many species of Cyanobacteria, Planctomycetes and Verrucomicrobia are known to thrive in higher temperatures and are tolerant to low soil pH [51][52][53].
The consistent dominance of the Proteobacteria phylum using both 16S rRNA and nifH genes in this study was in agreement with several other studies which reported their abundance in different soil ecosystems such as agricultural, forests, grasslands, saline as well as semi-arid soils [54]. This is attributed to their role in sustaining these environments through various biogeochemical processes. Cyanobacteria members are aquatic, but have also consistently been reported in soil crusts and rhizosphere soils in arid and semi-arid environments, particularly due to their role in atmospheric nitrogen fixation [51,54]. This was consistent with the nifH findings in this study. Their low measurements according to The average dissimilarity for the nifH gene was 91 ± 13% (Bray-Curtis) and 89 ± 6% (Jaccard) (Figure 4C,D). These indicated a low level of similarity between the samples. Similar to the 16S rRNA, nifH gene-based ordination analysis on diazotrophic communities revealed legumes specific clustering particularly for V. karroo ( Figure 5C,D). PERMANOVA analysis based on both matrices revealed that legumes and sites influenced structure and membership matrices (p < 0.01) (Table S4). For example, CAP analysis showed clear clustering patterns among legumes (e.g., L. diffusa, W. monoptera and V. karroo). In addition, nitrates and TN were found as the environmental properties that mainly influenced the nifH Beta diversity, as communities clustered in the direction of high nitrates and TN particularly those associated with V. karroo and some C. sericea species (Figure 5C,D).

Bacteria Composition in Rhizosphere Soils Used Is Associated with Environmental Conditions
Based on the metabarcoding of the 16S rRNA, several of the taxa identified (e.g., Proteobacteria, Actinobacteria, Bacteroidetes, and Cyanobacteria) have previously been shown to dominate rhizosphere and desert soils [7,50]. Most of them are reported to have inherent properties making them well suited to the SKB soils, such as adaptation to high tempera-tures and soil acidity as is the case for the SKB. For instance, many species of Cyanobacteria, Planctomycetes and Verrucomicrobia are known to thrive in higher temperatures and are tolerant to low soil pH [51][52][53].
The consistent dominance of the Proteobacteria phylum using both 16S rRNA and nifH genes in this study was in agreement with several other studies which reported their abundance in different soil ecosystems such as agricultural, forests, grasslands, saline as well as semi-arid soils [54]. This is attributed to their role in sustaining these environments through various biogeochemical processes. Cyanobacteria members are aquatic, but have also consistently been reported in soil crusts and rhizosphere soils in arid and semi-arid environments, particularly due to their role in atmospheric nitrogen fixation [51,54]. This was consistent with the nifH findings in this study. Their low measurements according to the universal 16S rRNA gene could be due to its low-resolution power compared to functional and group specific primers such as the nifH gene [55]. Firmicutes on the other hand, were less diverse and abundant in the considered rhizosphere soils using the both metabarcoding approaches. This was not surprising as members of this phylum are mainly associated with contaminated environments such as sludges and soil impacted by acid mine drainage [54,56].
Environmental factors strongly influenced the total and diazotrophic bacteria in the soils examined. This was particularly evident for soil pH and nutrients. Significantly positive correlations of >70% for pH and >40% for TN and TC, were observed with the richness and structure of total bacteria. These factors accounted for much of the variation within communities (at least 45% for pH, 20% for TN and at least 15% for TC). These factors further influenced the diazotrophs, as significant positive correlations of >40% were observed for species richness for all the soil properties except pH. TC and TN influenced the structure of the diazotrophs as >55% positive correlations were observed for the species structure, while TC showed 45% positive correlations with species evenness. These factors accounted for at least 18% of variations within the diazotrophs. The influence of environmental factors on microbes was also revealed for beta-diversity analyses. Similar effects of pH, TN, TC, and NO 3 − on soil bacteria were found in the bulk or rhizosphere soils of Acacia dealbata in the grassland biome in South Africa [57], and in an agricultural soil in Kenya [58]. Soil pH has also been reported as a main driver of changes in the structure and composition of total soil bacteria [58,59].
It is likely that the different groups of bacteria identified in this study form part of specialized communities, some of which interact with plants to facilitate or enhance their ability to colonize and become established in the harsh environmental conditions of the SKB. Ways in which this facilitation can occur is through nitrogen, carbon and phosphorus compounds cycling, thereby improving the plants' access to nutrients [60,61]. Many species of the Gemmatimonadetes, Cyanobacteria, Proteobacteria, Chloroflexi, Firmicutes, and Acidobacteria, for example, are known to facilitate carbon assimilation via (bacterio) chlorophyll-based photosynthesis [62,63], while certain Proteobacteria and Firmicutes are capable of phosphorus solubilization, thus making this element available to plants [64,65]. However, in the nitrogen-poor soils of the SKB, cycling of this element is perhaps equally important, and various members of the Proteobacteria [66,67], Planctomycetes [53], Bacteroidetes [68], Verrucomicrobia [69], Firmicutes [64], and Cyanobacteria [51,52] are known to be capable of fixing atmospheric nitrogen to provide ammonium and/or nitrates for plant growth.
The bacteria in the rhizosphere soils of the legumes investigated are involved in the promotion of the plants' growth and might even protect them against pathogens. Many Acidobacteria, Firmicutes and Proteobacteria are well known to produce indole-3acetic acids (IAAs), siderophores and the enzyme aminocyclopropane carboxylate (ACC) deaminase [64,70]. IAA is a plant growth hormone [13,64], and siderophores chelate iron and avail it to plants for growth, while making it unavailable to other organisms, particularly pathogens [64,71]. Under stressful conditions (e.g., drought and heat), the ACC deaminase promotes root growth by lowering ethylene levels in plants [68,[71][72][73].
Taken together, the diverse bacteria measured thus underscore the crucial role they play in supporting plants in the SKB.

Legume Species Influence Diazotrophs in the Succulent Karoo Biome
Metabarcoding with the nifH allowed a more in-depth analysis of the diazotrophic taxa occurring in the rhizosphere soils examined. Our findings revealed diverse OTUs which were separated into at least 19 families of diazotrophs. The distribution, diversity and abundance of these diazotrophs were however, modulated by legumes considered. For example, the rhizosphere soils of L. diffusa were mainly dominated by OTUs belonging to Phyllobacteriaceae, Rhizobiaceae, Nostocaceae, and Scytonemataceae. These families contain species known to fix atmospheric nitrogen. For example, Phyllobacteriaceae contained the rhizobial genus Mesorhizobium and Rhizobiaceae contained the genera Sinorhizobium/Ensifer among others. These findings support previous reports which demonstrated the symbiotic relationships of Mesorhizobium and Sinorhizobium/Ensifer with Lessertia species [17,66]. Nostocaceae and Scytonemataceae on the other hand, are known to contain species with free-living nitrogen fixation abilities [51]. The rhizosphere soils of W. monoptera mainly contained species belonging to Bradyrhizobiaceae, Methylobacteriaceae, Nostocaceae, Rhodospirillaceae, Acetobacteraceae, Microchaetaceae, Scytonemataceae, and Enterobacteriaceae families. Bradyrhizobiaceae and Methylobacteriaceae members have symbiotic rhizobial capabilities with these legumes, which is consistent with previous reports regarding the distribution of Bradyrhizobium and its association with a wide range of legumes [17,28,66]. Methylobacteriaceae members have also been reported to possess symbiotic associations with diverse legumes in the papilionoid tribe Crotalarieae [74,75]. Members from the other families contain free-nitrogen fixation capacities. For example, members from Rhodospirillaceae [76,77] and Microchaetaceae [51] as well as Enterobacteriaceae that occurred specifically in rhizosphere soils of W. monoptera [78]. Moreover, some diazotrophic members belonging to Acetobacteraceae promote plant growth potentially through phosphorus solubilization and antagonism against plant pathogens [65]. The rhizosphere soils of L. diffusa and W. monoptera therefore, contained diazotrophic communities with the ability to mainly fix atmospheric nitrogen and solubilize phosphorus for W. monoptera species. According to the alpha-diversity analyses, generally low species richness, structure and evenness was measured under the rhizospheres of W. monoptera and L. diffusa. Consistently, the communities under these two legumes clustered separately, being distinct from those under C. sericea and V. karroo. Both legumes are not widely distributed in the SKB, which explains the low diversity, abundance and specialization of their associated diazotrophs.
The rhizosphere of C. sericea was dominated by OTUs from Bradyrhizobiaceae, Rhodocyclaceae, Rhizobiaceae, Nostocaceae, Phyllobacteriaceae, Rhodospirillaceae, Acetobacteraceae, Sphingomonadaceae, Methylobacteriaceae, and Comamonadaceae families. Species from Bradyrhizobiaceae, Rhizobiaceae, Methylobacteriaceae and Phyllobacteriaceae have symbiotic characteristics with legumes in the papilionoid tribe Crotalarieae [66,[79][80][81]. Members of the Rhodocyclaceae, Nostocaceae, Rhodospirillaceae, Acetobacteraceae, and Comamonadaceae are free-living nitrogen fixers [82][83][84][85]. In addition, some members from Acetobacteraceae are known to promote plant growth through phosphorus solubilization and antagonism against pathogens [65]. Some members from Comamonadaceae promote plant growth through carbon cycling [86], while members from Sphingomonadaceae are known to promote plant growth through the production of IAAs [87]. These communities did not have specific clustering pattern with C. sericea according to the ordination analysis. Therefore, the diazotrophs associated with C. sericea are not only diverse and abundant, but also possess different plant growth promoting properties, ranging from atmospheric nitrogen fixation, phosphorus solubilization, production of IAAs and antagonism against pathogens. This could be the reason this legume is commonly distributed in the SKB and also in the neighboring Fynbos biome.
V. karroo species had higher diversity and abundance of diazotrophs than all the other legumes. Significantly higher species richness and diversity under V. karroo were observed using 16S rRNA-based metabarcoding. These were supported by ordination analysis, where beta-diversity indices showed clustering of members associated with V. karroo away from those associated with the other legumes. Additionally, higher pH, TC, TN, NH 4 + , and NO 3 − were measured in the V. karroo rhizosphere soils than for other legumes. Increased levels of nutrients and pH associated with V. karroo are believed to support its ability to attract diverse and abundant bacteria [57,59]. This gives it a competitive advantage in modifying the nutrient status of its rhizosphere soils more strongly than other legumes in the same habitat. In Acacia sensu lato, such advantages have been suggested to be as a result of nitrogen fixation, organic matter production via root exudates, and leave fall under their canopy [57]. These mechanisms drive the invasiveness of certain legumes [14,19,20] and may also play an important role in the encroaching behavior of V. karroo in the SKB.
The bacterial distribution patterns observed in the rhizosphere communities associated with V. karroo suggest that it modified its rhizosphere environment, favouring or enhancing the establishment of particular communities [19,20,24,88]. For example, OTUs from Geobacteraceae, Alcaligenaceae, Lachnospiraceae, Desulfovibrionaceae, Ectothiorhodospiraceae, and Pseudomonadaceae taxa were specifically associated with the rhizosphere soils of V. karroo. These diazotrophs are free-nitrogen fixers [61,[89][90][91][92][93][94][95]. Moreover, some members from Alcaligenaceae promote plant growth through production and control of IAAs levels [87,96]. OTUs such as those from Pseudomonadaceae promote plant growth through phosphorus solubilization, production of IAAs and siderophores [97][98][99]. Although this study is the first to report on the metabarcoding-based rhizosphere bacteria associated with V. karroo, different invasive trees such as Berberis thunbergii DC. (Japanese barberry) have been shown to attract different taxa including Pseudomonadaceae in their rhizosphere environment [88]. Acacia dealbata, an Australian invasive legume in South Africa, mainly attracted members from Bradyrhizobiaceae and Pseudomonadeceae among others in its rhizosphere [19,57]. Therefore, interactions between V. karroo and particular soil bacteria lead to significant changes in the composition of rhizosphere communities. Other free nitrogen-fixers which were associated with V. Karroo, but also found in the rhizosphere of the other legumes included members from Rhodocyclaceae, Nostocaceae, Rhodospirillaceae, Acetobacteraceae, Microchaetaceae, Sphingomonadaceae, and Comamonadaceae families [51,65,77,84,86,87]. In addition, some members from Acetobacteraceae, Sphingomonadaceae, and Comamonadaceae also promote plants growth through phosphorus solubilization, the production of IAAs, siderophore production for protection against pathogens, and carbon cycling [85,86,100].
With regards to possible rhizobial symbionts, V. karroo showed the ability to associate with diverse rhizobia. For example, several OTUs that were affiliated with Phyllobacteriaceae (Mesorhizobium), Bradyrhizobium, and a few with Rhizobiaceae (Rhizobium and Ensifer) formed part of its rhizosphere. These findings are consistent with other studies [14,17,66,101]. Undoubtedly, the ability of V. karroo to establish symbiotic interactions with a wide range of rhizobia also provides it with a competitive advantage over other legumes to colonize and encroach on new areas [25]. The high diversity of its potential rhizobial symbionts, combined with unique bacteria in V. karroo's rhizosphere soils in the SKB, support its involvement in reconstructing its microbiome, a strategy used by most invasive members of the mimosoid clade of legumes to invade and colonize new environments [19]. The rhizosphere of V. karroo therefore contained the most diverse and abundant diazotrophs that possess a wide range of plant growth promotion properties such as atmospheric nitrogen fixation, phosphorus solubilization, the production of IAAs and siderophores, carbon cycling as well as antagonism effects against pathogens.
Most studies that investigated the influence of plant species and environmental conditions on soil microbial communities have mainly focused on non-legumes versus legumes [55]. For example, legume species Stylosanthes guianensis (Aubl.) Sw., Trifolium pratense L., and Medicago sativa L. enriched soil microbial communities compared to grass species (Paspalum natatum, Festuca arundinacea L., Lolium perenne L.) in China [102]. In Nigeria, a legume species Pterocarpus erinaceus improved the soil physico-chemical properties and selected for ten bacterial species compared to a non-legume Anoigessus leiocarpa, which was associated with low values of soil properties and only six bacterial species [103]. This study further revealed that legume species, which fall under the same family Fabaceae, select for specific microbial communities and also share some microbial communities across their rhizospheres.

Conclusions
This study revealed diverse diazotrophs associated with legumes in the SKB in South Africa. Environmental factors such as soil pH, soil nutrients, and legume species influenced the microbial communities. A member of Mesorhizobium species was common in all rhizosphere soils considered. Other diazotrophs such as Bradyrhizobiaceae, Nostocaceae, among others were shared across several legume rhizosphere soils, while others such as Enterobacteriaceae and Pseudomonadaceae were specific to specific legume rhizosphere soils. Legumes such as L. diffusa and W. monoptera that are not widely distributed in the SKB haboured less diverse and abundant diazotrophs which mainly possess atmospheric nitrogen fixation properties compared to the widely distributed C. sericea and the encroaching V. karroo species. These legumes both contained diverse diazotrophic communities, associated with diverse plant growth promotion properties such as atmospheric nitrogen fixation, phosphorus solubilization, production of IAAs and siderophores, carbon cycling and plant protection against pathogens.
Since the metabarcoding approach of 16S rRNA and nifH genes only allowed us to identify communities that are potentially involved in important processes such as nitrogen fixation, metatranscriptomics analyses are recommended to study gene expression involved in crucial biological processes. Future research will also focus on root-bacteria nodulation and nitrogen-fixation with legumes and respective rhizosphere soils used in this study. If effective, such strains will be characterized and recommended for development as potential inoculum in the agricultural industry.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms10020216/s1, Table S1: Statistical analysis of soil chemical properties of the study sites showing main effects of different legume species and sites, Table S2: Diverse bacterial phyla in the Succulent Karoo biome in South Africa, their inherent properties and possible impacts on plant communities, Table S3: Families containing diazotrophic bacteria (nifH gene) in the rhizosphere soils examined and their possible roles in plant growth promotion of legume species in the Succulent Karoo biome in South Africa, Table S4: Analysis of variance (ANOVA) and Permutations analysis of variance (PERMANOVA) to test the effects of factors 'legume species', 'sites' and their interactions on the respective alpha and beta diversity matrices, using the 16S rRNA and nifH gene barcodes.