The role of microbiomes in cooperative detoxification mechanisms of arsenate reduction and arsenic methylation in surface agricultural soil

Sampling site characteristics and sample collection

The sampling site was an agricultural area located in Samphran District, Nakhon Pathom Province, Thailand. The land is mainly used for banana cultivation. Chemical fertilizers and pesticides are applied once per month. The study area is surrounded by factories. According to Thai regulations, the arsenic concentration in soil within residential areas and agricultural areas should not exceed 6 and 25 mg/kg, respectively. Since the study area is located in a residential zone and has the potential for arsenic contamination due to anthropogenic activities, an investigation into arsenic genes associated with the soil microbiome was conducted. Surface soil samples (0–15 cm depth) were collected from three plots across the study area. For each plot, approximately 100 g of surface soil were randomly collected from three locations. In order to make a composite sample, the collected soil samples from three locations from each plot were subsequently pooled and mixed on site. Each composite sample was subsequently analyzed for soil properties. The samples were collected in plastic bags and kept on ice during transportation. The soil samples were collected in February (T1) and August (T2) 2021, respectively representing soils during both the dry and wet seasons.

Analysis of soil properties

Soil samples were also collected for the analysis of soil properties which were conducted by Central Laboratory (Thailand) Co., Ltd., according to standard protocols. Soil pH and moisture were analyzed according to the method of analysis of chemical fertilizer B.E. 2559 method 1.02.01 and method 1.04.01, respectively. Total organic carbon (TOC) was measured according to the manual on organic fertilizer analysis, APSRDO. Total nitrogen (TN) and total phosphorus (TP) were analyzed by an in-house method based on AOAC (2019). Arsenic (As) and cadmium (Cd) concentrations were respectively analyzed by inductively coupled plasma optical emission spectroscopy (ICP-OES) and inductively coupled plasma mass spectrometry (ICP-MS), according to an in-house method based on AOAC (2019). Copper (Cu) and zinc (Zn) concentrations were analyzed by ICP-OES, according to an in-house method based on official methods of analysis of fertilizers, Japan (1987). The concentrations of arsenic, cadmium, copper, and zinc were analyzed due to their common prevalence in polluted soils and their toxicity to impacted living organisms (Gong et al., 2020). Soil texture was analyzed by mechanical analysis through the pipette method (Olmstead, Alexander & Middleton, 1930).

DNA extraction

Total genomic DNA was extracted using a DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany), following the manufacturer’s protocols. For each sample, multiple soil samples were separately extracted. The quantity and quality of extracted DNA were estimated by a NanoDrop™ One spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis. The extracted DNA from each sample were subsequently pooled and used for downstream processes.

Analysis of 16S rRNA gene sequencing

Library construction and sequencing of the V3–V4 hypervariable regions of the bacteria 16S rRNA gene were prepared using primers 341F and 806R. Triplicate samples were conducted for each sample (T1_1, T1_2, T1_3 and T2_1, T2_2, T2_3). The 16S rRNA gene libraries and sequencing were prepared by Novogene Co., Ltd. (Beijing, China), using an Illumina NovaSeq 6000 platform, following the manufacturer’s protocol. The Illumina platform generated 250 bp paired-end raw reads. Raw data of 16S rRNA gene amplicon sequences were deposited in the Sequence Read Archive (SRA) under the Bioproject number PRJNA980983.

For the data preprocessing step, the characteristics of the raw sequences (length, quantity, and quality score) were evaluated using FastQC version 0.11.8 (Andrews, 2010). Adapter sequences, including barcode and primer sequences, were subsequently trimmed using Cutadapt version 3.5 (Martin, 2011). The minimum length of the trimmed sequences was 200 bp. Amplicon sequence variants (ASVs) were identified using the Divisive Amplicon Denoising Algorithm 2 (DADA2) (Callahan et al., 2016). All analyzed sequences were truncated at 223 bp for forward reads and 219 bp for reverse reads. A classifier was generated for taxonomic assignment with minimum and maximum length of 380 bp and 470 bp, respectively. Taxonomic information was assigned to all ASVs based on the SILVA 138 SSU database (Quast et al., 2013). ASVs assigned as neither Bacteria nor Archaea were removed. Before the analyzed ASVs were normalized, ASVs with doubletons sequences were eliminated. Alpha rarefaction curves were conducted and alpha diversity indices (e.g., Chao1 and Simpson’s index) were calculated. All the 16S rRNA gene analyses were conducted using the Quantitative Insights Into Microbial Ecology2 (QIIME2) software version 2021.8 (Bolyen et al., 2019). To determine significant difference of microbial taxa between the dry and wet seasons, a two-tailed test was conducted (p-value <0.05). Principal Coordinates Analysis (PCoA) based on Bray–Curtis dissimilarities were generated by QIIME2 version 2021.8 and visualized by the RStudio program. To examine the significant difference between microbial groups, permutational multivariate analysis of variance (PERMANOVA) was conducted using QIIME2 version 2021.8.

Analysis of arsenic specific gene sequencing

Library construction and sequencing of the arsenic functional genes were prepared using specific primers, targeting the aioA, arrA, arsC, and arsM genes (Table S1). The libraries and sequencing were prepared by BTSeq™ Contiguous Sequencing Service, using an Enzymatic Preparation Kit (EP Kit) (Celemics, Korea), following the manufacturer’s protocol. An Illumina MiSeq platform was used to generate 150 bp paired-end raw reads. Raw data of arsenic specific gene amplicon sequences were deposited in the Sequence Read Archive (SRA) under the Bioproject number PRJNA980983.

The qualities of raw sequences (length, quantity, and quality score) were assessed using FastQC version 0.11.8 (Andrews, 2010). Pair-ended reads were joined using the Mothur software (Schloss et al., 2009). The joined sequences were subsequently compared against the previously collected arsenic functional gene sequences from the GenBank database using the blastx tool (Camacho et al., 2009). The merged sequences were evaluated based on the criteria of 80% coverages and 80% similarities. For each gene, the analyzed sequences were clustered into 97% operational taxonomic units (OTUs) (Zhang et al., 2015; Yang et al., 2020) using the CD-HIT program version 4.8.1 (Li & Godzik, 2006). Representative OTU sequences of each gene were selected and verified using the blastn tool (Camacho et al., 2009). The verified OTU sequences were included in phylogenetic analysis. Neighbor-joining trees were constructed using MUSCLE alignment (Edgar, 2004; Tamura, Nei & Kumar, 2004) and were generated with the maximum composite likelihood method with 1,000 bootstrap values, using the MEGA X software (Kumar et al., 2018).

Analysis of shotgun metagenomic sequencing

Library construction and sequencing of the DNA samples (T1_1 to T1_3 and T2_1 to T2_3) were undertaken by Novogene Co., Ltd. (Beijing, China), using an Illumina NovaSeq6000 platform, following the manufacturer’s protocol. The Illumina platform generated 150 bp paired-end raw reads. Raw metagenomic sequencing data were deposited in the Sequence Read Archive (SRA) under the Bioproject number PRJNA980983.

The length and quality scores of raw reads were assessed using FastQC version 0.11.8 (Andrews, 2010). After removing adapter sequences, the analyzed sequences were trimmed using Trimmomatic version 0.36 (Bolger, Lohse & Usadel, 2014). To obtain high-quality sequences, trimming was performed with a quality score cutoff of Phred score >= 20 and a minimum length of 100 bp. Obtained clean reads were assembled into contigs using the metaSPAdes software version 3.15.4 (Nurk et al., 2017). The assembly was performed using multiple k-mer lengths (21, 33, 55, and 77), and the k-mer length of 77 was selected for further analysis. The quality of the assembled contigs was assessed using the QUAST tool version 5.2.0 (Gurevich et al., 2013). The qualified contigs were then annotated using Prokka software version 1.14.6 (Seemann, 2014). The taxonomic profiles were created by mapping metagenomic reads against the standard Kraken2 database. Normalization was employed by scaling with the smallest number of total reads. The qualified contigs were aligned against the customized aioA, arsC, arrA, and arsM database. To construct a customized database of arsenic-related genes, complete protein sequences of bacteria for the aioA, arsC, arrA, and arsM genes were collected from the National Center for Biotechnology Information (NCBI). Subsequently, redundant bacterial sequences were discarded, retaining only the unique ones. To count gene abundance, the annotated contigs were subsequently mapped against the clean reads obtained from data preprocessing using the Bowtie 2 algorithm version 2.2.1 (Langmead & Salzberg, 2012) and the samtools software version 1.16.1 (Li et al., 2009). The abundance of arsenic functional genes was retrieved using the BEDTools software version 2.30.0 (Quinlan & Hall, 2010). The abundance of arsenic functional genes was then normalized using the Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) method (Trapnell et al., 2010). The abundances of each gene across two seasons were compared using a pair t-test.

Characteristics of agricultural soils

The agricultural soil samples were collected from both the dry (T1) and wet (T2) seasons which were 6 months apart. Although the proportions of sand, silt, and clay varied, the soil texture of both T1 and T2 samples was classified as clay (Table 1). Soil pH suggested that the soil remained neutral throughout the year. The concentrations of arsenic, cadmium, and copper in both soil samples were comparable, while the concentration of zinc was much higher in T1 than T2 (Table 1). Total organic carbon (TOC) and total phosphorus (TP) in both soils were detected at low concentrations, while total nitrogen (TN) was not detected (Table 1). Overall, the soil properties during both seasons were similar, except for the apparently elevated zinc concentration in the dry season.

Soil microbiomes of the dry and wet seasons

The 16S rRNA gene sequencing generated 152,505 to 189,101 raw reads per sample, and they were normalized to obtain a total of 82,109 high-quality reads for downstream analysis (Table S2). Rarefaction curves indicated that sequencing depth was optimal (Fig. S1). The numbers of ASVs in T1 and T2 samples were 2,000 to 2,100 and 1,600 to 1,900, respectively. Species richness and evenness of the dry and wet samples were comparable (Table S3). Principal coordinate analysis (PCoA) based on Bray–Curtis distance and PERMANOVA suggested that although soil microbiomes of each season were separately clustered (Fig. S2), they were not significantly different (p-value > 0.05). Overall, the soil microbiome was consistent throughout the year.

Although Archaea were detected in T1 (∼3%) and T2 (∼7%) samples, Bacteria were the dominant domain found in both T1 (∼97%) and T2 (∼93%). The major bacterial phyla found in T1 and T2 samples were Actinobacteriota, Proteobacteria, Firmicutes, Chloroflexi, and Acidobacteriota (Fig. 1A). Other detected phyla found in T1 and T2 soils were Myxococcota, Gemmatimonadota, Nitrospirota, Crenarchaeota, Bacteroidota, and Cyanobacteria. Although microbial profiles and their relative abundance at the phylum level were similar in both soil microbiomes, statistical analysis indicated that the abundance of Actinobacteriota, Proteobacteria, Firmicutes, Acidobacteriota, Myxococcota, Nitrospirota, Crenarchaeota, and Cyanobacteria in T1 and T2 were significantly different (p-value < 0.05) (Fig. 1A).

For better insight into the dynamics of agricultural soil microbiomes, microbial taxa and their abundance were analyzed at the genus level (Fig. 1B). A total of 27 genera present at greater than 1% of the total microbial abundance in at least one samples were identified. Due to the limitation of the database to identify microbial taxa at the genus level, those in which the genus remained unclassified were identified at the lowest taxonomic rank. Common bacterial taxa highly detected in both T1 and T2 soils were Bacillus (∼11%), uncultured_Roseiflexaceae (∼6%), and Streptomyces (∼4%) (Fig. 1B). Microvirga (∼2%) was also consistently present in soils of both the dry and wet seasons. The relative abundance of these taxa in soils from both seasons was not significantly different. Interestingly, the abundance of members of the archaeal phylum Crenarchaeota, including Candidatus Nitrocosmicus, Nitrososphaeraceae, and Candidatus Nitrososphaera, was significantly higher in T2 than in T1 (p-value < 0.05) (Fig. 1B).

Taxonomic profiles retrieved from metagenomics analysis revealed that top five phyla found in soils from both seasons were Actinobacteria, Proteobacteria, Firmicutes, Planctomycetes, and Bacteroidetes (Fig. S3). Among the detected genera, Streptomyces were the most abundant, followed by Nocardioides, Micromonospora, Bacillus, and Pseudomonas, respectively (Fig. S4).

Communities of arsenic-transforming bacteria

The arsenic-functional genes, aioA, arsC, arrA, and arsM, were investigated to study the communities of arsenite-oxidizing bacteria, detoxification arsenate-reducing bacteria, dissimilatory arsenate-reducing bacteria, and arsenic methylation bacteria, respectively. PCR screening of the aioA, arsC, arrA, and arsM genes showed that the arsC gene was undetectable. Instead of common culture-independent techniques, such as clone library construction, high-throughput sequencing of the aioA, arrA, and arsM genes was employed to explore the communities of arsenic-transforming bacteria. High-throughput sequencing of the aioA, arrA, and arsM genes produced 54,654 to 138,093 raw reads per sample (Table S4).

The phylogenetic tree of the aioA sequences revealed that the majority of the aioA sequences retrieved from both soil samples were closely related to uncultured aioA clones previously recovered from diverse environments, such as arsenic-contaminated soil, groundwater, and aquatic sediment (Fig. 2). While the main retrieved aioA sequences from T1 were closely related to uncultured aioA clones, a minor portion of those from T2 was closely related to aioA sequence belonging to Microvirga sp. (Fig. 2). Most analyzed arrA sequences from both soil samples were closely associated with uncultured arrA clones previously discovered in paddy soils, Mekong Delta sediments, a cache valley basin, and arsenic contaminated sediment (Fig. 3). None of the arrA sequences from T2 was associated with known dissimilatory arsenate-reducing bacteria, while those from T1 were closely related to the arrA gene belonging to Geotalea uraniireducens and Desulfuromonas sp. WB3 (Fig. 3). Phylogenetic analysis also demonstrated that the majority of the arsM sequences retrieved from both soil samples were closely associated with uncultured arsM clones previously reported (Fig. 4). The analyzed arsM sequences from T1 were closely related to those of uncultured clones found in estuary sediments, rhizosphere soils, paddy soils, and river water affected by acid mine drainage, while the analyzed arsM sequences from T2 were closely associated with those of uncultured clones previously retrieved from estuary sediments, rice rhizosphere soil, and paddy soil. The only known arsenic methylation bacteria closely related to the arsM sequences from both soil samples were Rhodopseudomonas palustris, which were isolated from soil near a tin mine (Fig. 4).

Overall, the results suggested that the communities of arsenic-transforming bacteria in both soil samples were broadly comparable. The analysis of arsenic-transforming bacterial communities in agricultural soils of both seasons demonstrated that the majority of arsenite-oxidizing bacteria, dissimilatory arsenate-reducing bacteria, and arsenic methylation bacteria found in this study was closely associated with uncultured arsenic-transforming bacteria recovered from various environments.

Prediction of arsenic-transforming pathways

Based on the resulting 16S rRNA gene sequences, arsenic-functional genes were predicted by the PICRUSt2 software against the KEGG database. Although arsenic functional genes were found in both T1 and T2, they accounted for only ∼0.3% of the total predicted genes. To compare the relative abundance of arsenic functional genes across all samples, only the detected arsenic functional genes were extracted and normalized. The results showed that the proportion of arsenic functional genes in the soils of both dry and wet seasons was relatively consistent (Fig. 5). A paired t-test revealed no significant difference (p-value >0.05) in the abundances of each arsenic functional gene between T1 and T2, suggesting comparable potential arsenic-transforming pathways in soil across both seasons. A large proportion of the ars operon (e.g., arsC, ARSC1, ARSC2, ACR3, arsB, and arsA), responsible for arsenic detoxification mechanisms, was found in the soils analyzed in this study. As for arsenic-functional genes responsible for energy metabolism, the aioA and aioB genes were detected at very low abundance in both T1 and T2 soils (Fig. 5). Based on the prediction of arsenic pathways through the 16S rRNA gene sequences, arsenic detoxification mechanisms predominated. Since the arsenic cycle in soils is associated with the nitrogen cycle (Feng et al., 2023), genes involved in nitrogen metabolism enzymes were also retrieved (Table S5). Nitrate reduction was highly detected across all soil samples.

Metagenomics sequencing generated a high volume of sequences which were subsequently quality-verified and trimmed of adapter sequences (Table S6). Thirteen arsenic-functional genes in the soil microbiomes were detected and their abundance were normalized by the FPKM method. The results indicated that the overall proportion of arsenic-functional genes in soils of the dry and wet seasons were comparable (Fig. 6). Both arsenic detoxification (e.g., arsC, arsA, and acr3) and metabolism (e.g., aioA, aioB, and arrA) pathways were identified: however, the relative abundance of arsenic detoxification genes was higher than arsenic metabolism genes. Among the four key functional marker genes, the relative abundance of the arsC gene was greatest, followed sequentially by the arsM, aioA, and arrA genes. While the metagenomics analysis of agricultural soil microbiomes identified all four key arsenic-functional genes, the PICRUSt2 prediction tool discovered only the arsC and the aioA genes (Figs. 5 and 6). Based on shotgun metagenomics analysis, detoxification arsenate reduction and arsenic methylation pathways may cooperatively mitigate arsenic levels in the soil.

Agricultural soil microbiomes and arsenic-transforming bacteria

The levels of arsenic, cadmium, and copper detected in the soils analyzed in this study were within the range of background levels found in soils elsewhere (Sabet Aghlidi et al., 2020; Liu et al., 2022). Previous studies indicated that zinc concentrations in agricultural soils ranged from 10 to 300 mg/kg (Barber, 1995; Noulas, Tziouvalekas & Karyotis, 2018). However, the extremely high concentration of zinc in the T1 sample cannot be clearly explained due to the lack of sufficient information. The co-occurrence of cadmium, copper, zinc, and arsenic could impact their availability in the soils due to competitive adsorption on soil particles (Lu & Xu, 2009; Gong et al., 2020). However, the specific interactions among these elements are complex, and their availability can vary depending on specific environmental factors and soil conditions. While iron is another element associated with arsenic mobility in soil, a prior study indicated that, in comparison to the clay fraction, both iron and organic matter have a marginal influence on the overall stability and accumulation of arsenic in agricultural soils (Dousova et al., 2016). To assess the impact of these metal elements, including arsenic, on human health, further analysis of their concentrations in banana products is warranted. Low levels of soil nutrients (i.e., TOC, TP, and TN) suggest that they may have been absorbed by crops. Additionally, the absence of nitrogen in both soil samples could be due to microbial nitrification and denitrification (Del Prado et al., 2006), as supported by the high abundance of nitrogen-related functional enzymes, including reductase and ammonia monooxygenase (Table S5). Nitrospira, nitrite-oxidizing bacteria, play key roles in nitrification, were also detected in the analyzed soils (Fig. 2A) (Daims & Wagner, 2018).

The concentrations of arsenic in both soil samples analyzed in this study were in the range of those found in background soils (∼2 to 20 mg/kg) (Xiao et al., 2016; Dunivin, Miller & Shade, 2018; Zhang et al., 2021). A previous study showed that coarse sand particles enriched arsenic, while clay particles were the primary source of increased arsenic in soils (Zou et al., 2023). Among different soil particle size fractions, clay harbored the highest abundance of arsenic functional genes, such as the arrA, arsC, and arsM genes (Zou et al., 2023). The high proportion of clay in our analyzed soils may have contributed to the release of arsenic from the soil, favoring the occurrence arsenic-transforming bacteria. Although the arsenic concentrations of the soils in both dry and wet seasons were relatively low, arsenic-transforming bacteria were commonly detected. Although it is complicated to compare taxonomic profiles resulting from the analysis of the 16S rRNA gene and shotgun metagenomics due to differences in technical issues and databases, the top three phyla highly detected by both analyses were Actinobacteria, Proteobacteria, and Firmicutes, respectively (Fig. 1A and Fig. S3). Arsenic-transforming bacteria are ubiquitous in terrestrial environments, especially paddy soils, and arsenate-resistant bacteria affiliated with Proteobacteria, Bacteroidetes, and Firmicutes have been widely isolated from soils without the detection of arsenic (Jackson, Dugas & Harrison, 2005; Xiao et al., 2016; Zhang et al., 2021). Moreover, the major bacterial taxa (i.e., Actinobacteriota, Proteobacteria, Firmicutes) found in this study were highly dominant in banana plantation soils, possibly contributing to the degradation of various organic compounds (Chou et al., 2017; Kaushal et al., 2022).

Agricultural soil microbiomes analyzed in this study were comparable to those previously reported. However, the proportion of dominant phyla varied across soils due to differences in soil characteristics, nutrient loads, and vegetation (Zhao et al., 2019; Liu et al., 2021). The abundance of Chloroflexi, Gemmatimonadota, and Bacteroidota showed no statistical difference between the two seasons, indicating their succession throughout the year. A previous study demonstrated that while the relative abundance of Gemmatimonadetes significantly increased, Chloroflexi and Bacteroidetes in soils were not significantly different during restoration stages (Liu et al., 2021). The phyla Chloroflexi and Bacteroidetes may constitute a core soil microbiome sustaining soil health.

Common bacterial taxa highly detected in both T1 and T2 soils were Bacillus, uncultured Roseiflexaceae, and Streptomyces (Fig. 1B and Fig. S4). Based on the presence of the arsC, arrA, arsM, and aioA genes, members of the genus Bacillus have versatile metabolisms in arsenic transformations (Afkar et al., 2003; Huang et al., 2018; Marwa et al., 2019; Mujawar, Vaigankar & Dubey, 2021). Bacillus sp. has been suggested for the purpose of arsenic bioremediation (Kabiraj et al., 2022). Although the association between Roseiflexaceae and arsenic transformation is not known, this genus is involved in carbon dioxide assimilation and has a positive impact on plant growth (Shi et al., 2020). Roseiflexaceae are one of the five keystone taxa found in a long-term fertilized soils used for maize cultivation, and are responsible for enhancing phosphorus flow and preventing the movement of toxic aluminum and manganese from the soil to the crops (Wang et al., 2022). Streptomyces sp., previously isolated from the rice rhizosphere, harbors the arsM gene, indicating its role in arsenic methylation (Kuramata et al., 2015). In addition, both Bacillus sp. and Streptomyces sp. were detected at high abundance in healthy banana plantation soils and they are associated with the suppression of banana pathogens (Zhou et al., 2019; Jamil et al., 2022; Kaushal et al., 2022), possibly facilitating banana growth. The genus Microvirga was also consistently detected across the soil samples (Fig. 1B and Fig. S4). Members of this genus are distributed across a broad range of environments, including those as diverse as polar soils and hot spring sediments, and they are reported to exhibit both arsenate reduction and arsenite oxidation ability (Tapase et al., 2017; Liu et al., 2020; Zhu et al., 2021). The genus Nocardioides was highly abundant, as detected by shotgun metagenomics (Fig. S4). Verification of its arsC gene indicated its ability in arsenate reduction (Bagade et al., 2016).

The abundance of Crenarchaeota, especially Candidatus Nitrocosmicus, Nitrososphaeraceae, and Candidatus Nitrososphaera, was significantly higher in T2 than T1 (Fig. 1). Members of Crenarchaeota are commonly found in soils, playing an important role in ammonia oxidation (Gornish et al., 2020; Behnke et al., 2021). These archaeal taxa potentially drive the nitrogen cycle in the agricultural soils subjected to routine fertilizer application. A previous study showed that among the detected archaeal taxa, Candidatus Nitrocosmicus and Nitrososphaeraceae were dominant across topsoil (0–10 cm depth) of ten years of restoring eucalypt woodland (Yan et al., 2020). Although Candidatus Nitrocosmicus, Nitrososphaeraceae, and Candidatus Nitrososphaera were also present in the top soils of our dry and wet season samples, they were more abundant in T2 than T1 (Fig. 1B). The increased relative abundance of these archaeal taxa was likely associated with the substantial reduction of zinc concentration in the T2 soil (Table 1). Increasing zinc application to 500 mg/kg in topsoil (0–15 cm depth) inhibited the transcriptional response of archaeal ammonia oxidation (Vasileiadis et al., 2012).

Among the four key arsenic-functional genes representing respective arsenic-transforming bacterial communities, the arsC gene was undetectable by PCR amplification possibly due to the limitation of specific primers used in this study. However, the presence of the arsC was retrieved from the analysis of 16S rRNA gene sequences through the PICRUSt2 software and shotgun metagenomics (Figs. 5 and 6). It is also possible that the arsC gene present in the analyzed soil is likely related to those of uncultured bacteria, making it undetectable using the available specific arsC primers.

The phylogenetic trees of aioA, arrA, and arsM revealed that the majority of aioA, arrA, and arsM sequences recovered from the analyzed soils were closely related to uncultured sequences previously found in diverse environments (Figs. 2, 3 and 4). This indicates the ubiquity of arsenic-transforming bacteria in the environment. However, due to the limitation of culture techniques, further exploration is still needed to gain in-depth knowledge of their association with arsenic metabolisms. Our results also showed that a few minor retrieved aioA sequences from the T2 sample were closely related to the known aioA sequence of Microvirga sp. The genus Microvirga, previously isolated from polluted soils, is capable of performing arsenite oxidation through the aioA gene (Tapase et al., 2017; Tapase & Kodam, 2018). Some minor detected arrA sequences from the T1 sample were associated with the known dissimilatory arsenate-reducing bacteria, Geotalea uraniireducens and Desulfuromonas sp. WB3, which were respectively isolated from paddy field soil and anoxic sediment (Osborne et al., 2015; Qiao et al., 2018). Rhodopseudomonas palustris was identified as the sole known arsenic methylation bacterium, exhibiting a close relationship with the arsM sequences from both T1 and T2 samples. Rhodopseudomonas palustris, well-known arsenic methylation bacterium harboring the arsM gene, was commonly found in paddy soils (Xiao et al., 2016; Afroz et al., 2019).

Potential arsenic-transforming pathways in agricultural soils

The presence of potential arsenic-transforming pathways in agricultural soils were predicted by both the analysis of 16S rRNA gene sequences and by shotgun metagenomics. Both analyses indicated that arsenic detoxification, as evidenced by the significant detection of a large proportion of the ars operon, played a pivotal role in driving the arsenic cycle in the agricultural soil (Figs. 5 and 6). In other arsenic-contaminated and uncontaminated soils (arsenic levels < 20 mg/kg) analyzed by metagenomics, arsenic detoxification overcame arsenic energy metabolisms (Xiao et al., 2016; Dunivin, Yeh & Shade, 2019). Based on the identified detoxification arsenic-functional genes, arsenate reduction conferred by the arsC gene was a dominant pathway. The arsC mechanism is the most extensive-studied arsenic detoxification and resistance. Arsenate exhibited maximal adsorption on clays under low pH conditions, but its adsorption capacity decreased as the pH increased (Goldberg, 2002). Given that arsenate serves as a substrate for arsenate reduction, it could be released from the soils analyzed in this study, enhancing the abundance and activity of the arsC gene. The properties of our analyzed soils, including pH and texture, could enhance the availability of arsenate for arsenate-reducing bacteria. The abundance and expression of the arsC gene in surface soils were generally high (Luo et al., 2014; Xiao et al., 2016). By analyzing the genome of aerobic arsenate-reducing bacteria isolated from arsenic-contaminated soil, together with the arsC gene, the arsH, arsB, and arsR genes controlling the regulation of arsenate reduction were identified (Tian & Jing, 2014). The detoxification arsenate reduction involves the cytoplasmic reduction of arsenate to arsenite by the arsC gene. Then, the produced arsenite is transported out of the bacterial cell through a specialized arsenite-pump which were mediated by the arsA/arsB or acr3 genes (Singh et al., 2021). The large proportion of arsC, arsA, arsB, and acr3 genes identified by both the PICRUSt2 program and the FPKM method supported that detoxification arsenate reduction within the soil was an important arsenic-transforming pathway. Cultivation-independent studies proposed that the arsC (grx) gene was less tolerant to high arsenic concentrations compared to the arsC (trx) gene, and increased levels of arsenic in a specific area led to a higher prevalence of the arsC (trx) gene (Escudero et al., 2013; Kurth et al., 2017). However, in this current study, the relative abundance of the arsC (trx) gene was higher in soils with low arsenic concentrations than the arsC (grx) gene (Fig. 6). The arsC (trx) gene is commonly found in Gram-positive bacteria (Mukhopadhyay et al., 2002). The high abundance of the genera Bacillus and Streptomyces, Gram-positive bacteria, possibly contributed to the higher relative abundance of the arsC (trx) than the arsC (grx) genes. The high abundance of arsR gene detected by the PICRUSt2 program possibly because the arsR genes across various bacterial strains are highly conserved (Xu et al., 1998). Overall results suggest that the detoxification arsenate reduction pathway is widespread in soils, both those with arsenic contamination (88.93–646.61 mg/kg) and those without (less than 15 mg/kg) (Xiao et al., 2016; Zhou et al., 2022).

Metagenomic analysis demonstrated that the second most prevalent arsenic-transforming pathway was arsenic methylation mediated by the arsM gene (Fig. 6). Biogeography of soil microbiomes revealed that the arsM gene was unexpectedly abundant, especially in cultivation-independent samples (Dunivin, Yeh & Shade, 2019). In paddy soils with low arsenic concentrations (less than 15 mg/kg) and polluted soils with extremely high arsenic concentration (88.93–646.61 mg/kg), arsenic methylation was also identified as the second most prominent process, following detoxification through arsenate reduction (Xiao et al., 2016; Zhou et al., 2022). Previous studies demonstrated that arsenic methylation was active under aerobic, anaerobic, and facultative conditions (Wang et al., 2014). In silico analysis also suggested that the arsM gene sequences are widespread in the environment and are conserved across various living organisms (Kabiraj et al., 2023). The arsM gene catalyzed the conversion of arsenite to volatile methyl arsenic to complete arsenic methylation process (Qin et al., 2006). Consequently, in surface soil, detoxification through arsenate reduction leads to the production of arsenite, which is subsequently methylated by arsenic methylation process. This suggests a cooperative mechanism between detoxification arsenate reduction and arsenic methylation in the surface soil from the agricultural area with low arsenic level. Significant coexistence of the arsC and arsM genes was previously found in paddy soils (Zhang et al., 2015; Yang et al., 2020). The mechanism of detoxification arsenate reduction and arsenic methylation was also identified in a plant body, primarily in its root (Ma et al., 2016).

The genes involved in methylation and reduction processes exhibited higher abundance compared to those associated with the oxidation process (Figs. 5 and 6). Arsenate is typically more abundant in soils compared to arsenite. Since arsenic methylation genes are commonly found in soils (Dunivin, Yeh & Shade, 2019), they may compete with arsenite oxidation genes for the arsenite produced during arsenate reduction. Additionally, arsenite is more mobile compared to arsenate. The arsenite produced through arsenate reduction could potentially leach into the aquatic environment. Consequently, due to the available arsenic substrate and the prevalence of arsenic methylation genes in soils, both arsenate reduction and arsenic methylation are likely to play roles in arsenic transformations in soils.

Since the abundance of genes mediating nitrate reduction mechanisms was highly detected (Table S5), it suggests the potential for arsenite oxidation coupled with nitrate reduction process. Due to the soil texture, the presence of clay could create anoxic conditions, facilitating the conditions for arsenite oxidation coupled with nitrate reduction. Previous studies suggested that co-occurrence of arsenite oxidation and nitrate reduction was an important process controlling arsenic availability in paddy soils (Li et al., 2019; Feng et al., 2023).

While shotgun metagenomics identified all four arsenic-transforming mechanisms: detoxification arsenate reduction (arsC), arsenite oxidation (aioA), arsenic methylation (arsM), and dissimilatory arsenate reduction (arrA), the analysis of 16S rRNA gene sequences found only the first two pathways. The use of shotgun metagenomics allows for broader detection of environmental patterns and variations in the abundance of arsenic-functional genes, as it targets all genes irrespective of primer coverage and specificity. Although the accuracy of PICRUSt2 in predicting microbial functions has improved, it is more suitable for human microbiome than soil microbiome samples (Langille et al., 2013; Sun, Jones & Fodor, 2020). Consequently, some functional genes, including those associated with arsenic, can be overlooked by using the PICRUSt2 prediction tool.