Influence of organic, synthetic and biofertilizers on the diversity of cassava rhizosphere microbiome in Northeastern Thailand

Experimental sites

Cassava was grown at two sites in Northeastern Thailand, including Nampong district, Khon Kaen Province and Seungsang district, Nakhon Ratchasima Province, both of which are owned by smallholder farmers who granted permission to sample soil on their properties. The Northeastern region of Thailand is the current dominant cassava growing region of the country. The soil type at the Nampong site (Latitude: 16.823750, Longitude: 102.979556) is a Nam Phong soil series (loamy, siliceous, isohyperthermic Grossarenic Haplustalfs) and the soil at the experimental field has a sandy texture (Kunlanit, Khwanchum & Vityakon, 2020). The soil type at the Seungsang site (Latitude: 14.490722, Longitude: 102.508167) is a Chokchai soil series (very fine, kaolinitic, isohyperthermic Rhodic Kandiustox) and the experimental field soil has a silty clay texture (Sirichuaychu et al., 2005). The soil chemical properties and growth season temperatures of the two experimental site regions were explained in our previous study (Thompson et al., 2025).

Experimental treatments

The investigation of changes in the rhizosphere microbiome under various fertilizer applications was conducted across eight experimental treatments, using the commercial Rayong 72 cassava cultivar, which is well-known for its high yield and resistance to drought (Sarakarn et al., 2000). The experiment was designed in a randomized complete block design (RCBD) with five replicates. Eight experimental treatments with different fertilizer management protocols consisted of (T1) no fertilizer application (control), (T2) recommended dosages of synthetic fertilizer based on soil analysis (RDCF) (Department of Agriculture, 2005), (T3) 3.12 t ha⁻¹ of chicken manure (CM), (T4) 937.5 L ha⁻¹ of swine manure extract (SME), (T5) 3.12 t ha⁻¹ of chicken manure + 937.5 L ha⁻¹ of swine manure extract (CM + SME), (T6) 50% of the recommended dose of synthetic fertilizers + 1.56 t ha⁻¹ of chicken manure (1/2 RDCF + 1/2 CM), (T7) 50% of the recommended dose of synthetic fertilizers + 468.8 L ha⁻¹ of swine manure extract (1/2 RDCF + 1/2 SME), and (T8) 3.12 ton ha⁻¹ of chicken manure + stalk soaking with plant growth-promoting rhizobacteria (PGPR-3) solution containing two types of bacteria: Azospirillum brasilense and Gluconacetobacter diazotrophicus (CM + PGPR). Details of application rate and macronutrient contents for each fertilizer application were explained in our previous study (Thompson et al., 2025).

Plant materials, DNA extraction and 16S rRNA amplicon sequencing

Soil samples from cassava tubers were collected during the growing season at three stages, including 2, 5 and 10 months after planting (MAP) at both Nampong and Seungsang field sites. Soil samples were collected from around the tubers of plants grown with all eight experimental treatments. For each stage, three cassava plants (replicates) were sampled per treatment, resulting in a total of 24 samples per site for each stage (48 samples for both sites at each cassava stage). A total of 144 samples were collected, from three stages of growth (2, 5, and 10 MAP) under the eight nutrient management treatments, at both study sites. For each sample, soil surrounding the tubers was collected by gently shaking and manually removing it with gloved hands. The soil of each sample was collected in plastic bags and was thoroughly mixed by repeatedly flipping the bags. For DNA extraction, 250 mg of soil from each sample was used per tube. The soil samples were stored at −20 °C before extraction. DNA extraction was carried out using the DNeasy PowerSoil Pro kit (Qiagen, Hilden, Germany). DNA quality and concentration were assessed using both 1% agarose gel and a Nanodrop ND-1000 Spectrophotometer. The extracted DNA was stored at −20 °C before proceeding to the sequencing step. For amplicon metagenomic sequencing, 200 ng of total DNA and more than 10 ng/µl of the DNA concentration was needed. The 16S rRNA amplicon amplification targeted the V3–V4 regions, with a fragment length of 470 bp. The primers used to amplify these regions were primer 341F, 5′CCTAYGGGRBGCASCAG 3′ and primer 806R, 5′ GGACTACNNGGGTATCTAAT 3′. The amplicon sequencing was conducted by an Illumina NovaSeq 6000 platform by Novogene Corporation (Singapore) using 2 × 250 bp pair-end reads. The output of total 16S rRNA amplicon sequences was used from microbiome analysis in the next step.

Microbiome analysis

After the Illumina sequencing step, all reads of the 16S rRNA amplicon sequences from all 144 soil samples were prepared for the microbiome analysis. Quality control for raw paired-end sequences was performed by trimming adapter and primer sequences and filtering low quality sequences (Phred Quality Score < 30) using the FASTP program (Chen et al., 2018). Cleaned paired-end sequences were subsequently joined into single-end sequences using FLASH version V1.2.7 (Magoč & Salzberg, 2011). Data analysis of 16S rRNA sequences was performed using QIIME 2 (version 2022.02) (Bolyen et al., 2019). The 16S rRNA sequences were quality filtered, denoised, and amplicon sequence variants (ASVs) were constructed using DADA2 (version 1.14) (Callahan et al., 2016) implemented in QIIME2. The minimum total frequency and number of observed sample parameters were defined as 100 and 2, respectively. Taxonomic classification was performed using the BLAST+ consensus taxonomy classifier against the SILVA database (released version 138.1) (Quast et al., 2013) by specifying the following parameters: percent identity cutoff = 0.9 and minimum fraction of consensus assignment = 0.5. Any sequences containing chloroplast or mitochondria according to taxonomic classification were removed before moving to the next step. Beta diversity analysis based on Bray-Curtis dissimilarity matrix was used to measure the similarity among samples from different treatments and was visualized by principal coordinate analysis (PCoA). The significant effects of treatments in shaping bacterial communities were calculated by permutational multivariate analysis of variance (PERMANOVA) using 999 as the number of permutations. The raw sequence data were deposited under NCBI BioProject Accession number PRJNA1202398. Due to the limitations of using only 16S rRNA amplicon data, which does not capture functional microbial potential, FAPROTAX (Louca, Parfrey & Doebeli, 2016) analysis was perform to translate our bacterial taxonomic profiles into putative metabolic functions, allowing us to infer the ecological roles of bacterial community across treatments.

Analysis of predominant bacteria and correlation coefficient with cassava production

The predominant bacteria were identified by comparing genus-level bacteria that showed significant differences across experimental treatments at the three stages (2, 5, and 10 MAP). Emphasis was placed on experimental treatments that contributed to the highest growth and yield of cassava at the Nampong and Seungsang sites. Comparisons were made using White’s non-parametric t-test in STAMP (Parks et al., 2014) with a significance level of p < 0.05. White’s non-parametric t-test was selected for its robustness in detecting differences in bacterial abundance in microbiome data, accommodating the non-normal distributions and small sample sizes. To identify the predominant bacteria at genus-level correlated with fresh tuber yield and starch content recorded at 10 MAP, the Spearman correlation coefficient was used with Scipy python package (Virtanen et al., 2020). Details of the fresh tuber yield and starch content for the experimental treatments are reported in a previous study (Thompson et al., 2025).

DNA extraction and 16S rRNA sequencing

DNA from a total of 144 soil samples, from the eight treatments at the Nampong and Seungsang sites, collected at three cassava stages (2, 5, and 10 MAP), was extracted and used for 16S rRNA amplicon sequencing using Illumina sequencing technology. The quality of the genomic DNA was found to be high and sufficient for sequencing. The genomic DNA bands from all soil samples were larger than 5,000 bp, with the 16S rRNA fragment approximately 1,500 bp in size. The DNA concentration in nearly all soil samples exceeded 50 ng/µl. The cassava plants, along with their roots and tubers with surrounding soil, collected at 2, 5, and 10 MAP from the Nampong and Seungsang sites, is shown in Fig. 1, along with the DNA patterns of soil samples.

Amplicon sequence variant derived from 16S rRNA sequences

To study the cassava rhizophere bacterial communities resulting from different fertilizer treatments, 16S rRNA amplicon sequencing was performed using the Illumina NovaSeq 6000 platform with paired-end 250 bp reads. The amplicons targeted the V3–V4 variable regions, with a fragment size of 470 bp. After sequencing DNA from the 144 soil samples, the initial analysis using MultiQC showed that the read lengths of all 144 samples ranged from 386 to 430 bp. Amplicon sequence variant (ASV) data for cassava rhizophere soil at 2, 5, and 10 MAP at the Nampong site is shown in Table S1. At 2 MAP, the non-chimeric ASV count ranged from 84,533 to 113,786, accounting for 43.76% to 60.27% of the input 16S rRNA sequences. At 5 MAP, the non-chimeric ASV count ranged from 89,748–156,981, accounting for 69.16–84.02% of the input 16S rRNA sequences. At 10 MAP, the non-chimeric ASV count ranged from 69,497–161,325, accounting for 70.99–87.03% of the total input 16S rRNA sequences. ASV data for cassava rhizophere soil at 2, 5, and 10 MAP at the Seungsang site is shown in Table S2. At 2 MAP, the non-chimeric ASV count ranges from 78,003–100,871, accounting for 42.42–51.7% of the total input 16S rRNA sequences. At 5 MAP, non-chimeric ASV count ranged from 91,829–143,529, accounting for 73.92–84.49% of the total input 16S rRNA sequences. At 10 MAP, non-chimeric ASV count ranged from 54,038–145,551, accounting for 41.75–85.68% of the total the input 16S rRNA sequences.

Taxonomic classification at phylum and genus levels

The microbial composition at the phylum level for all cassava rhizophere soil samples at 2, 5, and 10 MAP from the Nampong and Seungsang sites showed similar proportions, as illustrated in Fig. 2A. It was found that the microbial populations at the phylum level at 2 MAP were noticeably different from those at 5 and 10 MAP, with the populations at 5 and 10 MAP showing only slight differences. At 2 MAP, the Firmicutes phylum was the most abundant, followed by Proteobacteria and Actinobacteria. At 5 and 10 MAP, Actinobacteria was the most abundant phylum, followed by Proteobacteria and Firmicutes.

The analysis of microbial proportions at the genus level for all cassava rhizophere samples at 2, 5, and 10 MAP showed similar results at both Nampong and Seungsang sites, as illustrated in Fig. 2B. The microbial populations at the genus level at 2 MAP differed from those at 5 and 10 MAP, mirroring the differences observed at the phylum level. The microbial populations at the genus level were similar at both 5 and 10 MAP. At 2 MAP, the Bacillus genus was the most abundant. However, at 5 and 10 MAP, Bacillus significantly decreased, though it remained the most abundant genus compared to others, including Steptomyces, Mycobacterium, Conexibacter and Acidothermus.

Alpha diversity of microorganisms

The alpha diversity indicated that overall observed microbial features in the soil around the tubers at 2, 5, and 10 MAP were similar at both the Nampong and Seungsang sites (Table S3, Fig. 3). The highest abundance of observed microbial features occurred at 5 MAP, with a decline at 10 MAP. The number of observed features at 5 and 10 MAP was significantly higher than that observed at 2 MAP at both sites. At Nampong, observed microbial features decreased at 10 MAP, but no significant difference was found between 5 and 10 MAP. In contrast, at Seungsang, observed microbial features decreased at 10 MAP and showed a significant difference compared to 5 MAP.

The observed bacterial features around the tubers at both sites for each treatment at 2, 5, and 10 MAP are illustrated on Fig. 4 and the details from treatments that resulted in significant differences are shown in Table S4. At Nampong, no significant differences in observed microbial features were found across experiments at 2 MAP. However, by 5 MAP, the no-fertilizer control (T1) had fewer observed microbial features than treatments that included chicken manure, such as CM (T3), CM + SME (T5), 1/2 RDCF + 1/2 CM (T6), and CM + PGPR (T8), while no significant differences were observed among the control (T1), chemical fertilizer (T2), and SME (T4) treatments. At 10 MAP, the alpha diversity of CM (T3) remained significantly higher than that of the control (T1) and RDCF (T2).

At Seungsang, significant differences among treatments were observed as early as 2 MAP, with SME (T4) and CM + SME (T5) treatments showing higher bacterial diversity than the no-fertilizer control (T1) and chemical fertilizer RDCF (T2) treatments. At 5 MAP, CM (T3) and 1/2 RDCF + 1/2 CM (T6) treatments exhibited greater alpha diversity than did SME (T4), 1/2 RDCF + 1/2 SME (T7), and CM + PGPR (T8) treatments. At 10 MAP, significant differences were also observed, with the control (T1) and CM + PGPR (T8) treatments displaying higher alpha diversity than what resulted from the SME (T4) treatment.

Beta diversity of microorganisms

The beta diversity of microorganisms in the soil surrounding cassava roots at 2, 5, and 10 MAP in Nampong and Seungsang is illustrated in Fig. 5. The results indicate significant differences in beta diversity across the different growth stages within each site. Additionally, the beta diversity between the two sites also showed slight variations. At Nampong, the beta diversity at 2 MAP differed significantly from that at 5 and 10 MAP, while the diversity at 5 and 10 MAP was nearly the same. Similarly, at Seungsang, there was some overlap in beta diversity across 2, 5, and 10 MAP, but the diversity at 5 and 10 MAP was almost identical, as observed in Nampong. However, within each treatment, the beta diversity at 2, 5, and 10 MAP did not show significant differences for both sites (Table S5).

Analysis of metabolic functions of microorganisms and predominant bacteria resulting from different fertilizer applications

To compare bacterial genera resulting from different treatments, bacteria were classified from soil samples around cassava tubers grown in Nampong and Seungsang at 2, 5, and 10 MAP. At 2 MAP, 95 bacterial genera were identified, accounting for approximately 66.25–96.15% of the total bacteria found in the 48 soil samples collected from both sites. At 5 MAP, 92 bacterial genera were identified, accounting for approximately 56.86–77.34% of the total bacteria found in the 48 soil samples from both sites. At 10 MAP, 66 bacterial genera were identified, accounting for approximately 63.41–81.81% of the total bacteria found in the 48 soil samples from both sites. Details of bacterial genera are provided in Table S5. To infer putative metabolic functions from bacterial taxonomic profiles, we employed FAPROTAX to assess the ecological roles of bacterial communities across treatments, as presented in Fig. S1 and Table S6. The dominant predicted metabolic functions were chemoheterotrophy and aerobic chemoheterotrophy, observed at 2, 5, and 10 MAP. These functions peaked at 10 MAP at both the Nampong and Seungsang sites. However, there were no clear differences in the predicted metabolic functions among the different treatments.

Previous work (Thompson et al., 2025) has shown that CM (T3), CM + SME (T5), and CM + PGPR (T8) treatments had a significant impact on cassava growth and yield. Specifically, T3 and T5 had a pronounced effect at Seungsang, while T8 notably influenced cassava growth and yield at Nampong. Therefore, the predominant bacteria in treatments T3, T5, and T8 were considered to be target bacteria. In summary, these target bacteria, including Pseudomonas, Tumebacillus, Lysinibacillus, Paenibacillus, Dongia, Acidibacter, Sphingomonas, and Bacillus are likely to contribute to cassava productivity at 2, 5, and 10 MAP. A comparison of dominant bacteria, analyzed using White’s non-parametric t-test in STAMP, is presented in Table 1 and Fig. 6.

Correlation of predominant bacteria with cassava production

Correlation analysis was conducted to evaluate the relationship between target bacteria and cassava fresh tuber yield as well as between target bacteria and starch content at harvesting time, 10 MAP. Of the eight targeted genera, only six were identified at 10 MAP, including Tumebacillus, Bacillus, Paenibacillus, Lysinibacillus, Dongia, and Sphingomonas. Spearman correlation analysis indicated that Tumebacillus was the most prominent genera, showing a significant correlation with fresh tuber yield (p = 0.041 at Seungsang and p = 0.060 at Nampong) (Fig. 7, Table 2). Paenibacillus and Bacillus also exhibited a notable relationship with yield (p = 0.054 at Seungsang and p = 0.062 at Nampong) (Fig. 7, Table 2). However, Spearman correlation analysis found no significant relationship between the six target bacteria and starch content (Table S8).