The selected sampling sites in Lichtenburg (25°59′40.8″S, 26°31′46.6″E) and Randfontein (26°11′51.3″S, 27°33′18.6″E) are maize growing areas in South-Africa which have been under intensive cultivation continuously for more than 10 years. They are both characterized with shrubs and grasses being a semi-arid region, the maize cultivar WE 3128 was planted on both sites. The average annual temperature and rainfall is around 16.9°C and 601 mm respectively in Lichtenburg while in Randfontein, it is around 15.3°C and 742 mm respectively. The type of soil found on the two farm sites can be classified as the Luvisol type based on the World reference base for soil resources classification (Anjos et al., 2015).

After careful excavation of the maize plants, the roots were shaken by hand to remove soil loosely adhering to the roots, while the soil tightly attached to the roots was collected as rhizosphere soil. Bulk soils were collected from 0–15 cm depth, 10 m away from the cultivated area through a soil auger (Castañeda and Barbosa 2017; Enagbonma et al., 2019). Samples were collected in triplicates for both rhizosphere and bulk soils and were transported on ice (so as to prevent dehydration and degradation) to the laboratory immediately after collection for analysis. Soil samples were passed through a 2 mm sieve to remove unwanted particles and stored at −20°C until use for downstream applications.

The soil properties were measured using 20 g of air-dried soil from each sample, pH meter (Eutech^TM pH 150/Ph450 series) was used in measuring the soil pH in 1:2.5 soil water ratio (Springer 2014). A Gravimetric technique was used to determine the soil moisture as described by Shukla et al. (2014), and Morche’s (2008) method was used to determine sulfur concentration following extraction with 0.1 M hydrochloric acid (HCl). Total carbon and soil organic carbon were determined using the Santi et al. (2006) method, and the Walkey-Black procedures (Walkley and Black 1934) respectively.

The concentration of potassium in samples was measured following extraction with 1M ammonium acetate at pH of 7.0 (Jackson 1973). Bray and Kurtz’s protocol (Bray and Kurtz 1945) was used in measuring phosphorus concentration in samples, the described methods of Mussa et al. (2009) was used in determining soil nitrate, Potassium chloride (KCl) extraction method was used to analyze soil nitrate and ammonium (Reis et al., 2017), and the Loss-on-Ignition method was used to determine soil organic matter content following the procedure of Schulte and Hopkins (1996).

Following manufacturer’s protocol, QIAGEN DNeasy PowerSoil Isolation Kit was used for the extraction of total DNA using 0.5 g of each soil sample. Following standard protocol, the concentration of DNA in samples were first examined using the Life Technologies Qubit^® dsDNA HS Assay Kit. Nextera DNA Flex library preparation kit (Illumina) was used to prepare the libraries using about 50 ng DNA. Simultaneous fragmentation and addition of adapter sequences were done on the extracted DNA. The final concentrations of the prepared libraries were measured using the Qubit^® dsDNA HS Assay Kit (Life Technologies), and the average library size was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies). The libraries were pooled and diluted to 0.6 nM and sequenced paired end for 300 cycles using the HiSeq system (Illumina). Sequencing was done at the Molecular Research Laboratory (MR DNA) LP, Shallowater Texas, USA.

Unassembled DNA sequences were annotated with the Metagenomics Rapid Annotation (MG-RAST) pipeline version 3.3 following the prescribed lead of Meyer et al. (2008). Raw sequences were uploaded to the MG-RAST online server at https://www.mg-rast.org (Meyer et al., 2008). Following quality control, sequence annotation was done using BLAT algorithm against the M5NR database (Wilke et al., 2012) which provides a non-redundant incorporation of numerous databases such as KEGG (Kyoto Encyclopedia of Genes and Genomes), UniProtKB (UniProt Knowledgebase), and GenBank. Microbial taxonomic and gene profiling were performed using the M5NR and SEED level Subsystems respectively. BlastX with an e-value cutoff of E < 1 × 10⁻⁵, the minimum identity of 60% and a maximum alignment length of 15 base pairs was used to detect hits. Data normalization was applied on the MG-RAST pipeline so as to suppress the experimental noise effect. Sequences that failed annotation were left out in the analysis procedures. Carbon fixation, sulfur cycling, phosphorus utilization, nitrogen fixation, and other plant beneficial genes were categorized from the entire metagenomes obtained.

The average values of the relative abundance of genes obtained from the three replicate samples from each sampling site were used for data analysis. Following normalization of the dataset, the rarefaction curve was plotted using the analysis tool in the MG-RAST server. The differences in the soil properties between the soil samples were analyzed using one-way analysis of variance with Tukey’s pairwise comparison test for significance level (p < 0.05) (Gunsalus et al., 2016). Assumptions (normality and homogeneity of variance) were determined and met. PAST (Paleontological Statistics) statistical software version 3.2 was used to evaluate the alpha diversity indices (Evenness, Simpson and Shannon indices) and Kruskal-Wallis test was used to determine the significance level (Hammer et al., 2001). Beta diversity was analyzed using one-way analysis of Similarities (ANOSIM) via 999 permutations and further illustrated by principal coordinate analysis (PCoA) based on Euclidean distance-matrix. (Lee et al., 2016). Data representation on the Principal coordinate analysis (PCoA) plot was done using the Canoco 5.12 software. Circos online software (www.circos.ca) was used in generating the chart for the gene abundance (Krzywinski et al., 2009). Representation and visualization of microbial taxonomy at the family level was carried out with a graph made using the Microsoft Excel software, version 2013. Quality sequences are deposited and available from NCBI SRA dataset under the accession number PRJNA678469 and PRJNA678475.

As indicated in Table 1, the pH of R1 (Rhizosphere soil samples from Randfontein) and R2 (bulk soil samples from Randfontein) are 6.54 and 5.89 respectively, which differ significantly from that of F1 (rhizosphere samples from Lichtenburg) and F2 (bulk soil samples from Lichtenburg) which are 6.76 and 6.70 respectively. Sulphur and total carbon concentration did not differ significantly between R1 and R2, while a significant difference was seen in these two between F1 and F2. Potassium (K), N-NO₃ and N-NH₄ were significantly higher in the rhizosphere soils than their corresponding bulk soils. The concentration of organic carbon varied significantly between F1 and F2, while no significant difference was seen between R1 and R2. No significant difference was seen in the moisture content across the samples, however in the concentration of organic matter there was a significant difference between F1 and F2, while R1 and R2 did not differ significantly. There was no significant variation in the phosphorus content of F1 and F2, while in R1 and R2, these differ significantly.

The majority of the sequence read reached saturation, signifying that adequate data was obtained from the samples (Figure 1). The average number of sequences uploaded was 18, 942,494, 18,442,769, 15,666,716, and 15,045,869 for the rhizosphere soil samples from Lichtenburg (F1), bulk soil samples from Lichtenburg (F2), rhizosphere samples collected in Randfontein (R1), and bulk soil samples collected in Randfontein (R2) respectively. After quality control (QC), quality sequences that were retained in samples were 17,309,422, 16,779,302, 14,404,078, and 13,867,146 in F1, F2, R1, and R2 respectively with GC (Guanine-Cytosine) content of 65% in F1, F2, R1, and 66% in R2. Sequences with known protein functions were 5,391,241, 5,549,864, 4,763,335, and 4,954,381, for F1, F2, R1, and R2 respectively.

The dominant microbial class observed in the samples include Actinobacteria, Alphaproteobacteria, Betaproteobacteria, Deltaproteobacteria, Gammaproteobacteria, Planctomycetacia, Gemmatimonadetes, Solibacteres, Verrucomicrobiae, and Bacilli (Figure 2A). The dominant microbial families in the samples include Conexibacteraceae, Nitrospiraceae, Rubrobacteraceae, Nocardioidaceae, Planctomycetaceae, Nocardiaceae, Solibacteraceae, Rhizobiaceae, Burkholderiaceae, Pseudomonadaceae, Cytophagaceae, and Rhodobacteraceae (Figure 2B). No significant variation was seen in the relative abundance of the microbes in the samples (p > 0.05).

The metagenomic sequences obtained from the samples consist of genes involved in nutrient cycling (carbon fixation and degradation, nitrogen cycling, phosphorus utilization, and sulfur cycling) and other plant beneficial genes such as those involved in siderophore production, iron acquisition, phosphate solubilization, heat and cold shock (Figures 3, 4A,B). These genes aid plant growth and development and the presence of these in the maize rhizosphere indicated their contribution to maize growth.

The alpha diversity analysis (diversity within the habitats) of the microbial genes showed no significant difference (p > 0.05) in the rhizosphere and bulk soil (Figure 5). However, the beta diversity analysis showed significant variations between the habitats (p < 0.05), this was further illustrated by the principal coordinate analysis which shows distinct separation between the samples on the PCoA plot (Figure 6).

Among the genes detected in this category are those which code for Ribulose-1,5-biphosphate carboxylase/oxygenase (RuBPco/RUBISCO, which is one of the enzymes involved in Calvin cycle), Propionyl-CoA carboxylase (Pcc, which catalyzes the carboxylation reaction of Propionyl CoA), and carbon monoxide dehydrogenase (CODH, involved in the reductive acetyl-CoA pathway in carbon cycle which allows organisms to make use of carbon monoxide as a source of energy and carbon (IV) oxide as a source of carbon (Li et al., 2014d)).

Genes associated with Ribulose-1,5-bisphosphate carboxylase/oxygenase (cbbX, cbbRm), carbon monoxide dehydrogenase (coxF, coxG, coxM, cooC, coxE) and Propionyl-CoA carboxylase (pccB, pccX) were seen to be more abundant in the rhizosphere soil samples (F1 and R1) than the bulk soils (Figures 3, 4A). In the rhizosphere samples, it however predominates F1 than R1 signifying a higher capacity for carbon fixation in this sample (Figures 3,4A) (Li et al., 2014c). Mechanisms of carbon fixation have been identified in the past in different species of Rhodobacteraceae (Rhodobacter sphaeroides, Rhodobacteraceae bacterium KLH11), Thermoproteaceae (Thermoproteus tenax), Thermococcaceae (Thermococcus onnurineus), Halobacteriaceae (Haloquadratum walsbyi), Pseudomonadaceae, and Vibrionaceae, (Dijkhuizen and Harder 1984; Berg et al., 2010; Li et al., 2014c; Pujalte et al., 2014), therefore the prevalence of these organisms in F1 (Figure 2B) justifies the abundance of these genes observed in the sample.

As shown in Figures 3, 4A, Genes associated with inorganic polyphosphate degradation (ppx, ppgk), polyphosphate biosynthesis (pgk 1, ppk, thiD), dephospho-coenzyme A kinase (dpck), which catalyzes the ATP-dependent phosphorylation of dephospho-CoA in the biosynthesis of coenzyme A (CoA) were well represented in all the samples. The prevalence of these genes in F1 can be linked with the prevalence of Thermococcaceae, Rhizobiaceaea in this than R1 (Schut et al., 2014; Bechtaoui et al., 2019; Shimosaka et al., 2019). Also, the prevalence of ppx and ppgk in F1 suggests an intense polyphosphate transformation and more abundance of inorganic phosphorus (Li et al., 2014c) in the sample. Bacillariaceae was seen to predominate F1 rather than R1, some species of this bacteria family like the Bacillus cereus are known to be active in phosphorus metabolism. Bacillus cereus secretes different phospholipases which helps in phosphate retrieval metabolisms, it has been noted also for its ability to produce polyphosphate through the biodegradation of the herbicide glyphosate (Guddal et al., 1989; Acosta-Cortés et al., 2019).

The occurrence of siRA, suoX, papsS1, and dsrAB in the studied environments signifies the presence of microbes capable of metabolizing sulfur. siRA, genes were profound in F1 than other samples, thus indicating a high rate of sulfur transformation in F1 than R1 and other samples (Li et al., 2014c). Contrarily, genes involved with Sulfite oxidase and Sulfite reductase for oxidation and reduction of sulfur such as the suoX, papsS1 and dsrAB were more dominant in the R1 than F1 samples (Figure 4A).

Genes involved in nitrogen fixation like the nitrite and nitrate reductase genes (nifN, nifA, nirD, napH, napC, napF) and those responsible for nitrite and glutamase transport (nirC, gltk and gltj), nitric and nitrous oxide (norQ and nosX), as well as the nirM, nirH, nirL were also revealed in the metagenomes (Figures 3, 4A). This can be affiliated with the occurrence of Rhodobacteraceae, Rhizobiaceae, Nitrospiraceae, Methanocellaceae, and Burkholderiaceae in the samples (Figure 2B). Species of these bacteria families such as Burkholderia vietnamiensis, Burkholderia cepacia, Agrobacterium fabrum (strain C58), Bradyrhizobium diazoefficiens (strain JCM 10833) Azotobacter vinelandii, and Klebsiella pneumoniae are associated with nitrogen fixation. The richness of F1 in nifH and nodA genes could be linked to the abundance of Rhizobiaceae, Burkholderiaceae in this sample as some species of this bacteria family like the Rhizobium leguminosarum, Paraburkholderia xenovorans, Cupriavidus basilensis, Cupriavidus necator are known for fixing nitrogen (Bournaud et al., 2017; da Silva et al., 2012; Estrada-de los Santos et al., 2011).

Furthermore, the metagenomes contain genes pqqABCDEF and gcd which play active roles in the biosynthesis of pyrroloquinoline quinone (PQQ) and transformation of insoluble phosphate to dissolved phosphate (soil phosphate solubilization), thus making it available for plant use (Otieno et al., 2015). According to Li L. et al. (2014) the above mentioned genes perform important roles in plant growth and biocontrol of crown gall disease caused by Agrobacterium tumefaciens.

Phosphate limitation in the soil triggers the autophosphorylation of histidine residue on sensory kinase, thus the presence of phoH and phoU complex which is involved in response to phosphate starvation and transport helps to control the metabolism and uptake of phosphate in response to phosphate limitation to active target sites (Akinola et al., 2020a). Occurrence of phoHU genes in all the metagenomes shows the ability of the resident soil microbes in enhancing limited soil phosphate. This is in line with the study of Myo et al. (2019) which showed that the effect of Streptomyces lividans TK24 in reshaping phosphate limited soil was as a result of the pho complex possessed by the organism. The occurrence of high affinity phosphate transporter pstA gene was also discovered in the metagenome.

Another hormone which is significant and essential for proper growth in plants is the Indoleacetic acid (IAA) (Myo et al., 2019). The ppdC and ipdC genes which code for phenylpyruvate decarboxylase and Indole-3-pyruvate decarboxylase which are the enzymes that produce Indole acetic acid from tryptophan (Jijón-Moreno et al., 2015) were found in the soil samples. Indole acetic acid assists in modulating numerous physiological processes like cell division, geotropisms, phototropism, root initiation, apical dominance and growth rate in plants (Singh and Dubey 2018; Carrión et al., 2019). The occurrence of the ipdC genes in the metagenomes can be associated with Enterobacter cloacae, a species of the Enterobacteriaceae family which occurred in the samples.

In the relationship between plant and bacteria, ACC deaminase reduces plant ethylenes notable for hindering the nodulation process (Kang et al., 2019). Genes that encode for ACC deaminase enzymes such as rimM were found more in the rhizosphere than bulk soil samples, this agrees with previous findings in which abundance of this gene was discovered in the rhizosphere than the bulk soil (Li et al., 2014c; Franco 2015).

Furthermore, stress alleviating genes such as the cspA, cspB, cspC, cspD, cspE, cspF, cspI, responsible for reducing cold-shock and groES, grpE, yegD, yciM, dnaJ which mitigates heat-shock were also discovered with more abundance in the rhizosphere soils. The higher abundance of heat shock genes in R1 shows that its microbiota demonstrates higher resistance to heat and stress caused as a result of drought. Drought is a major abiotic stress for plant especially in dry/arid regions like South Africa, it reduces plants potentials for growth and productivity (Nicolaes et al., 2014; Aydin et al., 2019). Moreover, the higher abundance of cold shock genes observed in the F1 than R1 sample shows the higher ability of its inherent microorganisms in cushioning plants against the harsh effect of downshift in temperature, thus aiding plant cells to adapt (Keto-Timonen et al., 2016). The occurrence of these stress reducing genes in our metagenomes shows the potentials of the sampling sites in reducing/alleviating the threat posed by heat and cold stress on plant/crop.

The fhuA gene which is a Ferrichrome-iron receptor in Escherichia coli and also a family of outer membrane proteins was seen in the samples. This gene together with the TonB-ExbB-ExbD which are the energy transducing complex of the cytoplasmic membrane are important in mediating the active transport of vitamin B_12, siderophores and heme across the outer membrane of bacteria (especially gram-negative bacteria) (Zeng et al., 2013). In addition, the fhuA gene acts as a receptor for numerous bacteriophages, antibiotics and toxins which are harmful to bacterial cells. Rhizobium leguminosarum which belongs to the Alphaproteobacteria class are known to possess this gene (Carson et al., 2000; Akinola et al., 2020b), therefore the greater abundance of this class of bacteria in F1 than R1 and other samples, can be associated with the larger abundance of this gene in F1 than other samples.

The sbnF and sbnC genes are important in staphyloferrin B biosynthesis and are both associated with Staphylococcus aureus. The sbnF catalyzes the condensation of L-2,3-diaminopropionic acid (L-Dap) and citryl-diaminoethane to form L-2,3-diaminopropionyl-citryl-diaminoethane while the sbnC catalyzes the condensation of L-2,3-diaminopropionyl-citryl-diaminoethane and 2-oxoglutarate to form staphyloferrin B. The prevalence of these genes in F1 than R1 can be linked to the prevalence of Bacilli in this sample, since Staphylococcus aureus belongs to this bacteria class (Barka et al., 2016)

Environmental variables such as pH and the concentration of phosphorus, sulfur, total carbon and potassium, explains the distribution of plant beneficial microbes with plant growth promoting/beneficial genes in soil samples (Finn et al., 2020). As seen in this study, there were differences in the pH values of F1, F2, R1, and R2 (Table 1), sulfur concentrations varied significantly between F1 and F2, phosphorus differ significantly between R1 and R2, while between F1 and F2, significant difference was seen in potassium concentration and also between R1 and R2. The concentration of N-NO₃ in the rhizosphere soils (8.52 and 17.09 in R1 and F1 respectively) showed that the biological processes which occur in these region are capable of improving growth in maize (Tale and Ingole 2015).

The substantial concentration of phosphorus and potassium in the rhizosphere samples indicated that immense physiological processes that are vital in the metabolism and biosynthesis of vital nutrients necessary for optimum growth in plants occur in their environments (Peralta et al., 2013; Akinola et al., 2020a).

The dissimilarities observed in the rhizosphere samples examined in this study could be linked to the peculiarities of the sampling location. Lichtenburg, where the F1 samples were collected is an agricultural area, while Randfontein, where R1 samples were collected comprises engineering industries and a gold mine. Thus, waste and residues such as heavy metals could have infiltrated into the soil, therefore, altering the physicochemical properties, community structure, and diversity of the soil microbes. Uranium extracted from gold ores is known for negatively impacting the soil microbial community structure and affecting genes such as those responsible for carbon and nitrogen cycling as seen in this study (Sutcliffe et al., 2017). Chromium is another metal associated with gold mines and has been known to alter the soil’s bacterial diversity (He et al., 2016). These heavy metals affect the soil microorganisms and interfere with soil microbial metabolisms.

The knowledge of soil microbial metabolisms is important in understanding interactions which enhance the ecosystem and improve plant growth and production (Jiménez et al., 2016; Igiehon and Babalola 2017; Enagbonma and Babalola 2020). This study revealed an ecologically diverse microbiome with nutrient cycling/plant growth-promoting traits. In the metagenomic reads, microbes such as Rhizobiaceae, Burkholderiaceae, Frankiaceae, and Pseudomonadaceae, Vibrionaceae were discovered. Certain strains of these bacteria families harbor genes which carry out nitrogen fixation, sulfur cycling, production of phytohormones (such as cytokinins, auxins, and gibberellins which help plants in nutrient uptake), inhibition of harmful pathogens and other growth and development processes in plants (Li et al., 2014d; Omotayo et al., 2019).

The alpha diversity indices showed that the microbial genes did not differ significantly within the habitats as revealed by Kruskal-Wallis test (p = 0.99). However, significant variation (p = 0.01) was seen in the beta diversity analysis which was determined using analysis of similarities (ANOSIM). This was further illustrated on the principal coordinate plot (Figure 6) on which there was a distinct separation between the soil samples. The strength of the separation (R) is 0.58, which depicts a clear difference between the samples. A notable distance was seen between the two rhizosphere soils (F1 and R1) which were also clearly distinguished from their corresponding bulk soils. The distance between each sample as revealed by the PCoA plot with a cumulative variation of 97.35% between the two axes (Figure 6) showed that the metagenomes are distinctly different in genetic structure.

The metadata obtained from this study showed and supports the notion that maize rhizosphere is a hotspot of genes accountable for converting labile and organic carbon, sulfur, nitrogen, and phosphorus compounds (Li et al., 2014c; Omotayo and Babalola 2020; Omotayo et al., 2021). This can be linked to the release of root exudates into the rhizosphere, as these rhizodeposits have a major influence on the rhizosphere microbiome, thus impacting its microbial diversity and functional profile, consequently affecting the microbial genes (Hu et al., 2018). Root exudates also possess the capacity to enhance latent soil organisms and activate their potentials. A wide range of compounds such as amino acids, lipids, and enzymes are present in the root exudates which influences the availability and wealth of potent genes in this region (Hu et al., 2018).

More potential to catabolize complex aromatic compounds was shown in the rhizosphere soils, especially R1 due to the abundance of Pseudomonadaceae in this sample (Jiménez et al., 2010; Li et al., 2014c). Strains of bacteria belonging to this family such as the P. putida KT2440, P. putida CSV86 are known for its ability to colonize the rhizosphere of various plants and degrade aromatic compounds (Basu et al., 2006; Jiménez et al., 2010).