Assessing the Impact of Rice Cultivation and Off-Season Period on Dynamics of Soil Enzyme Activities and Bacterial Communities in Two Agro-Ecological Regions of Mozambique

Abstract

Soil ecosystem perturbation due to agronomic practices can negatively impact soil productivity by altering the diversity and function of soil health determinants. Currently, the influence of rice cultivation and off-season periods on the dynamics of soil health determinants is unclear. Therefore, soil enzyme activities (EAs) and bacterial community compositions in rice-cultivated fields at postharvest (PH) and after a 5-month off-season period (5mR), and fallow-fields (5-years-fallow, 5YF; 10-years-fallow, 10YF and/or one-year-fallow, 1YF) were assessed in two agroecological regions of Mozambique. EAs were mostly higher in fallow fields than in PH, with significant (p < 0.05) differences detected for β-glucosidase and acid phosphatase activities. Only β-glucosidase activity was significantly (p < 0.05) different between PH and 5mR, suggesting that β-glucosidase is responsive in the short-term. Bacterial diversity was highest in rice-cultivated soil and correlated with NO₃⁻, NH₄⁺ and electrical conductivity. Differentially abundant genera, such as Agromyces, Bacillus, Desulfuromonas, Gaiella, Lysobacter, Micromonospora, Norcadiodes, Rubrobacter, Solirubrobacter and Sphingomonas were mostly associated with fallow and 5mR fields, suggesting either negative effects of rice cultivation or the fallow period aided their recovery. Overall, rice cultivation and chemical parameters influenced certain EAs and shaped bacterial communities. Furthermore, the 5-month off-season period facilitates nutrient recovery and proliferation of plant-growth-promoting bacteria.

1. Introduction

Rice (Oryza sativa) is one of the most cultivated cereal grains globally, only behind maize (Zea mays) and wheat (Triticum spp.) [1]. It is a staple food in several developed and developing nations and constitutes a strategic trade commodity globally [2,3]. It is estimated that by the year 2027, the global rice market will be a 274 billion dollars venture [4]. Such market projections are based on continued government support, technological advancements, and a sustained global expansion of land areas cultivated with rice. However, the likelihood of expanding rice-cultivated land areas is threatened by a scarcity of water resources, competition from other land-use demands given rapid urbanization and industrialization, as well as the exorbitant cost of establishing new cropping fields suited for rice production, especially in Sub-Saharan Africa and Latin America [5,6,7]. Furthermore, factors such as climate, poor agricultural practices, soil health and quality, seed quality, water availability and limited scientific knowledge constitute a major constraint to global rice production [6,8,9,10,11].

Among such factors, soil health, which is a vital component for crop productivity and proper soil functioning, is often neglected during local rice cultivation [6,12]. Soil health embodies the interaction and balance between soil’s physical, chemical, and biological components [13,14]. The soil biological component comprises microflora and fauna, which synergistically sustain proper soil ecosystem functioning by driving organic matter decomposition, soil carbon sequestration and nutrient mobilization, and supporting soil food webs [15,16,17,18,19]. In particular, soil microbial communities and their interactions in the plant rhizosphere are responsible for the massive global nutrient transfer, including carbon and nitrogen mineralization [16,19,20], soil stabilization through bio-binding (aggregation) [21,22], suppression of pathogens through competitive exclusion [23,24,25] and modulation of plant immunity to pathogens [26]. Such functions of soil microbial communities make them important components for improving crop yield. Unfortunately, studies have shown that soil microbial communities and, consequently, their vital roles in promoting soil health and crop productivity, can be negatively influenced by agricultural practices, such as continuous cropping [27,28,29], long-term fertilization [30], organic amendments [31,32], irrigation and other soil management practices (e.g., no-to and conservation agriculture) [33,34,35,36,37].

In Mozambique, smallholder farmers commonly cultivate rice without rotation under rainfed farming systems in predominantly clayey and hydromorphic soils. Such rainfed rice farming systems represent a suitable model to study fundamental aspects of microbial ecology, such as diversity, structure, and dynamics of microbial communities and relationships between microbial groups [38,39,40]. However, there is a paucity of in-depth studies investigating differences in rice field soil microbial communities between geographical locations, seasonal time points (e.g., harvest vs. pre-planting) and cultivation status (e.g., uncultivated, fallow and actively cultivated). Such studies are important for understanding soil health in rice-cultivated lands, gaining insights into how land management and cultivation practices affect soil microbial communities and improving rice productivity through biofertilization.

Therefore, the present study investigated soil enzyme activities and bacterial communities in selected rice fields in two agro-ecological regions of southern Mozambique, namely Chokwe and Umbeluzi [41]. We hypothesized that: (1) soil chemical properties, enzymes and bacterial communities are affected by rice cultivation (i.e., fallow or active cultivation), and (2) that a five-month off-season period between harvest and planting aids the recovery of soil health determinants. The study aimed to understand the effect of rice cultivation on soil bacterial communities, determine possible temporal changes in soil bacterial communities, and elucidate the relationships between soil physicochemical properties and bacterial communities.

2. Materials and Methods

2.1. Study Area Description

Rice-growing fields in two Experimental Stations of Mozambique’s Agronomic Research Institute located in Umbeluzi, Maputo (26°03′ S 32°23′ E) and Chokwe, Gaza (24°32′ S 33°00′ E) were selected for the study. The study areas area is a shrub savannah located at an altitude of approximately 200 m above sea level. Daily temperature for the study areas ranges between 23 and 26 °C in the rainy/hot season (October to March) and between 17 and 23 °C in the dry/cold season (April to September). The climate of Umbeluzi is characterized as semi-arid, with an average annual precipitation of 679 mm and irregular rainfall patterns during the rainy season [41]. Chokwe has a semi-arid climate with an average yearly precipitation of about 620 mm and an enormous precipitation variation along the year (precipitation occurs essentially from November to March) [41]. The soils in Umbeluzi are of alluvial and basalt origin, while the soil textural class is clay loam with moderate drainage [42]. In contrast, soils in Chokwe are predominantly clayey soils [43].

The rice fields in both locations were cultivated with Makassane rice variety and regularly irrigated. During the land preparation period, fertilization with NPK 12–24–12 at a rate of 30–50 kg/ha was applied, while during the rice-growing period, urea (46% of nitrogen) was applied at a rate of 90 kg/ha. For weed control, chemical pesticides were used. Each rice planting season begins in October and ends in May of the following year, with a five-month period between harvest and planting for the next rice-growing season.

2.2. Soil Sampling

Soil sampling was conducted in the 2015/2016 season (Table 1) on sites with different rice cultivation statuses or time points (RCST) in the two study locations. Soil samples were collected from rice-grown fields at postharvest (PH) and from adjacent (to the rice-cultivated field) areas that had been fallow for at least five years in Chokwe (herein referred to as “5YF”) and for at least ten years in Umbeluzi (referred to as “10YF”). In Umbeluzi, a one-year fallow field (herein referred to as “1YF”) was also sampled (additional fallow sites were unavailable in Chokwe). The fallow fields were approximately 300 m and 100 m to the north of the rice-cultivated fields in Umbeluzi and Chokwe, respectively. In both locations, the predominant vegetation type in the fallow fields were shrubs (up to 40%), and the rest comprised grasses and broad leaves. For an investigation into the possible recovery of soil health during the five months between harvest and planting of the following (2016/2017) season, additional soil samples were collected from the rice grown fields shortly before planting in October 2016 (designated as “5mR”) (Table 1).

Bulk soil samples were randomly collected from a depth of 10 cm using a sterile auger. Three rice-cultivated plots were sampled in each location. Three representative composite samples comprising three simple random soil cores were obtained from each of the three plots. A total of nine bulk samples were collected per field. Soil subsamples for bacterial community analysis were stored in sterile 50 mL polypropylene (falcon) tubes, immediately placed on ice and frozen (−20 °C) before DNA extractions. Soil subsamples for physicochemical analyses were air-dried, ground and passed through a 2 mm sieve before analyses.

2.3. Chemical and Enzyme Activities in Bulk Soil

Chemical properties, including pH (H₂O), electric conductivity (EC), soil organic carbon (SOC), moisture, Olsen available phosphorus (P-Olsen), extractable potassium (Egner Rhiem K), nitrate (NO₃⁻), ammonium (NH₄⁺) and total nitrogen were determined using standard protocols [44,45,46].

Soil enzyme activities, including β-glucosidase, acid- and alkaline -phosphatases, as well as urease activities, were assessed from moist soils using standard methods. Briefly, β-glucosidase, acid phosphatase and alkaline phosphatase (phosphomonoesterases) activities were determined according to the protocol described by Tabatabai [47] using p-nitrophenyl-β-D-glucopyranoside (pNG) as the substrate for β-glucosidase and p-nitrophenyl phosphate (pNP) as the substrate for both phosphomonoesterases. The urease activity was determined according to the method described by Kandeler [48], using the unbuffered option. For all the enzyme activities, absorbance values of the extracts were determined in a segmented flow analyzer system with a dialysis pre-step to remove interferences of color and microparticles. The activities of β-glucosidase, acid- and alkaline -phosphatases were expressed in µg p-nitrophenol h⁻¹ g⁻¹ dry soil, while urease activity was expressed in mg N–NH₄⁺ 2 h⁻¹ kg⁻¹ dry soil.

2.4. Bacterial Community Analyses of Rice-Field Soils

2.4.1. DNA Extraction and Partial 16S rRNA Gene Sequencing

Total community DNA was extracted from 0.25 g soil samples using the power soil DNA isolation kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. DNA integrity and concentrations were determined by electrophoresis on 0.5% agarose gel and fluorometric quantification using a fluorometer (Qubit 2.0, Invitrogen, CA, USA), respectively. Following DNA extraction, amplicon libraries of the hypervariable V4 region of the 16S rRNA gene were prepared according to the Illumina MiSeq protocol using the primer pair 515 F (5′-GTGYCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACNVGGGTWTCTAAT-3′) [49]. Each forward and reverse primer contained Illumina forward and reverse overhang adapters, respectively (Illumina Inc., Foster City, CA, USA). Briefly, the library preparation protocol entailed a first “amplicon PCR” step involving 12.5 ng DNA, 0.2 µM of each forward and reverse primers, 12.5 µL of double-strength KAPA HiFi HotStart ready Mix [0.5 U DNA polymerase, 0.3 mM dNTPs, 2.5 mM MgCl₂] (Kapa Biosystems, Wilmington, MA, USA) and nuclease-free water (Thermo Fisher Scientific, Waltham, MA, USA) in a final reaction volume of 25 µL. PCR was performed in a thermal cycler (SimpliAmp, Applied Biosystems, Foster City, CA, USA). The PCR condition was an initial step of 95 °C for 3 min, followed by 25 cycles of 98 °C for 20 s, 55 °C for 10 s and 72 °C for 10 s, as well as a final extension of 72 °C for 1 min. PCR amplicons were purified with AMPure XP magnetic beads (Beckman Coulter, Brea, CA, USA) and indexed using Nextera XT primers (Illumina Inc., Foster City, CA, USA). Indexed amplicons were purified, quantified, normalized, pooled in equimolar concentrations and denatured in 0.2 N NaOH. Denatured libraries, along with PhiX sequencing control, were then loaded onto the MiSeq flow cell using the MiSeq v3 reagent kit for a paired-end (2 × 300 bp) sequencing run. Sequencing was performed at the Biotechnology Platform of the Agricultural Research Council, Pretoria, South Africa.

2.4.2. Bioinformatic and Diversity Analyses

Paired-end sequence reads were demultiplexed using the MiSeq reporter software (Illumina Inc., CA, USA) and checked for quality using FastQC v.0.11.9 (Babraham Institute, Cambridge, UK). Paired-reads were trimmed at both 5′ and 3′ends to eliminate poor quality nucleotides, denoised, merged, depleted of chimeric sequences and clustered into amplicon sequence variants (ASV) (operational definition for a species) by using the DADA2 denoiser [50] incorporated into QIIME 2 [51]. The resulting ASV count table was depleted of singletons, and representative sequences taxonomically classified using a trained classifier (V4 region of the 16S rRNA) of the SILVA reference (Release 132) [52]. After this, contaminant sequences (e.g., archaea, mitochondria and chloroplast) were removed, and the resulting ASV table normalized to 11,000 sequences per sample before computing diversity analyses in QIIME 2 and visualizations using ggplot2 v.3.2.1 [53], vegan v.2.5-5 [54] and labdsv v.1.8-0 [55] packages of R software v.4.0.2 [56].

2.5. Statistical Analyses

Except otherwise indicated, all statistical analyses were performed using R software v.4.0.2 [56]. Comparisons between RCST were performed only within locations. Comparisons between locations were not performed due to differences in the number of sampled plots, soil type, differing fallow periods between locations and underlying variations in agronomic practices, among others. To test the first hypothesis, data (physicochemical properties, enzyme activities and alpha diversity metrics) were compared between 10YF, 1YF and PH in Umbeluzi and between 5YF and PH in Chokwe using the Kruskal–Wallis rank-sum test and the Mann–Whitney U test, respectively. To test the second hypothesis, the Wilcoxon rank-sum test was used to compare values between PH and 5mR in both locations. Spearman’s rank-order correlation was also performed to test the relationship among soil properties, enzyme activities and alpha diversity.

For beta diversity analyses of microbiome data, relative ASV counts were first log₁₀-transformed (log₁₀ (x) +1, where x > 0) [57] using the “decostand” function in the vegan R package. Thereafter, differentiation among RCST in multivariate space was determined based on Bray–Curtis dissimilarities and visualized by both principal coordinates analysis (PCoA) and dendrograms using ape v5.4.1 [58] and dendextend v1.12.0 [59] R packages, respectively. Determination of significant differences in multivariate space was based on permutational multivariate analysis of variance (PERMANOVA) using the vegan package. Further permutational test for homogeneity of multivariate dispersion (PERMDISP) was performed to test the significance of within-group differences. A post hoc test on significant PERMANOVA (FDR-adjusted p ≤ 0.05) was also performed using the function “pair-wise.adonis” in the vegan package.

To determine differentially abundant phylotypes (false discovery rate (FDR)-adjusted p < 0.05, least discriminant analysis (LDA) score > 4.0) among RCST at the phylum and genus taxonomic ranks, the least discriminant analysis (LDA) effect size (LEfSe) [60] was performed using MicrobiomeAnlayst [61]. To determine the physicochemical properties that influence bacterial communities of these soils, a constrained redundancy analysis (RDA) was performed on Hellinger-transformed environmental and bacterial data using ampvis2 and vegan R packages. The significance of the constrained RDA axis was determined based on a permutation test using the function “anova.cca ()” in the vegan package.

Furthermore, we constructed a co-occurrence network using FlashWeave [62]. Prediction of direct associations was based on a graphical model inference (local-to-global learning approach) [63] using a modified semi-interleaved HITON-PC algorithm [62]. False discovery rate (FDR) adjustment, compositionality correction of the dataset using the centered log-ratio (clr) approach (genus level phylotypes and environmental variables) and data filtering to remove phylotypes occurring in less than 30% of the samples, and with relative abundance <1% were performed before network generation. Thereafter, the generated networks were visualized using Gephi (https://gephi.org/; accessed on 2 March 2021), while the topological properties of the networks were determined using igraph [64]. The network nodes (phylotypes and chemical properties) were clustered into modules (communities) using a multilevel modularity optimization algorithm [65]. The nodes with the highest alpha and betweenness centrality scores (top 10 nodes present in both centrality metrics) were considered network and/or module hubs (keystone species). Furthermore, to determine if the relative abundance of a phylotype (network nodes) significantly affects their importance in the ecological network, the centrality metrics were z-score transformed and subjected to permutational multivariate analysis of variance using the “adonis()” function in the vegan package [54].

2.6. Deposition of Sequences

Sequence data generated for this study are available in the Sequence Read Archives (SRA) of the National Centre for Biotechnological Information (NCBI) under the accession number PRJNA511657.

3. Results

3.1. Chemical Characteristics of Bulk Soil in Relation to Rice Cultivation and Fallow Periods

The pH of rice-cultivated field soils was higher than those of fallow field soils in Umbeluzi and Chokwe (Table 2). In Umbeluzi, significant (p < 0.05) differences in pH, electrical conductivity (EC), ammonium nitrogen (NH₄⁺–N) and nitrate (NO₃⁻–N) were observed between fields with different cultivation statuses (Table 2). NO₃⁻–N and NH₄⁺–N were significantly (p < 0.05) higher in the ten years fallow field (10YF) than in both one-year fallow (1YF) and rice cultivated field at postharvest (PH). However, trends in pH, total N and P-Olsen indicated higher values in the rice-cultivated field (PH) than in 10YF and 1YF; values of pH, total N and P-Olsen in 1YF were intermediate between 10YF and PH (Table 2). In Chokwe, significantly (p < 0.05) higher SOC and NH₄⁺–N were observed in 5YF than in PH, whereas NO₃⁻–N was significantly (p < 0.05) higher in PH compared to 5YF (Table 2).

Between PH and 5-months off-season recovery (5mR) time points in Umbeluzi, a significant (p < 0.05) increase was observed in only EC, NO₃⁻–N and K, whereas a significant (p < 0.05) reduction was observed in the only NH₄⁺–N and P-Olsen. In Chokwe, EC, C, NO₃⁻–N, and total N increased significantly (p < 0.05) between PH and 5mR sampling time points, whereas pH and moisture content reduced significantly (p < 0.05) during the five months (Table 2).

3.2. Enzyme Activities and Relationship with Sol Properties

With the exception of urease activity in Chokwe, enzyme activities were mostly higher in the 10-year fallow soils compared to the rice-cultivated (PH) soil in Umbeluzi and Chokwe (Table 3). However, significant (p < 0.05) differences were detected in only β-glucosidase and acid phosphatase activities among RCST in Umbeluzi and Chokwe, respectively (Table 3). The β-glucosidase activity was significantly (p < 0.05) higher in 1YF than in PH in Umbeluzi, while acid phosphatase activity was significantly (p < 0.05) higher in 5YF compared to PH in Chokwe. Only β-glucosidase activities were significantly different between PH and 5mR in both Umbeluzi and Chokwe, suggesting that of all enzymes assayed in this study, β-glucosidase is the most responsive (dynamic) following the five months off-season period (Table 3).

For the identification of soil chemical properties that drive the significant differences in β-glucosidase and acid phosphatase activities observed in both study locations, a Spearman rank-order correlation was performed. In Chokwe, a positive significant (p < 0.05) correlation was observed between pH and β-glucosidase activity in the PH site (rho = 0.73, p = 0.031), whereas a negative association was observed between pH and β-glucosidase activity in the corresponding 5YF site (rho = −0.83, p = 0.008) (Figure 1a,e). In addition, significant (p < 0.05) positive associations between β-glucosidase activities and SOC were observed for PH in Umbeluzi and 5YF in Chokwe (Figure 1b,g). Interestingly, in the 10YF field in Umbeluzi, β-glucosidase activity correlated negatively with EC (rho = −0.95, p < 0.001), but correlated positively with EC in PH (rho = 0.87, p = 0.005) (Figure 1c). Similarly, acid phosphatase activity was positively correlated with P-Olsen in 1YF soils in Umbeluzi (rho = 0.78, p = 0.017) (Figure 1d), but negatively correlated with P-Olsen for 5mR in Chokwe (rho = −0.80, p = 0.014) (Figure 1h). In Chokwe, a negative correlation was observed between acid phosphatase activities and pH in soils from 5YF (rho = −0.73, p = 0.031) (Figure 1f), while a positive correlation was observed between acid phosphatase activity and soil electrical conductivity in soils from 5mR (rho = 0.83, p = 0.008) (Figure 1i).

3.3. Bacterial Diversity and Community Structure in Fallow and Rice-Field Soils

3.3.1. Bacterial Diversity and Community Differentiation of ASVs

Following quality control (i.e., read trimming, low-quality filtering, denoising and chimeric sequences removal), a total of 3,606,164 high-quality sequence reads were clustered into ASVs (data not shown). After taxonomic classification of ASVs and elimination of non-bacterial phylotypes, rarefaction of sequence reads to 11,000 reads per sample was sufficient to describe the bacterial diversity in all soil samples (Supplementary Figure S1). In both Umbeluzi and Chokwe, the number of ASVs unique to the rice-cultivated soil (i.e., PH) was higher than in fallow fields (10YF and 5YF) (Figure S2a,b). In Umbeluzi, the number of ASVs in the sampled soils increased in the order of 10YF < 1YF < PH (Figure S2a). Interestingly, a reduction in total ASVs and unique ASVs was observed between PH and 5mR time points in both locations (Figure S2c,d).

The observed ASVs and Chao1 revealed higher species richness in PH compared to fallow fields in both locations, with significant (Bonferroni-adjusted p < 0.05) differences only between PH and 10YF in Umbeluzi (Figure 2a,b,e,f). Similarly, in Umbeluzi, ASV diversity (based on Shannon–Wiener index) was significantly higher in PH than in 10YF with no significant differences (p > 0.05) observed for PH and 5YF in Chokwe (Figure 2c,g). In addition, in Umbeluzi, the proportionality of species (evenness) in the soil bacterial community of PH was significantly higher than in 10YF, with species evenness in 1YF intermediate, but not significantly (p > 0.05) different from 10YF and PH (Figure 2d). This observation suggests that rice cultivation does not exert selective pressure on bacterial species compared to sites that have been fallow for more than one year. Except for significant (Wilcoxon, p < 0.05) differences in species evenness between PH and 5mR in Chokwe, no significant (p > 0.05) differences in species richness, diversity and proportionality were observed between postharvest and pre-planting time points (5 months period) in both locations (Figure 2).

Due to the significant difference observed in species richness and diversity metrics between RCST in Umbeluzi, we performed a Spearman’s rank correlation to obtain insight into the relationship between soil properties and species diversity in Umbeluzi. Correlation analysis revealed that species diversity and richness are directly related with selected soil properties, including NO₃⁻–N (r =0.48, p = 0.0115), NH₄⁺–N (r = 0.68, p = 0.0001) and electrical conductivity (r = 0.49, p = 0.012).

Bray–Curtis dissimilarity between soil bacterial community structure of RCST in Chokwe and Umbeluzi revealed differentiation between fallow fields (10YF, 5YF or 1YF) and rice-cultivated fields (PH or 5mR) in ordination space (Figure 3). While differentiation between different sampling timepoints (PH and 5mR) did not show distinct differentiation in multivariate space in Umbeluzi (Figure 2a,b), differentiation was observed between RCST in Chokwe (Figure 3c); this was also largely supported by the UPGMA cluster dendrogram in Figure 3d. The observed differences in multivariate space were significant in both Umbeluzi (PERMANOVA R² =14.21%, p = 0.001; PERMDISP p = 0.106) and Chokwe (PERMANOVA R² = 12.97%, p = 0.001; PERMDISP, p = 0.223). Pair-wise comparisons (post hoc PERMANOVA) revealed that the bacterial community structures are significantly different (FDR-adjusted p < 0.05) between RCST in both geographical locations (Table S1). Interestingly, the observed significant differences were not only between rice-cultivated (PH) and fallow soils (10YF, 5YF or 1YF) but also between postharvest (PH) and pre-planting (5mR) time points (Table S1).

3.3.2. Taxonomic Diversity and Differentially Abundant Phylotypes

The ASVs detected in all RCST spanned 47 phyla (Umbeluzi = 46 phyla; Chokwe = 43 phyla). Proteobacteria was relatively more abundant in all, but the ten-year fallow (10YF) site in Umbeluzi, where Actinobacteria was the most abundant phylum (Figure 4). Other relatively abundant (>10%) phyla in at least one RCST across both Umbeluzi and Chokwe include Acidobacteria, Bacteroidetes and Chloroflexi. Notably, in Umbeluzi, Actinobacteria was relatively more abundant in 10YF than in 1YF and the cultivated field (PH or 5mR). In contrast, Bacteroidetes were notably higher in the cultivated field (PH) in Umbeluzi than in the corresponding 10YF (Figure 4a). During the 5-month off-season period in Umbeluzi, the relative abundance of Actinobacteria, Chloroflexi, Germatiomonadetes, and Rokubacteria phyla increased, whereas those of Bacteroidetes, Verrucomicrobia markedly reduced in the same period (Figure 4a). In Chokwe, higher relative abundances of Proteobacteria and Acidobacteria were observed in the cultivated PH soil compared to 5YF, whereas markedly higher relative abundances of Actinobacteria and Chloroflexi were observed in 5YF compared to the cultivated site (PH). During the 5-month off-season period in Chokwe, increased relative abundances of Chloroflexi and Actinobacteria were observed, whereas the relative abundances of Bacteroidetes and Proteobacteria reduced in the same period (Figure 4a).

Sixteen and seven phylum-level features were differentially (FDR-adjusted p < 0.05, LDA > 4.0) abundant among RCST in Umbeluzi and Chokwe, respectively (Figure 4b,c). In Umbeluzi, Actinobacteria and Firmicutes were differentially most abundant in 10YF (Figure 4b). In contrast, Proteobacteria, Bacteroidetes, Verrucomicrobia, Cyanobacteria, Spirochetes and Fibrobacteres were differentially abundant in rice-cultivate soils at the PH time point, while Acidobacteria, Chloroflexi, Rokubacteria, Latescibacteria, Nitrospirae and Zixibacteria were differentially more abundant in the rice-cultivated field at the 5mR time point compared to other RCST (Figure 4b). In Chokwe, only Patescibacteria was differentially more abundant in the fallow field (5YF) compared to other RCST (Figure 4c), while 57% (4/7) of the differentially abundant phyla were most abundant in the rice-cultivated field at the postharvest (PH) time point (Figure 4c). In contrast, Actinobacteria and Armatimonadetes were differentially most abundant in 5mR (Figure 4c).

At the genus taxonomic rank, 24 classifiable phylotypes were relatively abundant (≥1%) in soils from at least one RCST (Figure 5a). Of these, Anaeromyxobacter, Bacillus, Haliangium, Marinobacter, MND1 (family: Nitrosomonadaceae), Nocardioides, Pseudomonas, RB41 (family: Pyrinomonadaceae), Rubrobacter and Sphingomonas were among the relatively more abundant phylotypes across RCST in both locations. In Umbeluzi, the genera Bacillus, Sphingomonas and Pseudomonas were more abundant in 10YF, 1YF and PH, respectively, whereas, in Chokwe, MND1 (family: Nitrosomonadaceae) was the most dominant in both 5YF and PH (Figure 5a). Between the PH and 5mR sampling time points (5-months postharvest fallow period), increased abundance of Marinobacter (+0.27% and +4.97%) and RB41 (+6.15% and +6.6%) was observed in both Umbeluzi and Chokwe locations, while an increase in Sphingomonas (+8.51%) and Geodermatophilus (+4.76%) was observed only in Umbeluzi and Chokwe, respectively. In contrast, a reduction in the relative abundance of MND1 (Nitrosomonadaceae) (−0.08% and −19.24%), Pseudomonas (−7.72% and −11.29%), Haliangium (−2.8% to 6.4%), Ellin6067 (family: Nitrosomonadaceae) (−4.23% and −4.98%), Sideroxydans (−0.82% and −4.77%), Rubrobacter (−0.07% and −3.4%) and Geobacter (−1.7 and −3.3%) was observed over the 5-month off-season period for both Umbeluzi and Chokwe locations (Figure 5a).

Following differential abundance testing, a total of 103 and 33 differentially abundant (Kruskal–Wallis FDR-adjusted p < 0.05, LDA > 4.0) genera were observed between all RCST in Umbeluzi and Chokwe, respectively (Tables S2 and S3). The first 10 (based on LDA score) differentially abundant genera in Umbeluzi and Chokwe are shown in Figure 5b,c, respectively. Notably, most of these top differentially abundant genera, such as Agromyces, Bacillus, Desulfuromonas, Gaiella, Lysobacter, Micromonospora, Norcadiodes RB41, Rubrobacter, Solirubrobacter and Sphingomonas, were more associated with the fallow (10YF, 5YF and 1YF) and 5mR soils. In contrast, only Pseudomonas and MND1 were differentially abundant in PH for Umbeluzi and Chokwe, respectively.

3.4. Relationship between Soil Properties and Bacterial Communities of RCST

The redundancy analysis (RDA) model for the association between bacterial species and chemical properties (environmental cues) in both Umbeluzi and Chokwe was significant (p = 0.001). The constrained variables explained only 12.72% and 11.21% of the bacterial community variations in Umbeluzi (Figure 6a) and Chokwe (Figure 6b), respectively. Permutational test of the significance of the environmental parameters in the RDA model revealed that ammonium nitrate (NH₄⁺–N), nitrate (NO₃⁻–N), pH and EC significantly (p < 0.05) influenced the RCST bacterial communities in Umbeluzi (Table S4). Whereas in Chokwe, nitrate (NO₃⁻–N), total nitrogen, P-Olsen, moisture content, SOC, pH and electrical conductivity (EC) significantly (p < 0.05) influenced the RCST bacterial communities (Table S4). In Umbeluzi, ammonium (NH₄⁺–N) and nitrate (NO₃⁻–N) strongly correlated with the centroid of the 10YF soil samples, whereas in Chokwe, nitrate (NO₃⁻–N) and EC were strongly correlated to the centroid of the 5mR soil samples. Furthermore, EC appeared to be the environmental variable differentiating the bacterial community of PH and 5mR soil from other RCST in Umbeluzi (Figure 6a). In contrast, SOC, NH₄⁺–N and P-Olsen mostly shaped the bacterial community of 5YF in Chokwe (Figure 6b). Notably, most individual bacterial phylotypes (genus taxonomic rank) in RCST were largely uninfluenced by the measured soil chemical parameters (Figure 6a,b). However, in Umbeluzi, Pseudomonas and RB41 correlated with EC, while Gaiella correlated with NH₄⁺–N and NO₃⁻–N. In Chokwe, Ellin6067, Masillia and Agromyces correlated with EC, pH and NO₃⁻–N (Figure 6b).

3.5. Co-Occurrence Networks of Microbes and Environmental Variables Interaction in Fallow and Rice-Cultivated Soils

Ecological networks depicting interactions among microbes and environmental variables for rice-cultivated and fallow fields in Umbeluzi and Chokwe followed a power-law distribution (Kolmogorov–Smirnov p < 0.05) (Figure 7 and Table S5). A total of 150, 414, 104 and 462 significant (FDR-adjusted Fisher’s z-transformation p < 0.05) direct associations (synergistic and competitive) were observed for rice-cultivated fields in Umbeluzi (RCU), rice-cultivated fields in Chokwe (RCC), fallow fields in Umbeluzi (FU) and fallow field in Chokwe (FU), respectively (Table S5). Furthermore, the centrality metrics were independent of the relative bacterial abundance in all the ecological networks (RCU: p = 0.07, RCC: p = 0.45, FU: p = 0.12, FC: p = 0.06) (Figure 7).

Multiple network topological properties indicate that the bacterial co-occurrence pattern in the rice-cultivated fields differed profoundly from those of the fallow fields and that differences in sites also influenced the topological properties of the networks (Table S4). This is particularly true for some features used to infer network complexity (average path length, network diameter, number of edges, average degree). For instance, the network average path length and diameter were higher in the rice-cultivated fields than in the fallow fields (10YF and 5YF) for both Umbeluzi and Chokwe (meaning rice cultivation reduced the network complexity). In contrast, the number of interactions (edges) and average degree was higher in Chokwe than Umbeluzi irrespective of cultivation status (network complexity is not entirely independent of site differences) (Table S4).

Notably, environmental variables did not significantly (p < 0.05) influence all but the RCU and RCC co-occurrence networks. In the RCU network, high centrality metrics were observed for Pseudomonas, Lysobacter, Stenotrophobacter, Nitrospira, Sphingomonas and RB41, suggesting that these phylotypes constitute potential keystone taxa. On the other hand, NO₃⁻–N, EC and pH were the soil chemical parameters that were significantly associated with some bacterial phylotypes (Figure 7a). In the RCC co-occurrence network (Figure 7b), potential keystone (with high centrality metrics) taxa included Geothermobacter, Dechloromonas, Desulfomicrobium, Hydrogenophaga, Haliangium, Chthonomonas, RB41, Bacillus, Bryobacter, Candidatus Koribacter and Nitrospira. Interestingly, iron-oxidizing and reducing bacteria, such as Sideroxydans, Geobacter and Anaeromyxobacter, were among the hubs in three different modules of the RCC network. Apart from SOC, other soil chemical parameters, such as NO₃⁻–N, EC, and pH, were consistent in their influence on the bacterial communities of the rice cultivated fields in both Umbeluzi and Chokwe.

Overall, the bacterial community of RCC seemed to partition into interconnected niches for iron-, nitrogen-, carbon- and sulfur-cycling as SOC, total N and P-Olsen associated significantly with several bacterial species (Figure 7b). Network hubs in FU included Methyloceanibacter, Zeaxanthinibacter, Ohtaekwangia, Skermanella, Desulfovirga and Aridibacter (Figure 7c), while the network and module hubs in FC included MND1, Chthoniobacter, Mycobacterium, Pedomicrobium and Haliangium (Figure 7d).

4. Discussion

Soil management practices in agroecosystems impact the balance between soil physicochemical properties, soil biota diversity and function, and soil health [66,67]. In this study, we investigated the impact of rice cultivation on some soil health indicators by comparing rice-cultivated soils (PH) to fallow fields (10YF, 5YF or 1YF) in two agroecologically distinct regions of Mozambique. Furthermore, the impact of the short 5-month off-season duration between harvest and pre-planting on these soil health indicators was explored.

The soil ecosystem supports several nutrient geochemical cycles, which are important for plant productivity. Of all chemical properties analyzed in this study, only ammonium nitrogen was significantly (p < 0.05) higher in fallow soils than in rice-cultivated fields. This observation suggests that nitrogen cycling is a key feature during rice cultivation and partly supports well-established findings that nitrogen is the most limiting factor for rice production in paddy soils [68,69,70,71]. Specifically, the higher NH₄⁺ contents in fallow soils are probably a consequence of crop absence, with N-uptake by the natural vegetation very reduced. In such situations, N demand is lower, and NH₄⁺ mineralized from organic matter tends to accumulate until its oxidation to NO₃⁻. Differing trends observed for Umbeluzi (e.g., differences in electrical conductivity and NO₃⁻–N were significant) and Chokwe (e.g., differences in SOC was significant) on comparisons of nutrient levels in their respective fallow and rice-cultivated fields are most likely due to wide variation in soil fertility status, soil physical properties and/or soil management practices, such as irrigation, tillage, fertilization types and application rate between the two regions. Studies have shown that different soil management practices and soil fertility status account for variations in nutrient availability in the soil and rice productivity across rice-growing regions [72,73,74,75,76,77,78,79,80]. Thus, the transfer of technological initiatives to farmers across rice-growing regions of Mozambique is crucial for minimizing inherent soil-based productivity restrictions in lower rice productivity areas and for the improvement of technical efficiency among farmers. Importantly, knowledge regarding how tillage practices influence the dynamics and intricacies of the interactions among the physical, chemical and biological components of rice field soils are crucial to achieving soil health and increased rice yields [79,81,82,83,84].

In agreement with the study of Li et al. [85] in paddy rice soil, it was observed that pH increased in the rice-cultivated soil compared to the fallow soils, with differences in pH between cultivated (PH) and fallow fields significant (p < 0.05) only in Umbeluzi (Table 2). This observation, taken together with the observed non-significant (p > 0.05) difference in total N and significantly (p < 0.05) higher nitrate and ammonium levels in fallow fields (10YF and 5YF) compared to PH, suggests high N-mineralization and nitrification (oxidation of ammonium to nitrates) in fallow fields with a concomitant increase in acidity [70,86]. In a study by Yu [87], it was observed that soil pH increased in rice-cultivated paddy soils and that the magnitude of increase was related to the amount of organic matter. However, in the present study, relative to the fallow field, although an increase in pH was observed in the rice-cultivated soils in both Chokwe and Umbeluzi, pH increased, while SOC (from which organic matter is often derived) reduced in the cultivated soil in both locations (Table 2). The disparity in observations between this study and that of Yu [87] mirrors the constant biochemical flux of soil properties in relation to location-specific differences in SOC (or organic matter) availability and initial soil pH, among other factors [87,88]. Importantly, the pH of the rice-cultivated fields at both postharvest and pre-planting time points in both locations are above the optimum (5.5–7.0) for rice cultivation [80,87]; therefore, nutrient availability and crop productivity may be hampered. Thus, pre-planting addition of appropriate soil ameliorant or fertilizer to lower soil pH before rice planting is required [68,69,70].

Indeed, fertilization and irrigation control are crucial factors that influence the mobility and utilization efficiency of nutrients, such as nitrogen and phosphorus, in paddy soils [68,89]. Given that fertilizer was applied in the rice-cultivated fields during the prior planting season and that the fallow fields were not fertilized and irrigated for at least 5–10 years, the comparable levels of total N, P and K across both fallow and cultivated soils at the postharvest time of sampling suggests loss of nutrients in the rice-cultivated soils. Such loss of nutrients may be due to biochemical processes, such as denitrification, ammonia volatilization and clay fixation, as well as nutrient losses in water runoff, deep percolation and uptake by the rice cultivars [90,91,92]. Moreover, the observed significant reduction in ammonium (NH₄⁺–N), nitrate (NH₄⁺–N), and K between postharvest (PH) and pre-planting (5mR) in both Umbeluzi and Chokwe also suggest further leaching, nutrient adsorption (e.g., insoluble phosphates due to binding on clay particles, Al and Fe hydroxides/oxides) and/or oxidation of nutrients during the off-season. These observations agree with previous studies in which significant reduction in the availability of nutrients, such as nitrogen, phosphorus and potassium in paddy soil, were attributed to water runoff, biological activity, leaching, binding and nutrient uptake by plants, among other factors [68,90,93,94].

Enzymatic activities in soils are important for the decomposition of organic matter and mineralization of nutrients; they are useful indicators of soil biological activity in arable soils and deductively soil health [17,95,96]. In general, the enzyme activity values obtained in rice-cultivated soils are comparable to values obtained in rice paddy soils elsewhere [97]. The observed higher enzyme activities in mostly short-term (1YF) and long-term fallow (5YF and 10YF) soil samples compared to the rice-cultivated PH, and 5mR in both study locations suggests higher biological (metabolic) activity in the fallow soils compared to the rice-cultivated soils. This observation is in agreement with an earlier study by Dinesh et al. [98], who reported a decline in microbial activity, enzyme synthesis and accumulation as a result of reduced C turnover and nutrient availability in long-term cultivated and deforested wet tropical forest soils. In contrast, Raiesi and Beheshti [99] observed that the activities of enzymes, such as urease, alkaline phosphatase and acid phosphatase, were not affected by paddy rice cultivation or land-use changes in Northern Iran. However, when expressed in relation to microbial biomass carbon or soil organic carbon (i.e., specific enzyme activity), the authors observed an increase due to rice cultivation [99]. Thus, the lower activities of beta-glucosidase and phosphatase in rice-cultivated soils in both Umbeluzi and Chokwe suggest that rice cultivation, along with the associated reductive effect on nutrient availability, contribute to a reduction in the stoichiometry of these enzyme activities involved in the mineralization of carbon and phosphorus [95,97,98,100,101]. In the case of urease, we observed higher mean activity in the rice-cultivated field in Chokwe compared to the corresponding fallow field (5YF), but the opposite was true for Umbeluzi. Urease catalysis the hydrolysis of urea into ammonia and CO₂; therefore, the higher mean urease activity in the rice cultivated field (at postharvest and pre-planting time point) compared to the fallow-site in Chokwe may be due to high residual urea following its (urea) application during the growing season. Taking the contrasting observation in Umbeluzi into consideration, altogether, the pattern in urease activity and associated soil chemical properties suggest a complex interaction between soil components and the influence of other factors not accounted for in the current study.

Nevertheless, the pattern of changes (i.e., from long-term fallow fields to cultivated fields) in the respective SOC and P-Olsen contents in Umbeluzi and Chokwe, respectively, explains the significant (p < 0.05) differences detected in β-glucosidase and acid phosphatase activities among RCST in Umbeluzi and Chokwe, respectively. Indeed, the suitability of enzyme activities as soil fertility indicators lies in their sensitivity to soil management effects and the often direct relationship with vital soil-quality parameters, such as organic matter, soil physicochemical properties, soil microbial activity and biomass [102,103,104,105,106]. Thus, the observed significant (p < 0.05) changes in only β-glucosidase activities between PH and 5mR soils in both study locations suggest that β-glucosidase is the most responsive (dynamic) of all other enzymes assayed in this study. A similar sensitivity of β-glucosidase was observed in anthropogenically disturbed topsoil in South Africa [107]. Such sensitivity of β-glucosidase is probably due to the need for soil microbes to scavenge for C through the secretion of extracellular β-glucosidase as carbon is an important element for all life forms, including soil biota. As a note, we observed high standard deviation values for enzymatic activities in the study sites. A similar observation has also been reported for undisturbed soils in South Africa [17]. Such an observation is most likely due to spatial heterogeneity consequent on dispersal limitations of nutrients and biomolecules in soil ecosystems [108,109]. Essentially, this observation emphasizes the need for extensive and representative soil sampling in soil-related research.

In this study, next-generation sequencing of the bacterial 16S rRNA gene was used to explore the bacterial diversity and community composition in fallow and rice-cultivated soils. A higher number of unique ASVs were present in the rice-cultivated field than in the fallow fields, and a reduction in the number of unique ASVs was observed between PH and 5mR time points. Altogether, these observations suggest that rice cultivation and associated treatments, such as fertilization, encourage the proliferation of certain species, which are uniquely adapted and functionally relevant in mineralizing nutrients in the rice-cultivated soil environment compared to the fallow soils and to the 5mR soils, which were not actively managed for the short five months off-season duration. Our observations of fluctuating bacterial species richness and diversity in fields with different rice cultivation status or time points is similar to findings from studies investigating bacterial communities in rice agroecosystem during rice plant developmental stages [110,111], under different farming conditions [112] and for wild and cultivated rice cultivars [113].

The observed significantly higher species richness, diversity and proportionality in the rice cultivated field compared to the fallow fields in Umbeluzi (Figure 2a–d) may be due to soil properties and factors mentioned earlier. Similarly, Huang et al. [114] reported that a 19-year continuous application of inorganic fertilizer improved bacterial richness and diversity in rice-cultivated paddy soil. In a related study, Kumar et al. [115] also reported that bacterial diversity (Simpson index) increased with rice cultivation in soils amended with nitrogen fertilizer. According to Jiao et al. [116], the dry-wet conditions during rice cultivation promote the growth of aerobic and anaerobic bacteria, thus improving overall alpha diversity. On the contrary, observed non-significant differences in species richness, diversity and evenness between cultivated and fallow fields in Chokwe (Figure 2e–h) suggest a deviation from observations in Umbeluzi and the aforementioned studies. Such deviations may be due to soil properties and management factors unique to the sampled areas in Chokwe.

The significant differences among bacterial communities of RCST in multivariate space (Figure 3 and Figure 6) indicate that bacterial community functions in the fallow soils differed significantly from those of the rice-cultivated soils. The likely factors that must have resulted in the observed bacterial community compositional changes between the rice-cultivated and fallow soil are the differences in fertilizer management and agronomic practices. Studies have shown that soil bacterial community composition in agroecosystems is influenced by land use and management practices [27,30,31,117]. For example, Wang et al. [30] observed that the stages of rice cultivation (tillering, booting and ripening stages) have a significant impact on bacterial community structure. Though the influence of rice cultivation time points on bacterial community structure has been demonstrated, it is important not to neglect the role of site differences on bacterial communities. For instance, the importance of spatial differentiation is particularly highlighted in the pattern of clustering of the PH and 5mR time points in Umbeluzi (PH and 5mR clustered together) and in Chokwe (PH and 5mR clustered separately). This observation suggests that beyond the influence of cultivation status, differences in soil-type and soil chemical properties between the two study sites are important factors that determine how bacterial species assemble into highly diverse communities. Investigation of bacterial species-environment relationship revealed that ammonium (NH₄⁺), nitrate (NO₃⁻), total nitrogen, P-Olsen, moisture, SOC, pH and electrical conductivity (EC) significantly (p < 0.05) influenced the bacterial communities (Table S4). Changes in physicochemical parameters have been severally reported to have a strong influence on bacterial community structure in diverse environments [17,84,115,118]. The impact of some of these physicochemical terms and the microbial community changes they drive in rice-cultivated soils have been extensively reviewed [94].

Proteobacteria was the most abundant phyla across most RCST, while others with a relatively high abundance were Bacteroidetes, Acidobacteria and Actinobacteria. This finding is consistent with previous studies that profiled the abundant bacterial phyla in rice-cultivated mudflats [118,119]. An important observation in the distribution of bacteria across all RCST was the dominance of Proteobacteria and Acidobacteria in the rice cultivated PH soil in Chokwe and also the dominance of Bacteroidetes in the rice cultivated soil in Umbeluzi. Proteobacteria comprises several species that are considered to be fast-growing copiotrophs and thus proliferate in nutrient-rich soils [120]. They are also known to be rice-straw degraders [121], and this explains their relatively high abundance in the rice-grown fields in both Umbeluzi and Chokwe. Meanwhile, Bacteroidetes and Acidobacteria are among the predominant phyla that have been reported in the rice agroecosystems [120,122,123]. They are reportedly involved in the degradation of complex polysaccharides, including cellulose, hemicellulose, starch, pectin and xylans [124,125,126]. Their potential role in the rice-cultivated soil is to mobilize nutrients necessary for plant growth through active participation in the degradation of complex organic matter. In contrast to the rice cultivated fields, Actinobacteria and Chloroflexi dominated the fallow (10YF and 1YF) fields in both locations and increased significantly in the 5mR time point. This observation suggests that both phyla play important roles in recovering rice cultivated fields. Bacterial species of the phylum Actinobacteria contribute directly to plant growth and yield through phosphate solubilization, nitrogen fixation and iron acquisition [127]. They also contribute to plant resistance to biotic and abiotic perturbations through the release of biocidal or inhibitory chemical substances [127].

Bacterial composition at the genus taxonomic rank revealed that distinct taxa dominated at different time points. Interestingly, Pseudomonas and MND1 (Nitrosomonadaceae) were biomarkers in the rice-cultivated fields in Umbeluzi and Chokwe, respectively, while Agromyces, Bacillus, Lysobacter, Micromonospora, Norcadiodes, Rubrobacter, Solirubrobacter and Sphingomonas were differentially abundant in the 5mR and fallow (10YF, 5YF and 1YF) fields. This observation suggests that rice cultivation inhibits the growth and potential contribution of these organisms to soil health or that the allowance of a fallow period supports their growth and proliferation. Sphingomonas are rhizosphere-associated bacteria and have been reported to improve plant yield through siderophore production [128,129] and herbicide metabolism in the agroecosystems [130]. Both Bacillus and Lysobacter comprise several species that are established biocontrol agents, supporting plant growth through the release of exopolysaccharides that prevents the growth of plant pathogenic microorganisms [131,132]. Siderophores and polyamines producing Micromonospora have been demonstrated to significantly improve the length of roots, shoots and dry weight of Salicornia bigelovii plants [133]. Furthermore, Micromonospora, Norcadiodes, Rubrobacter and Solirubrobacter are within the phylum actinobacteria. In the agroecosystems, Actinobacteria contribute to nutrient cycling and the degradation of organic compounds, including pesticides and herbicides [127,134,135]. This panoply of established plant growth-promoting bacterial species demonstrates that the 5-months postharvest fallow period is important for both soil nutrient recovery and plant pathogen suppression before the next planting cycle. Overall, our findings indicate that rice cultivation significantly impacts the compositional structure and potential function of microbial communities.

The co-occurrence network enabled the detection of ecological rearrangement that occurs as a result of rice cultivation. Co-occurrence networks have been demonstrated to be an important tool for underpinning microbial community assembly patterns in their natural environments and also in response to environmental stressors or disturbances [136,137,138,139,140]. Similar to real-world networks like human social relationships [141] and the world wide web [142], the networks for both the rice-cultivated and fallow fields followed a power-law distribution. This implies that few bacterial genera and chemical properties have many interactions, while others either have very few or no associations.

The high modularity in all RCST is a direct indication of niche partitioning potential by the co-occurring microbes [143]. The formation of broad ecological niches is reportedly tied to the high diversity of bacterial species involved in several metabolic processes, including the cycling of key nutrients [116]. Also, niche separation reduces competition while promoting the co-evolution of specific bacterial phylotypes within the microbiome [144]. The topological properties used to infer network complexity revealed that the fallow fields had a more complex and stable bacterial community than rice-cultivated fields. This finding is similar to a previous study by Niu et al. [145], where they reported that the complexity of microbial networks decreased with tobacco cultivation while the fallow fields maintained a complex microbial interaction. Microbial community co-occurrence networks with reduced complexity are considered more susceptible to environmental disturbances [146]. Factors such as intensified cropping and frequent fertilizer and herbicide application have been shown to reduce microbial community co-occurrence network complexity [147]. Summarily, the overall network complexity (i.e., the microbial communities’ efficiency) is a function of both rice-cultivation status and site differences. In contrast, the high modularity value (ranging from 0.8–0.9) in all networks (Figure 7a–d) indicates the presence of distinct bacterial ecological niches in both the rice-cultivated (PH) and fallow (5YF and 10YF) fields in Umbeluzi and Chokwe.

Keystone species in the rice-cultivated field in Umbeluzi included Pseudomonas, Lysobacter, Stenotrophobacter, Nitrospira, Sphingomonas and RB41. In the rice-cultivated field in Chokwe, Nitrospira and RB41, among other genera, were also found to be keystone species, while Geobacter, Stenotrophobacter, Sideroxydans, Anaeromyxobacter and Candidatus Nitrososphaera were among module hubs. Network hubs are important because they control the flow of information across the entire network and, therefore, determine the community structure [148]. In this study, most network hubs in the rice-cultivated fields are known to be involved in energy-cycling processes and suppressing plant pathogens. Lysobacter preys on pathogenic microorganisms by releasing a vast array of lytic enzymes and antibiotics, making them good biocontrol agents [149]. They are also known to be in symbiotic relationships with rhizosphere-associated microbial groups [150], where they contribute to plant productivity by suppressing the activity of pathogenic fungi, bacteria, and algae [151]. The genus Stenotrophobacter grows optimally at neutral, slightly acidic, or slightly basic pH conditions and is known to secrete exoenzymes (e.g., alkaline and acid phosphatases) involved in the metabolism of phosphate-containing compounds [152]. Nitrospira is an important nitrogen-fixing bacteria, and some strains have been reported to exhibit complete nitrification [153], while another keystone species, Sphingomonas, produces siderophores and contributes to nitrogen fixation [154]. The detection of these potential plant beneficial microorganisms as network hubs in the rice-cultivated fields supports our hypothesis that the PH and 5mR time points are important soil recovery periods, particularly considering nitrates, total nitrogen, P, and K increased significantly between PH and 5mR time points in Chokwe. In addition, the detection of iron-oxidizing and reducing bacteria in association with carbon and several established sulfur and nitrogen metabolizers suggests the existence of a coupled metabolism of iron, carbon, nitrogen and sulfur in the rice-cultivated field in Chokwe. Already, there is documented evidence about the syntrophic relationships between iron-oxidizing bacteria and iron-reducing bacteria and between iron-oxidizing bacteria and sulfate-reducing bacteria, and how their overall interaction affects biogeochemical cycles in freshwater and marine environments [155,156,157,158,159]. Furthermore, important is the fact that iron-oxidizing bacteria have been previously identified as keystone species that determine the structure of microbial communities [148] irrespective of their abundance in the environment.

5. Conclusions

This study revealed that rice cultivation and differences in soil physicochemical parameters were the most important factors that influence alpha diversity and bacterial assemblages into diverse ecological niches. The compositional assessment revealed that the relative abundance of Actinobacteria significantly increased between the postharvest and 5-months off-season fallow period in both study locations, while the abundance of Bacteroidetes reduced significantly. The 5-month off-season fallow period is also important for soil recovery, considering the panoply of established plant growth-promoting bacterial species that were biomarkers at this time point. The positive change in soil health indicators between the postharvest and 5-months off-season fallow period further established this fact. The bacterial co-occurrence network revealed that the ecological rearrangement of the soil microbiome due to rice cultivation led to a reduction in the efficiency of the bacterial communities. Overall, this study demonstrates that rice cultivation and differences in physicochemical parameters are important factors that determine bacterial alpha diversity, structure and the diversity of their environmental functions.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/agronomy11040694/s1, Figure S1: Rarefaction curve of Amplicon sequence variants sampled to a depth of 11,000 sequences per sample, Figure S2: Shared and unique ASVs of RCST in Umbeluzi and Chokwe. Venn diagram was computed using the DeepVenn online tool (http://www.deepvenn.com/; accessed on 10 January 2021), Table S1: Post hoc-adonis (PERMANOVA) for soil bacterial community structure for RCST, Table S2: Differentially abundant phylotypes in Umbeluzi, Table S3. Differentially abundant phylotypes in Chokwe, Table S4: Goodness of fit of the environmental variable with significant influence on bacterial communities in Chokwe, Table S5: Topological properties for the ecological networks in Figure 6.

Author Contributions

Conceptualization, A.M.P.M., A.I.R.-B. and D.F.; methodology, O.T.E., B.A.M., I.F., M.N.T., A.N. and J.C.; validation, J.C. and R.A.A.; formal analysis, O.T.E., C.C.O. and I.F.; investigation, B.A.M., A.M.P.M., A.I.R.-B. and D.F.; resources, D.F., A.I.R.-B., R.A.A., A.N. and A.M.P.M.; data curation, O.T.E., D.F. and J.C.; writing—original draft preparation, O.T.E., C.C.O., V.N.A.N. and B.A.M.; writing—review and editing, A.M.P.M., A.I.R.-B., R.A.A., A.N., J.C., I.F. and D.F.; visualization, O.T.E. and C.C.O.; supervision, A.M.P.M., A.I.R.-B., R.A.A., A.N., J.C. and D.F.; funding acquisition, A.M.P.M., A.I.R.-B. and D.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Portuguese Foundation for Science and Technology (UID/AGR/04,129/2020 (LEAF); UIDB/00239/2020 (CEF); UIDB/00616/2020 (CQ-VR) and UIDB/04033/2020 (CITAB)), the Camões, Instituto da Cooperação e da Língua and Fundação para a Ciancia e a Tecnologia and the International Rice Research Institute (A-2014-98). Ana M. P. Melo received funds from BioData.pt (LISBOA-01-0145-FEDER-022231) and ELIXIR-CONVERGE (INFRADEV-03-2019/Ref: 871075), funded by Horizon 2020.

Data Availability Statement

Sequence data generated for this study are available in the Sequence Read Archives of the National Centre for Biotechnological Information under the accession number PRJNA511657.

Acknowledgments

The Centre for High-Performance Computing (CHPC), South Africa, is hereby acknowledged for providing computational resources to this research (Project CBBI1183).

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Relationship between enzyme activities and soil chemical properties in Umbeluzi (a–d) and Chokwe (e–i). Rho (ρ) and p-values embedded in the plot are rounded off to 2 d.p. and 3 d.p., respectively. 10YF, ten years fallow; 5YF, five years fallow; 1YF, one-year fallow; PH, postharvest; 5mR, five months off-season time point. “10YF/5YF” denotes fallow fields in Umbeluzi or Chokwe, respectively. SOC, Soil organic carbon.

Figure 2. Observed ASVs, Chao1, Shannon–Wiener index of diversity and Pielou’s evenness in Umbeluzi (a–d) and Chokwe (e–h). Black-filled dots depict mean values on the boxplot. Plots for 10YF, 1YF or PH with different superscript letters for each location are significantly different (Bonferroni adjusted p < 0.05) based on Kruskal–Wallis or Mann–Whitney test. 5mR is only compared to PH. Plots for 5mR with asterisks (*) is significantly (Wilcoxon rank-sum test p < 0.05) different from PH.

Figure 3. Multivariate differentiation of bacterial communities in rice cultivation statuses or time points (RCST). (a,b) Umbeluzi (c,d) Chokwe. Principal coordinate analysis (PCoA) plot and unweighted pair group method with arithmetic mean (UPGMA) cluster dendrogram are based on Bray-Curtis dissimilarities. Dotted lines in the PCoA plot show the distance of every sample to its group centroids in multivariate space, while ellipses show 95% confidence intervals (standard error) in multivariate space around group centroids. 10YF, ten years fallow; 5YF, five years fallow; 1YF, one-year fallow; PH, postharvest; 5mR, five months off-season time point. “10YF/5YF” denotes fallow fields in Umbeluzi or Chokwe, respectively. The differences in multivariate space among RCSTs in Umbeluzi and Chokwe are significant (PERMANOVA p < 0.05; see Table S1).

Figure 4. Dominant and differentially abundant (false discovery rate (FDR)-adjusted p < 0.05, least discriminant analysis (LDA) > 4.0) phylotypes at phylum taxonomic rank. (a) Mean relative abundance (b) differentially abundant phyla in Umbeluzi and (c) Differentially abundant phyla in Chokwe. 10YF, ten years fallow; 5YF, five years fallow; 1YF, one-year fallow; PH, postharvest; 5mR, five months off-season time point.

Figure 5. Relative abundance and differentially abundant phylotypes at the genus taxonomic rank. (a) Mean relative abundance. (b) Top 10 differentially abundant phylotypes (FDR-adjusted p < 0.05, LDA > 4.0) in Umbeluzi. (c) Top 10 differentially abundant phylotypes (FDR-adjusted p < 0.05, LDA > 4.0) in Chokwe. In Figure 4a, phylotypes with relative proportion less than 0.01 (1%) and phylotypes unclassified at the genus taxonomic level were excluded. Family names are in parenthesis for some phylotypes without standard genus names. See Tables S2 and S3 for other differentially abundant phylotypes. 10YF, ten years fallow; 5YF, five years fallow; 1YF, one year fallow; PH, postharvest; 5mR, five months postharvest recovery or pre-planting.

Figure 6. Constrained redundancy analysis (RDA) explaining bacterial species-environment interaction for the different RCST in Umbeluzi (a) and Chokwe (b). 10YF, ten years fallow; 5YF, five years fallow; 1YF, one year fallow; PH, postharvest; 5mR, five months off-season recovery period. “10YF/5YF” denotes fallow fields in Umbeluzi or Chokwe, respectively. Gray dots on the plot represent bacterial species.

Figure 7. Co-occurrence network of bacterial communities in rice-cultivated fields in Umbeluzi (a) and Chokwe (b) and in the fallow fields in Umbeluzi (c) and Chokwe (d). The nodes represent bacterial genera or environmental variables and are colored according to modules. Nodes with <3 connections are colored gray. The size of each node is proportional to their betweenness centrality. The connections (edges) stand for significant (FDR-adjusted p < 0.01) associations, while their sizes are proportional to ρ. Positive associations (edges) are colored green, while negative associations are colored red. Rectangles denote environmental variables, while ellipses denote bacteria phylotypes.

Table 1

Study site description.

Location(Region)	Site Code	Description (Rice Cultivation Status or Time Point, RCST)	Sampling Date	Sampling Time (Planting Season)
Umbeluzi(Maputo)	10YF	At least ten years without rice cultivation	May 2016	Fallow field
	PH	Cultivated fields with Makassane rice variety under an irrigated system. The fields were completely flooded. Samples were collected immediately after rice harvest (i.e., postharvest)	May 2016	Postharvest (2015/2016)
	1YF	One year without growing rice. Cultivated with Makassane rice variety one year prior	May 2016	Fallow field
	5mR	Same site as PH. Samples were collected five months postharvest but before pre-planting preparation for the next season. The field was completely flooded at the time of sampling	October 2016	Pre-planting (2016/2017)
Chokwe(Gaza)	5YF	At least five years without rice cultivation. The fields were completely flooded at the time of sampling	May 2016	Fallow field
	PH	Cultivated fields with Makassane rice variety in an irrigated system. The fields were completely flooded. Samples were collected immediately after rice harvest (i.e., postharvest)	May 2016	Postharvest (2015/2016)
	5mR	Same site as PH. Samples were collected five months postharvest but before pre-planting preparation for the next season. The field was completely dry at the time of sampling	October 2016	Pre-planting (2016/2017)

Table 2

Selected physicochemical properties of bulk soil.

Properties	Umbeluzi				Chokwe
	10YF	1YF	PH	5mR	5YF	PH	5mR
pH (H2O)	7.17 ± 0.93 b	7.67 ± 0.48 b	8.25 ± 0.21 a	7.99 ± 0.49	7.76 ± 0.50 a	8.08 ± 0.29 a	7.11 ± 0.69 *
EC (mS/cm)	5.18 ± 6.35 a	0.41± 0.05 b	0.55 ± 0.13 b	4.45 ± 0.66 *	1.00 ± 0.32 a	0.63 ± 0.31 a	5.42 ± 3.06 *
Moisture (g H2O kg−1)	30.50 ± 9.18 a	27.84 ± 8.33 a	29.71 ± 12.18 a	22.75 ± 21.12	53.35 ± 2.31 a	56.67 ± 5.82 a	39.47 ± 11.32 *
SOC (g kg−1)	16.09 ± 5.08 a	14.00 ± 2.25 a	15.24 ± 1.85 a	15.92 ± 1.44	20.68 ± 4.16 a	14.44 ± 1.21 b	16.23 ± 0.91 *
NH4+–N (mg N kg−1)	33.28 ± 58.39 a	2.13 ± 0.83 c	6.96 ± 2.34 b	4.73 ± 0.90 *	6.36 ± 1.38 a	4.34 ± 0.93 b	10.22 ± 15.1
NO3− –N (mg N kg−1)	29.17 ± 51.46 a	0.10 ± 0.31 b	0.06 ± 0.19 b	1.87 ± 4.62 *	0.00 ± 0.00	0.00 ± 0.00	22.52 ± 25.12 *
Total N (g kg−1 dry soil)	0.70 ± 0.34 a	0.75 ± 0.19 a	0.77 ± 0.15 a	0.91 ± 0.12	1.34 ± 0.63 a	0.77 ± 0.11 a	1.14 ± 0.08 *
P-Olsen (mg kg−1)	28.30 ± 8.86 a	26.44 ± 6.59 a	29.87 ± 8.22 a	20.81 ± 1.65	15.42 ± 3.91 a	9.67 ± 3.58 a	9.39 ± 1.67
Egner Rhiem K (mg kg−1)	240.70 ± 45.21 a	251.37 ± 35.77 a	239.25 ± 33.92 a	328.59 ± 60.59 *	463.61 ± 104.05 a	524.98 ± 60.14 a	492.85 ± 249.28

Values are mean ± SD (N = 9). Mean for 10YF, 5YF, 1YF or PH with different superscript letters across rows for each location (i.e., within location comparison) are significantly different (Bonferroni adjusted p < 0.05) based on Kruskal–Wallis or Mann–Whitney test. 5mR is compared to the only PH. * 5mR is significantly (p < 0.05) different from PH based on the Wilcoxon rank-sum test. Comparisons between locations were not performed due to differences in the number of sampled sites, soil type and varying agronomic practices. EC, electrical conductivity; SOC, soil organic carbon.

Table 3

Enzyme activities in soils.

Location	RCST	β-Glucosidase(µg p- Nitrophenol g−1 h−1)	Acid Phosphatase(µg p- Nitrophenol g−1 h−1)	Alkaline Phosphatase(µg p- Nitrophenol g−1 h−1)	Urease(mg N–NH4+ kg−1 2 h−1)
Umbeluzi	10YF	41.18 ± 25.55 ab	303.54 ± 201.14 a	218.68 ± 104.83 a	49.59 ± 32.97 a
	1YF	43.99 ± 14.5 a	188.31 ± 57.79 a	170.91 ± 31.69 a	42.67 ± 8.97 a
	PH	25.83 ± 3.68 b	127.15 ± 31.99 a	148.68 ± 39.31 a	25.11 ± 11.33 a
	5mR	30.50 ± 6.69 *	122.98 ± 38.36	121.43 ± 45.49	24.34 ± 11.64
Chokwe	5YF	60.70 ± 23.99 a	351.9 ± 111.1 a	240.51 ± 69.41 a	71.22 ± 19.47 a
	PH	58.51 ± 9.97 a	198.37 ± 44.24 b	172.99 ± 43.16 a	82.00 ± 21.10 a
	5mR	74.26 ± 12.91 *	220.03 ± 27.04	191.61 ± 27.48	82.30 ± 20.67

Values are mean ± SD of 9 bulk samples. Means for 10YF, 1YF or PH with different superscript letters across columns for each location (Umbeluzi or Chokwe) are statistically different (p < 0.05). * Significantly (p < 0.05) different from PH-based on the Wilcoxon rank-sum test (paired-sample t-test statistic) for each location. Comparisons between locations were not performed due to differences in the number of sampled sites, soil type and varying agronomic practices.