Metagenomic Study of the Community Structure and Functional Potentials in Maize Rhizosphere Microbiome: Elucidation of Mechanisms behind the Improvement in Plants under Normal and Stress Conditions

: The community of microbes in the rhizosphere region is diverse and contributes signif-icantly to plant growth and crop production. Being an important staple and economic crop, the maize rhizosphere microbiota has been studied in the past using culture-dependent techniques. However, these limited culturing methods often do not help in understanding the complex community of microbes in the rhizosphere. Moreover, the vital biogeochemical processes carried out by these organisms are yet to be fully characterized. Herein, shotgun metagenomics, which enables the holistic study of several microbial environments, was employed to examine the community structure and functional potentials of microbes in the maize rhizosphere and to assess the inﬂuence of environmental variables on these. The dominant microbial phyla found in the soil environments include Actinobacteria, Microsporidia, Bacteroidetes, Thaumarchaeota, Proteobacteria and Firmicutes. Carbohydrate metabolism, protein metabolism and stress metabolism constitute the major functional categories in the environments. The beta diversity analysis indicated signiﬁcant differences ( p = 0.01) in the community structure and functional categories across the samples. A correlation was seen between the physical and chemical properties of the soil, and the structural and functional diversities. The canonical correspondence analysis carried out showed that phosphorus, N-NO 3 , potassium and organic matter were the soil properties that best inﬂuenced the structural and functional diversities of the soil microbes. It can be inferred from this study that the maize rhizosphere is a hotspot for microorganisms of agricultural and biotechnological importance which can be used as bioinoculants for sustainable agriculture. This study provided insight into the structure and functional potentials of the community of microbes in the maize rhizosphere. The understanding of the microbes which help in the functioning of the rhizosphere ecosystem would assist in the development of strategies which can be used to increase agricultural production through manipulation of these soil microorganisms. This study also revealed that the maize rhizosphere is endowed with microorganisms such as Myxococcales, Burkholderia, Haloferax and Aquiﬁcales which are important in plant growth, helping plants to survive in low-nutrient soil, and in the phytoremediation of heavy metals, and production of iron-chelating compounds and phy-tohormones. These microorganisms, when isolated, can serve as important bioinoculants which can be used for increased crop production and sustainable agriculture. Moreover, this study shows the importance of nitrogen, phosphorus, potassium, organic matter and N-NO 3 , as these are the major environmental variables that shaped the microbial community structure and also inﬂuenced the metabolic activities of microbes in the maize rhizosphere microbiome. Thus, the soil environment cannot afford to be lacking the accurate proportion of these soil nutrients if there would be optimum production.


Introduction
The presence of various microorganisms in the rhizosphere region makes it an important interface in the exchange of nutrients and resources between plants and the soil [1][2][3]. The impacts of microbes in this zone are manifold, ranging from the decomposition of organic matter to the disintegration of metabolic by-products which enhances the availability of nutrients and essential elements [4]. Several metabolic processes such as protein metabolism, sulfur cycling and phosphorus and potassium metabolism, which occur in the rhizosphere, contribute to the overall wellbeing of plants [5].
The secretion of organic compounds through the root exudates into the rhizosphere causes complex metabolic interactions and an interplay among the organisms that could be beneficial, neutral or harmful to the plants [6] (microbes present in the rhizosphere region form a symbiotic association with plants, and they are connected together by diverse pathways [7]). Plants release some of their photosynthetically fixed carbon as carboxylic acids, sugars, amino acids and various secondary metabolites into the rhizosphere through the root exudates [2]; these are then used as energy sources for the soil microorganisms which in turn provide growth, development and sustainability for the plant [8].
Maize (Zea mays) is a commonly cultivated staple crop, with billions of tons produced annually globally. The wide availability and consumption of this staple crop, its substantial contribution to the economy and its relevance as a raw material in the production of various commodities make it indispensable [9]. To match the high demand of this crop and the growing global population, there is a need to increase its production. Yet, to achieve sustainable maize farming, it is essential to understand the taxonomic/structural and functional potentials of maize rhizosphere microbes, especially those that aid the functioning of the soil ecosystem, plant growth and development, as this will help to develop effective techniques to increase agricultural production through the manipulation of soil microorganisms.
Efforts have been made in the past to characterize the rhizosphere microbiome of certain crops such as barley, rice, soybean and wheat through culture-dependent methods [10,11], but these do not provide a full understanding of the rhizosphere microbiome and its functional traits. The advent of next-generation sequencing/high-throughput sequencing techniques (shotgun metagenomics) has provided an advantage over the use of traditional culturing procedures in understanding the metabolic activities/processes occurring in the rhizosphere. These techniques allow the use of bioinformatics tools to analyze the soil microbial taxonomy and function [12]. Through shotgun metagenomics, more direct assessments and broader insights into the activities and functional attributes of resident rhizosphere organisms can be achieved than through optimizing culture-dependent procedures. Moreover, studies which profile the microbial communities (bacteria, fungi and archaea at different levels) in maize rhizosphere soil (especially of semi-arid regions), their functional potentials and the effect of environmental variables on these remain underexplored in South Africa. Thus, using metagenomics, this study seeks to:

•
Determine the structure and function of microbial communities in maize rhizosphere and bulk soils (of two major maize agroecosystems in semi-arid South Africa); • Determine the effect of environmental variables on the microbial communities.
We assumed that the maize rhizosphere soil would be endowed with beneficial microbes and functional activities which are of agricultural importance, and that the structure and functional potentials of these microbes would be influenced by environmental variables. This study provides insight into the structure and functional capabilities of microbial communities in maize rhizosphere soils.

Site Description and Soil Sampling
The selected sites for this study were those for maize production, being located in Lichtenburg (25 • Figure 1). The sampling locations are renowned regions for maize cultivation and maximum production of maize in South Africa. Both locations, which are semi-arid regions, are characterized by shrubs and grasses, the average annual temperatures being 16.9 • C in Lichtenburg, and 15.3 • C in Randfontein. Rainfall occurs from September to May in both locations, the annual amounts being approximately 601 mm and 742 mm in Lichtenburg (F) and Randfontein (R), respectively. Based on the world reference base for soil resources classification, the type of soil found on the two farm sites can be classified as the Luvisol type [13]. At the time of sampling, the maize cultivar planted on both farm sites was WE 3128, and the maize plants were at the flowering stage. In the two fields, conservation tillage (specifically no till) has been the adopted method, and artificial irrigation systems are used in support of rainfall.
Sampling was conducted in March 2019, which is the summer season of the Republic of South Africa. In order to collect the rhizosphere soils, maize plants were excavated carefully, and the soil loosely attached to the root was removed, while that bound to the root was collected. Bulk soils were also collected from 0 to 15 cm depth at a distance of 10 m away from the cultivated area (this depth was chosen because major microbial activities take place in this region [14,15]). Samples were collected in triplicates for both rhizosphere and bulk soils and transported on ice to the laboratory within 24 h. These were thereafter pre-processed (passed through a 2 mm sieve to remove pebbles and non-soil particles) and stored at −20 • C until use for downstream applications.

Analysis of Soil Properties, DNA Extraction and Metagenomic Sequencing
Using various methods, the soil properties (physical and chemical) were measured using 20 g of soil from each sample. A pH meter was used in measuring the hydrogen ion concentration in a 1:2.5 soil/water ratio [16]. The gravimetric technique was used in determining the soil moisture, as described by Shukla et al. [17]; the sulfur concentration in samples was measured following extraction with 0.1 hydrochloric acid (HCl) using the Morche [18] method; and total carbon and soil organic carbon were determined using the Santi et al. [19] method, and the Walkey-Black procedures [20], respectively. The concentration of potassium in samples was measured following extraction with 1M ammonium acetate at a pH of 7.0 [21], and phosphorus was measured using Bray and Kurtz's protocol [22]. The potassium chloride (KCl) extraction method was used to analyze soil nitrate and ammonium, absorbance was measured at 220 nm and 260 nm [23] and the loss-on-ignition method was used to determine the soil organic matter content following the procedure of Schulte and Hopkins [24].

DNA Extraction and Sequencing
Using the QIAGEN DNeasy PowerSoil Isolation Kit, whole microbial DNA was extracted from 5 g of each sampled soil following the manufacturer's instructions. After assessing the concentration and quality of the extracted DNA with the use of a Nano Drop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA), it was subjected to wholemetagenome shotgun metagenomic sequencing on the Illumina Hiseq platform, which was carried out at the Molecular Research Laboratory, Shallowater, Texas, USA. With shotgun metagenomics on the Illumina Hiseq, large datasets with whole genome information can be obtained (since our target is to obtain the structural and functional information of microbial communities), and the high accuracy and sensitivity of this method are a major advantage. This method stands out from other next-generation sequencing method such as Roche 454 because it has the biggest output and lowest reagent cost, and it provides unbiased genome coverage and high-quality pairwise alignment [25,26]. Preparation of the libraries was conducted using the Nextera DNA Flex library preparation kit (Illumina) following the standard protocol, and the Qubit ® dsDNA HS Assay Kit (Life Technologies) was optimized to evaluate the initial concentration of DNA, after which a total of 50 ng DNA was used to prepare the libraries. The samples were passed through the simultaneous fragmentation and addition of adapter sequences, which were used during a limited-cycle PCR when unique indices were added to the sample. The final concentration of the libraries was 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 then pooled, diluted (to 0.6 nM) and paired-end sequenced for 300 cycles using the NovaSeq system (Illumina).

Bioinformatics Analysis of Sequences
Using default settings, unassembled DNA sequences were annotated on the Metagenomics Rapid Annotation pipeline (MG-RAST) version 3.3, following the method of Meyer et al. [27]. The sequences were uploaded to the MG-RAST online server at https://www.mg-rast.org (MG-RAST online server at https://www.mg-rast.org, accessed on 10 June 2020), quality control (QC) was conducted to remove adapter and artificial sequences (that were formed by sequencing artifacts) and the ambiguous bases (that is, removal of sequences with >5 ambiguous base pairs with 15 phred score cutoff) were filtered. Sequences with a length of >2 standard deviation from the mean were removed, and reads of low quality were discarded. Following quality control (QC), quality sequences were annotated using the BLAST-like alignment tool (BLAT algorithm) [28] against the M5NR database, which provides non-redundant incorporation of various databases [29]. Taxonomic classifications were conducted using the SEED Subsystem, and assigning the functional categories was performed using the SEED Subsystem level 1 and level 2 databases (the subsystems include anabolic and catabolic processes at the highest level of organization, while at the lowest level, they include specific pathways). These were conducted with an e-value of 1 × 10 −5 , a minimum identity of 60% and a maximum alignment length of 15 base pairs. The data normalization tool was applied on MG-RAST to subdue the influence of experimental noise. Following an independent analysis of the sequences in MG-RAST (the triplicate samples for each sampling site), the average values of their relative abundance were used for statistical analysis. The whole genome raw sequences of the soil samples were deposited and are available from the NCBI SRA dataset under the accession numbers PRJNA678469 and PRJNA678475.

Statistical Analyses
The differences in the soil properties between the soil samples were analyzed using one-way analysis of variance with Tukey's pairwise comparison. A taxonomic abun-dance heat map was generated through Heatmapper and the Morpheus online tool at https://software.broadinstitute.org/morpheus/ (accessed on 21 August 2020) [30,31]. Diversity indices (Pielou evenness, Simpson and Shannon-Wiener indices) were evaluated and compared across the samples using the Kruskal-Wallis test on PAST (Paleontological Statistics) statistical software version 3.2 [32]. Based on the Euclidean distance matrix, beta diversity was analyzed and illustrated using principal coordinate analysis (PCoA). Oneway analysis of similarities (ANOSIM) via 999 permutations was used to test for differences in the functional categories of the sampling sites [33]. Principal component analysis (PCA) was performed to show the distribution of the functional categories across the samples. Canonical correspondence analysis (CCA) was conducted to find the experimental variables that best explained the functions, and forward selection of environmental variables and a test of significance using the Monte Carlo test with 999 random permutations were performed. All examined soil properties [34] (Table 1) were included as explanatory variables in the canonical correspondence analysis. Data representation in the PCoA, PCA and CCA was conducted using the Canoco 5.12 software (default settings were used in all software packages) [30,31].

Analysis of Soil Properties (Physical and Chemical Analysis of Soil)
As it is indicated in Table 1, variations were observed in the pH of the soil samples, which were acidic in nature. The pH and P of R1 and R2 differ significantly from those of F1 and F2. Sulfur and total carbon concentrations did not differ significantly between the samples. K and N-NH 4 were significantly higher in the rhizosphere soils than their corresponding bulk soils. The concentration of organic carbon varied significantly between R1 and R2, while no significant difference was seen between F1 and F2. Although the moisture content in the samples varied, it did not differ significantly, and no significant difference was seen in the concentration of organic matter in samples.

Metagenomic Sequence Information
For the rhizosphere soil samples from Lichtenburg (F1), the average number of sequences uploaded was 18,942,494, while for the bulk soil samples (F2), it was 18,442,769. The average number of uploaded sequences from the rhizosphere samples collected in Randfontein (R1) was 15,666,716, while for the bulk soil samples (R2), it was 15,045,869 on average. Following quality control, the quality sequences which were retained in F1 totaled 17,309,422, with a GC content of 65%; F2 had 16,779,302, with a GC content of 65%; and 14,404,078 sequences were retained in R1, with a GC content of 65%, while 13,867,146 sequences, with a GC content of 66%, were found in R2. Sequences with known protein functions totaled 5,391,241 (F1), 5,549,864 (F2), 4,763,335 (R1) and 4,954,381 (R2) (Table S1).

Distribution of Major Microbial Phyla Across the Sites
Metagenomic analysis using the SEED Subsystem database revealed twenty-eight (28) phyla in all samples ( Figure 2, Table S2). The major observed phyla in samples include Proteobacteria, Actinobacteria, Gemmatimonadetes, Bacteroidetes, Firmicutes, Chloroflexi, Planctomycetes, Crenarchaeota and Basidiomycota, as shown in Figure 2. The majority of the microbial phyla such as Proteobacteria, Firmicutes, Chloroflexi, Planctomycetes, Verrucomicrobia, Chytridiomycota and Cyanobacteria, among others, were found to occur more in F1, while Actinobacteria, Basidiomycota, Thaumarchaeota and Microsporidia occurred more in R1 than in F1. Cyanobacteria, Actinobacteria, Firmicutes, Chloroflexi and Planctomycetes were found to be more predominant in F2 than in R2, while Verrucomicrobia and Gemmatimonadetes occurred more in R2 than in F2. F1 has a higher abundance of the microbial phyla, and a considerable difference was seen in the distribution of these microbes across the examined sites; however, within the habitat (each environment), the Kruskal-Wallis test conducted showed no significant difference (p > 0.05).
A heat map, which is a data visualization technique, was used to represent the relative abundance of taxa in this study from the phylum to the genus level. Being color-shaded matrices where the color of each cell portrays its value, the data are in intervals, with each interval having a color from a divergent color scheme (red-green). The heat map presents the abundance of the studied taxa from the lowest to the highest, as illustrated in the legend. Clustering in heat maps ( Figure 2) was used in recognizing groups of correlated samples [35], and these were revealed as structures which are block-like along the diagonals (dendrogram). Analysis carried out on Heatmapper (http://www.heatmapper. ca/expression/, accessed on 21 August 2020) was conducted using default parameters, and column clustering and a dendrogram were applied to show the similarities in the sites with regard to the microbial community structure using the average linkage and Euclidean distances.
The dendrograms, which help to combine heat maps with hierarchical clustering, arranged the sites based on the similarity/distance between them. The clustering seen in Figure 2 between the rhizosphere soils (F1 and R1) showed that the microbial communities in these samples share some similarities, hence the close distance between them. This could be associated with the root exudates released by the same type of maize cultivar which was planted on the two sampling sites. Furthermore, the clustering observed between the two rhizosphere soils and the bulk soils (F2 and R2) could be linked to the fact that rhizosphere soil organisms are recruited from the bulk soil microorganisms.
Principal component analysis (PCA) was carried out to determine the distribution of the microbial phyla between the sampling sites ( Figure 3). Axes 1 and 2 of the PCA plot explained 77.81% and 9.09% of the variance, respectively. Thus, PCA1 and PCA2 contribute 86.9% of the total variation in the dataset. The vector length of the PCA shows the phylum which dominates each of the habitats. From Figure 3, the vectors show the microbial communities which were distributed and predominate each sample. Variables which correlate with each other are shown by vectors pointing in the same direction, and the samples that are close to each other on the PCA show they share similarities. It can be seen from Figure 3 that most of the variables (microbes) contributed to PCA 2, with the higher values in the variables moving the samples to the left of this plot. The angle between the variable vectors which is small shows that the variables are highly associated.

Diversity Estimation for Microbial Communities in the Rhizosphere and Bulk Soils
Alpha diversity indices (Shannon and evenness) at the phylum, class, family, order and genus level were within the same range across the samples (Table 2), with no significant difference (Kruskal-Wallis p > 0.05). However, a significant difference was observed in the beta diversity analysis conducted using one-way ANOSIM (p = 0.01). The R value of 0.58 obtained from the beta diversity analysis implied a clear difference/dissimilarities in the samples, and this is shown in the PCoA plot ( Figure 5). From the conducted principal coordinate analysis, which was used to assess/compare the samples, no definite clustering was seen in the examined environments. Samples F1a, b and c were distinct and located away from R1a, b and c, and F2a and b were clearly different from R2a and b, implying that the microbial community structures of these samples varied from each other. However, samples F2c and R2c were closer to the F1 samples, depicting that these environments share some similarities. The close distance between F2c and F1a-c shows that the structure of the community of microbes in F2c shares more similarities with F1a-c than F2a and F2b. Likewise, the close distance shown between R2c and R1a-b depicts that the microbial community structures in this environment (R2c) are similar ( Figure 5). Thus, the closeness of F2c and R2c with the rhizosphere soils (F1a-c and R1a-c) confirms that rhizosphere microorganisms are related to those in the surrounding soils and that the surrounding soil microbes could impact the structure of the microbial communities in the rhizosphere [36,37].

Influence of Soil Properties on Microbial Community
To determine the relationship between the measured soil properties/environmental variables and the distribution of bacteria, fungi and archaea phyla, canonical correspondence analysis was performed using all the measured physical and chemical soil properties ( Table 1). The CCA plot showed that the microbial community composition was dependent on the soil's properties. The vector length of phosphorus (P) showed a positive correlation with Actinobacteria, Ascomycota, Basidiomycota, Bacteroidetes, Thaumarchaeota, Proteobacteria, Crenarchaeota and Firmicutes, and the vector length of nitrate (N-NO 3 ) was seen to correlate positively with Chloroflexi, Korarchaeota, Blastocladiomycota, Euryarchaeota, Gemmatimonadetes, Neocallimastigomycota and Deinococcus-Thermus, while the vector length of organic matter showed a positive correlation with Nitrospirae, Cyanobacteria, Verrucomicrobia, Chlorobi and Planctomycetes ( Figure 6). The forward selection of environmental variables and the Monte Carlo permutation test with 999 random permutations were applied to determine the environmental factor that best explains the variation in the microbial composition in samples, with phosphorus significantly contributing 56.9% of the variation (p-value 0.04) ( Table 3).

Functional Analysis of Samples
The SEED Subsystem database revealed twenty-eight (28) major functional categories at level 1 for all samples (Figure 7, Table S3). The main categories in the rhizosphere and bulk soil samples included amino acids and their derivatives, protein metabolism, RNA metabolism, cofactors, vitamins, prosthetic groups, pigments, miscellaneous, carbohydrates, fatty acids, lipids and isoprenoids and respiration, as shown in Figure 7. The significantly different functional categories included iron acquisition and metabolism, membrane transport, nucleosides and nucleotides, cell wall and cell capsule, phages, prophages, transposable elements, plasmids, fatty acids, lipids and isoprenoids, respiration, amino acids and derivatives, photosynthesis and nitrogen metabolism. Others such as carbohydrate metabolism, metabolism of aromatic compounds, secondary metabolism, miscellaneous, regulation and cell signaling, phosphorus metabolism and potassium metabolism were not significantly different between the environments.  In F1, the dominant functional groups included protein metabolism, carbohydrates, sulfur metabolism, phosphorus metabolism, iron acquisition and metabolism and nitrogen, secondary and potassium metabolism, while in F2, photosynthesis, dormancy and sporulation, stress response, nucleosides and nucleotides, amino acid derivatives, clustering-based subsystems and fatty acids, lipids and isoprenoids were prevalent. RNA metabolism, respiration and metabolism of aromatic compounds were among the dominant functional groups in R1, while in R2, motility and chemotaxis, virulence, disease and defense, cell wall and cell capsule were observed to be dominant.
Principal component analysis (PCA) was performed in order to determine the distribution of the functional categories across the samples (Figure 8). The first and second axes of the PCA plot explained 98.9% and 0.49% of the variance, respectively, with the PCA vector length showing the functional category that dominates each habitat.  Figure 8, the vectors show the functional/metabolic activities which were distributed and predominate each sample. Variables which correlate with each other are shown by vectors pointing in the same direction, and the samples that are close to each other on the PCA show they share similarities. The angle between the variable vectors which is small shows that the variables are highly associated. Functional annotation at SEED Subsystem level 2 showed potassium metabolism to be the most dominant. The relative abundances of this across the samples were 14.86% (F1), 14.6% (F2), 14.55% (R1) and 14.2% (R2). The unknown function category was the next dominant group across the samples, with a relative abundance of 6.44% in F1, and 6.01%, 6.05% and 6.04% in F2, R1 and R2, respectively (Figure 9)

Alpha and Beta Diversity Indices of the Functional Categories Present in Samples
At the level 1 subsystem, the alpha diversity indices showed no significant difference (p > 0.05) in the functions within the habitats. The Simpson index which was the same for all samples (0.92), showing high diversity within the environments, while the Shannon and evenness indices were within the same range (Table 4), where the evenness index ranged from 0.59 to 0.61, indicating that the distribution of functional categories in the habitats was moderately even. However, the beta diversity analysis showed that significant differences exist in the functional categories between the samples (p = 0.01), with the strength of separation between each environment being substantial (R = 0.58). The principal coordinate analysis performed also showed a distinct separation between the samples, as seen in Figure 10. The location of F1 on the PCoA plot was distinctly far from R1, which depicts a clear difference in the functional categories in the samples, with a difference also being observed between F2 and R2, the samples being separated from each other on the PCoA plot. However, closeness was seen between samples F1b and F1c, as well as R1a, R1b and R1c, depicting similarities in these samples.

Influence of Soil Properties/Environmental Variables on the Microbial Functional Diversities in Samples
To determine the relationship between the environmental variables/soil properties and the microbial functional diversities, canonical correspondence analysis was performed. Forward selection of environmental variables that best explain the functional diversities was conducted, and N-NO 3 , P and K were the parameters selected based on the significance test. The result shows that there were positive correlations between N-NO 3 and DNA metabolism, nitrogen metabolism, dormancy and sporulation, photosynthesis, carbohydrates, fatty acids, lipids and isoprenoids, RNA metabolism, cell division and cell cycle. Potassium (K) correlated positively with clustering-based subsystems, stress response, amino acids and derivatives, metabolism of aromatic compounds, motility and chemotaxis, respiration, sulfur metabolism, potassium metabolism and membrane transport. Likewise, a positive correlation was seen between phosphorus (P) and secondary metabolism, virulence, disease and defense, phosphorus metabolism, cell wall and capsule, nucleosides and nucleotides, iron acquisition and metabolism and protein metabolism ( Figure 11, Table 5). Phosphorus significantly contributed 34.5% of the variation in the samples, and K contributed 38.9% of the total variation, while N-NO 3 explained 26.6% of the variation.

Structural Analysis
From the alpha diversity analysis, no significant difference was seen in the microbial community structure within the environments (p > 0.05-Kruskal-Wallis test). As it is represented in Table 2, the Shannon diversity index (used to characterize species diversity) is greater than 1 in all samples, and this depicts the high diversity of the microbial communities. The evenness index of the microbial community showed a good distribution of the metagenomes, especially at the order level [38]. The principal component analysis showed that each of the sites has distinctive microbial phyla, which accounts for the combined 86.9% variation between the samples (Figure 3). The vector arrows in the PCA depict the phyla that dominate each site. Thus, the microbial group that is more prevalent in each of the sampling sites can be deduced.
The PCoA plot ( Figure 5) showed variations in the microbial community structure of the two rhizosphere soil samples, with clear differences being noted between the rhizosphere samples and their corresponding bulk soil samples. This was buttressed by the conducted analysis of similarities (ANOSIM) which indicated a significant difference (p = 0.01) and separation (R = 0.58) between the samples. This agrees with our assumption that there would be variation in the microbial community structure of the examined soil environments. The separation observed between F1 and R1 in the PCoA plot denotes disparity in the microbial communities, which could be linked to differences in the soil properties [39][40][41].
Proteobacteria are involved in a broad range of activities such as nitrogen, carbon and sulfur cycling necessary for nutrient cycling in plants [47]. Alphaproteobacteria can survive in harsh environments such as those with a low level of nutrients while carrying out their activities, with Caulobacterales, which belongs to this phylum, being capable of inducing nitrogen fixation despite the poor nutrient conditions [48]. Similarly, myxobacteria/Myxococcales, which belong to the Proteobacteria phylum, are capable of surviving in low-nutrient soil by aggregating into fruiting bodies, and until adequate nutrients are available, they can thrive and carry out activities such as protecting the plants from phagocytosis and toxins [49].
Firmicutes are an important bacteria phylum, which can be used as biofertilizers due to their vital role in plant growth promotion, phytoremediation of heavy metals and biocontrol of plant pathogens. These important traits exhibited by the members of this phylum are relevant in crop production and position the organisms in this phylum as those that can be exploited as microbial inoculants and biofertilizers for sustainable agriculture [50]. Actinobacteria also possess attributes that make them important for use as bioinoculants/biofertilizers, and this includes their ability to solubilize phosphorus, zinc and potassium, and their production of iron-chelating compounds and phytohormones such as indole acetic acid. These plant growth-promoting and soil health-enhancing attributes exhibited by Actinobacteria contribute to their relevance for use as biofertilizers [51].
Observed in the metagenomes from this study is the presence of Rhizobium, Bacillus, Burkholderia, Pseudomonas, Haloferax and Aquificales, which occurred more in F1 (maize rhizosphere soil from the first site). These microorganisms are important in instilling a certain degree of tolerance into plants towards drought and other abiotic stresses, such as high temperature, metal toxicity, salinity and chilling injury, which makes them exploitable as microbial inoculants [52,53].
The differences observed in the concentration of total carbon and moisture in the samples could be a major factor that influenced the variation in the microbial structure, as these soil properties have been shown to influence soil microbial abundance and diversity [54,55]. In their study, Liu et al. [54] found total soil carbon and soil water content as part of the main contributors to variations in bacteria and fungi (Actinomycetes and arbuscular mycorrhizal fungi) diversity in the soil. Similarly, Li et al. [55] discovered the significant impact of soil carbon on the diversity of soil bacteria (Deltaproteobacteria and Gammaproteobacteria) and other microbes.
Adeboye et al. [56] noted that soil organic carbon and pH influence the abundance of soil microorganisms, and Table 1 shows that a significant difference (p < 0.05) occurred in the pH of F1 and R1, and between R1 and R2. This difference is capable of influencing the variation in the relative abundance of microbes in the samples according to the report of Zheng et al. [57]. The role of soil carbon in stability and soil fertility is crucial and cannot be overlooked [58]. Carbon has been known to influence the microbial community structure [54], and the higher carbon content observed in the rhizosphere soils compared to their corresponding bulk soils could be linked to the presence of more beneficial microbial phyla such as Firmicutes, Planctomycetes, Verrucomicrobia, Crenarchaeota and Chlorobi observed in these samples.
The canonical correspondence analysis ( Figure 6, Table 4) showed that phosphorus explained 56.9% of the microbial diversity in the samples, and organic matter explained 26.20%, while N-NO 3 contributed about 16.9%. The result of the CCA showed that phosphorus (P), organic matter and N-NO 3 mostly shaped the microbial communities in the samples. Phosphorus (P) is an essential component of plants' "energy unit", being important for the general health and vigor of plants as it converts other nutrients into building blocks with which plants grow. P is also a vital nutrient for plants and the organisms in the soil and is usually correlated with pH, while maintenance of the soil organic matter is important to phosphorus availability. The high pH in F1 and the higher concentration of phosphorus in this sample confirms that pH influences the P availability for soil microorganisms [46]. The high pH in F1 can be linked to the higher phosphorus concentration present in this sample and, hence, its increased abundance of microbes specifically. This study is in line with that of Wang et al. [59], who stated that the addition of phosphorus to soil impacts the abundance and diversity of soil organisms. Similarly, the higher organic matter content in F1 can also influence the abundance of microorganisms in this sample as organic matter allows an increased water storage capacity which provides the level of higher humidity required for the expansion of certain microorganisms [60,61].

Functional Potentials in Maize Rhizosphere
Understanding microbial functions in biogeochemical processes is important in developing and advancing environmentally safe approaches to enhance production in agriculture/maize production for food security by manipulating the soil organisms. In this study, metabolic/functional processes in maize rhizosphere and bulk soils were examined. No significant difference was observed in the alpha diversity indices of the functional categories in the habitats (p > 0.05). As it can be seen in Table 4, the Shannon-Wiener and evenness indices of the rhizosphere and bulk soils, which ranged between 2.80 and 2.83, show high functional diversity in the samples, while the evenness index (0. 59-0.61) shows that in all samples, the functional categories were well distributed [62].
The principal coordinate analysis conducted showed distinct differences between the habitats, which was illustrated by the position of the samples on the PCoA plot ( Figure 10). A significant difference was noted in the rhizosphere soils (F1 and R1) and the bulk soils (F2 and R2). The large distance between locations of the rhizosphere samples (F1 and R1) illustrated dissimilarities in their functional groups, with variation also being indicated in the bulk soil samples (F2 and R2), which was represented by the distance between them. Axes 1 and 2 of the PCoA explained 82.85% and 10.72% of the variation across the samples, respectively.
Similarly, the beta diversity analysis conducted using one-way analysis of similarity (ANOSIM) showed significant differences between the samples (p = 0.01), and the R value of 0.58 which depicts the strength of separation between the samples showed a strong dissimilarity between them, which agrees with our assumption that there will be differences in the functional categories across each environment.
Principal component analysis was performed to test the variations in the microbial functions from each site, the result showing that each sample possesses a distinct functional/metabolic profile (Figure 8). PCA axes 1 and 2 explained 98.9% and 0.49% of the variation, respectively, in the functional categories. The length of the vector arrows showed the strength of dominance of the functions in the samples (the functions on the longest vector length are those which are majorly performed by the microbes in each environment). On the PCA plot, functional categories such as disease virulence and defense, nitrogen metabolism, secondary metabolism, sulfur metabolism and potassium and phosphorus metabolism placed sample F1's microbiome as distinct from the microbiome of F2, R1 and R2 (Figure 8). This implies that its associated microbiome majorly helps with these functions. This is supported by the prevalence of Aquificales, Bacillales, Myxococcales, Pseudomonadales, Sordariales and Bdellovibrionales in the sample [63] ( Figure 4A). These organisms are noted for the above-mentioned functions and other significant contributions to plant growth [47,64].
The functional groups involved in R1 included motility and chemotaxis, respiration and metabolism of aromatic compounds, which can be linked to the abundance of Acidobacteriales, Flavobacteriales and Methanomicrobiales and occurrence of Pseudomonadales in the samples [65][66][67][68][69]. Functional groups such as cell wall and cell capsule, photosynthesis, cofactors, vitamins, prosthetic groups, pigments and DNA metabolism distinguished F2 from other samples (Figure 8), and this is evident in the higher composition of Bacillales, Caulobacterales and Clostridales in this sample. RNA metabolism, clustering-based subsystems and amino acids and derivatives distinguished R2 from other samples in the PCA, which can be linked with the higher composition of Chromatiales, Ktedonobacteriales, Nostocales and Thermotogales in these samples [70][71][72][73][74][75].
Higher functional diversity was seen in microbes in the rhizosphere, which is in line with the study of Garcia-Fraile et al. [76], in which bacteria cells obtained from the rhizosphere contained a significant higher level of alkaline phosphatase, β-glucosidase and dehydrogenase activities than those isolated from the bulk soil. This is because plants alter the diversity of microbes within the rhizosphere to provide a suitable and healthy environment for sprouting and growth. The larger amount of carbon immobilized in microbial biomass indicates that soil organic matter, which is observed to be higher in rhizosphere soils, provides greater levels of more labile carbon in these environments than in the bulk soils [77].
In order to determine the relationship between the functional diversities of microbes in the habitats and the soil properties/environmental variables, canonical correspondence analysis was performed ( Figure 11). The CCA result shows that K, P and N-NO 3 were the soil properties that best explain the diversities in the functional categories in the samples, where K and P explained 38.9% and 34.5% of the variation, respectively, while N-NO 3 explained 26.6% of the variation (Table 5, Figure 11). Soil properties are major factors that drive the diversity and structure of soil microbes, thus influencing the functional potentials in the soil ecosystem [8,15].
The results from the soil property analysis show heterogeneity in the physical and chemical attributes of the soil samples ( Table 1). The pH of samples R1 and R2 was below the ideal range of 6-7.5 [78], while that of samples F1 and F2 lies within the normal range, thus indicating the availability and balance of other important soil nutrients [79]. The higher pH observed in F1 in our study can be linked to the higher abundance of other nutrients (such as phosphorus, potassium and sulfur) observed in this sample as compared to others, hence the dominance of more functional categories in F1 than other samples.
Other soil properties such as N-NO 3 , K and N-NH 4 , as observed in the samples, were within the range needed for sufficient plant growth, with these findings being in line with those of Babalola et al. [80]. The higher abundance of organisms such as Pseudomonales, a phosphate-solubilizing bacteria, in F1 could be linked with the raised phosphorus content in the sample [81,82].
The larger amount of carbon immobilized in the microbial biomass indicates that soil organic matter, which was observed to be higher in the rhizosphere soils, provides greater levels of more labile carbon in this environment than in the bulk soils [77]. From this study, higher functional categories were shown in the rhizosphere samples (Figure 7), which is in line with our hypothesis. This is expected due to the higher concentration of important soil nutrients in this region and the root exudates being released into the rhizosphere, which serve as a source of energy for soil microbes to carry out diverse anabolic and catabolic processes.

Conclusions
The use of the shotgun metagenomics approach in this study provided an advantage over the limitations of culture-dependent procedures. It provided a comprehensive and simultaneous study of the structure (at the phylum, class, order, family and genus levels) and functional/metabolic activities of bacteria, archaea and fungi in the maize rhizosphere. Moreover, with the use of this culture-independent technique, bioinformatics tools were applied to determine the environmental variables that best influenced the structure and function of these microbes.
This study revealed the taxonomic composition and metabolic activities of microbes in maize rhizosphere and bulk soils. An abundance of plant growth-promoting microbes and other beneficial bacteria which function in the resistance to abiotic stress was discovered across the sites, although with more prevalence in the rhizosphere soil obtained from the first sampling site. Due to the observed correlation between certain soil properties such as phosphorus and the higher microbial abundance in the rhizosphere samples, especially F1, it can be deduced that the physical and chemical properties of the soil influence the soil microbial abundance and diversity. A distinct difference was observed in the functional categories and nutrient pathways of the soil samples. The canonical correspondence analysis performed showed that soil properties such as P, N-NO 3 and organic matter influenced the microbial diversity, with P explaining the major variation. This agrees with our assumption that there would be variation in the community structure of the examined soil samples as a result of the differences in the soil properties. The presence of unclassified sequences, which were more abundant in F1, signifies the presence of vital novel species in the sites. Further studies are hereby recommended in mining these microorganisms which are yet to be identified, as they could lead to the discovery of important microbes lurking in the maize rhizosphere and their potentials.
Moreover, the studied maize rhizosphere and bulk soils consist of similar functional categories with varying relative abundances. Though the rhizosphere soils have higher functional hits, the alpha diversity analysis showed no difference in the functional categories in the habitats, while the beta analysis showed significant differences across the samples. The influence of the soil physical and chemical properties on the abundance of the functional groups across the habitats was also observed. Hits with unknown functions were discovered in each environment, with higher abundance in F1. Hence, further investigation into the functional categories in maize rhizosphere soils is recommended in order to unveil the unknown groups and their ecological roles which could be relevant/crucial not only in the biogeochemical cycle but also in industrial and biotechnological processes.
This study provided insight into the structure and functional potentials of the community of microbes in the maize rhizosphere. The understanding of the microbes which help in the functioning of the rhizosphere ecosystem would assist in the development of strategies which can be used to increase agricultural production through manipulation of these soil microorganisms. This study also revealed that the maize rhizosphere is endowed with microorganisms such as Myxococcales, Burkholderia, Haloferax and Aquificales which are important in plant growth, helping plants to survive in low-nutrient soil, and in the phytoremediation of heavy metals, and production of iron-chelating compounds and phytohormones. These microorganisms, when isolated, can serve as important bioinoculants which can be used for increased crop production and sustainable agriculture.
Moreover, this study shows the importance of nitrogen, phosphorus, potassium, organic matter and N-NO 3 , as these are the major environmental variables that shaped the microbial community structure and also influenced the metabolic activities of microbes in the maize rhizosphere microbiome. Thus, the soil environment cannot afford to be lacking the accurate proportion of these soil nutrients if there would be optimum production.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/su13148079/s1, Table S1: Sequence information, Table S2: Diversity evaluation of bacteria (at phyla level) in soil samples, Table S3  Informed Consent Statement: Necessary permissions were obtained and granted from field owners from which samples were taken. This study was also approved by the Faculty of Natural and Agricultural Sciences, and the ethics committee of the North-West University (ethical clearance no: NWU-01191-19-A9).

Data Availability Statement:
The raw sequence data were deposited and are publicly available from the NCBI SRA dataset under the accession numbers PRJNA678469 and PRJNA678475. Acknowledgments: O.O.B. acknowledges the National Research Foundation for grants (UID123634 and UID 132595), which funds several research activities in her laboratory.

Conflicts of Interest:
The authors declare no conflict of interest.