High Arsenic Levels Increase Activity Rather Than Diversity or Abundance of Arsenic Metabolism Genes in Paddy Soils

doi:10.21203/rs.3.rs-142737/v1

Download PDF

Research

High Arsenic Levels Increase Activity Rather Than Diversity or Abundance of Arsenic Metabolism Genes in Paddy Soils

https://doi.org/10.21203/rs.3.rs-142737/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 28 Sep, 2021

Read the published version in Applied and Environmental Microbiology →

Version 1

posted 13 Jan, 2021

You are reading this older preprint version

Read the latest preprint version →

Background

Arsenic (As) is a toxic metalloid pervasively present in the environment. Microorganisms have evolved the capacity to metabolize As, and As metabolism genes are ubiquitously present in the environment even in the absence of high concentration of As. However, the As metabolism genes diversity, relative abundance, and transcriptional activity in response to different As levels remain unclear, limiting our understanding of the microbial activities that control the fate of an important environmental pollutant. To address this issue, we applied metagenomics and metatranscriptomics to paddy soils showing a gradient of As concentrations to investigate As resistance genes (ars), including arsR, acr3, arsB, arsC, arsM, arsI, arsP and arsH as well as energy-generating As respiratory oxidation (aioA) and reduction (arrA) genes.

Results

Somewhat unexpectedly, the relative DNA abundances and diversity of ars, aioA and arrA genes were not significantly different between low and high (~10 vs ~100 mg kg^-1) As soils. By comparison to available metagenomes from other soils, geographic distance rather than As levels drove the different composition of microbial communities. Arsenic significantly increased ars genes abundance only when its concentration was higher than 410 mg kg^-1. In contrast, between low and high (~10 vs ~100 mg kg^-1) As soils, metatranscriptomics revealed a significant increase in transcription of ars and aioA genes, which are induced by arsenite, the dominant As species in paddy soils, but not arrA genes, which are induced by arsenate. Co-occurrence patterns of arsR, acr3, and arsM genes were revealed by network analysis corroborating that the arsR, acr3 and arsM genes are usually organized in a single ars operon. The transcriptome level response appeared to be community-wide as opposed to taxon-specific.

Conclusions

High As levels increased the activity of As metabolism genes rather than their abundance or diversity in paddy soils. These findings advance understanding of how microbes respond to high As levels and the diversity of As metabolism genes in paddy soils, and indicated that future studies of As metabolism in soil, and likely other environments, should include the function (transcriptome) level.

General Microbiology

Arsenic resistance genes

Arsenic respiratory genes

Abundance and diversity

Transcriptional activity

Network

Paddy soil

Arsenic (As) is a toxic metalloid, widely distributed in Earth’s crust. Relatively high amounts of As are introduced into the environment from nature sources such as volcanic eruptions, weathering of rocks and geothermal activities as well as anthropogenic activities such as mining, pigment production and application of As-based pesticides in agriculture [1, 2]. Because of its pervasive presence in the environment, microorganisms have evolved the capacity to metabolize As [3], and this capacity is hypothesized to be an ancient mechanism that emerged at least 2.72 billion years ago [4, 5].

Microbial mediated As biotransformation include As detoxification to mitigate toxicity and As respiration to generate energy. Several mechanisms have evolved to detoxify both inorganic and organic As. The inorganic As compounds can be detoxified by cytoplasmic arsenate [As(V)] reduction, which is catalyzed by a reductase (ArsC) [6], followed by arsenite [As(III)] efflux outside the cell, which is catalyzed by two evolutionarily unrelated As(III) efflux permeases, i.e. ArsB and Acr3 [7]. Additionally, As(III) can be methylated by an As(III) S-adenosylmethionine (SAM) methyltransferase (ArsM) to the more toxic organic methylarsenite [MMAs(III)] [8], followed by the organoarsenical detoxification pathways, which have been recently shown to be present in various microbial genomes. The latter detoxification pathways include MMAs(III)’s oxidization to methylarsenate [MMAs(V)], which is catalyzed by the MMAs(III)-specific oxidase ArsH [9], and actively transported outside the cell by the MMAs(III) efflux permease ArsP [10], or MMAs(III)’s demethylation to the less toxic As(III) by the ArsI lyase, which cleaves the carbon-As bond in MMAs(III) [11]. These As resistance genes are usually organized in an ars operon, which is controlled by one of three known ArsR transcriptional repressors (ArsR1, ArsR2 and ArsR3) regulated selectively by As(III) and one (ArsR4) by MMAs(III), depending on the exact version of the operon [12, 13].

The As respiratory pathways include oxidation of As(III) either to detoxify it to less toxic As(V), or couple it to ATP production [14], and dissimilatory reduction that couples As(V) reduction to anaerobic heterotrophic growth [15]. As(III) oxidation is catalyzed by As(III) oxidase (Aio) whose large subunit (AioA) has been well characterized in both bacteria and archaea [16]. On the other hand, As(V) reduction is catalyzed by respiratory As(V) reductase (Arr), whose large catalytic subunit (ArrA) serves as a reliable marker for the process [17].

Due to the anaerobic conditions typically prevailing in paddy soils, As(III) is commonly the predominant form of As, which is more mobile and bioavailable to microbes for biotransformation [18–20]. Indeed, microbial genes involved in As resistance such as arsC and arsM, and genes involved in As respiratory such as aioA and arrA have been detected to be widely presented in the paddy soils with As contamination [21], and their activity to biotransform As has also been experimentally demonstrated [22–24]. Recent studies of soils with low As levels, i.e. typically less than 15 mg kg^− 1 (or ppm), also showed relative high abundance and diversity of microbial As metabolism genes; most notably, arsB, acr3, arsH, arsR, arsC, arsM, aioA and arrA genes in paddy soils [25] and arsC, arsM, aioA and arrA genes in estuarine sediments across southeastern China [26]. These results reveal high prevalence of As metabolism microbes in the environment even in absence of high concentration of As. However, these previous studies were carried out at the DNA level and thus, the activity of the As metabolism genes detected remains speculative. Further, these DNA-based results showing prevalent As metabolism microbes in both high and low As level environments are somewhat surprising because gene functions and species selected by a specific factor tend to be more abundant in ecosystems where the factor (selective pressure) is stronger as this has been shown for various soil factors such as pH, total organic carbon (TOC), total nitrogen (TN), phosphorus (TP) [27] and organic pollutants [28].

Herein, we applied both metagenomics and metatranscriptomics to investigate the resident (at the DNA level) and active (at the RNA level) As-associated microbial genes in paddy soils of different As concentrations, aiming to address the following questions: (1) Do As metabolism genes in high As levels paddy soils show higher relative abundances and/or transcriptional activities than their counterparts in low As levels paddy soils to cope with the higher As level present? (2) Which are the key As metabolism genes and pathways coping with high As levels? (3) What are other factors than As levels driving microbial communities compositional and gene transcriptional differences in soils with As contamination?

Sites description and sample collection

Soil samples were collected from six distinct paddy fields in Hunan province, China, including three sites with high As levels (67.3–104.0 mg kg^− 1) in Cili (CL_AsHig), Chenzhou (CZ_AsHig) and Shuikoushan (SKS_AsHig), and three sites with low As levels (2.5–10.8 mg kg^− 1) in Hanpu (HP_AsLow), Lianhua (LH_AsLow) and Yuchangping (YCP_AsLow). Eighteen soil samples, three replicates from the same field for each site, were collected from the soil surface (0–20 cm) in July 2016 during the rice-growing period (flooded). Paddy soil samples were placed in sterile plastic bags and transported to the laboratory on ice for DNA and RNA extraction and soil physiochemistry analysis. Soil properties, including pH, Eh, total concentration of carbon (TC), nitrogen (TN), phosphorus (TP), sulfur (TS), arsenic (As), lead (Pb), cadmium (Cd), chromium (Cr), antimony (Sb) and iron (Fe) were determined on a composite sample of the three replicates following standard methods [21]. Details of the underlying methods are provided in the Supplementary methods.

Soil DNA and RNA isolation, library preparation and sequencing

Soil DNA and total RNA was extracted from 0.5 g and 2 g of soil using the FastDNA SPIN kit (MP Biomedicals) and E.Z.N.A Soil RNA Midi kit (Omega, Bio-Tek Inc., USA), respectively, according to the manufacturer’s protocol. For each site, DNA or total RNA extracted from the triplicates were pooled together in order to obtain enough high quality material for each site and to overcome high sample-to-sample heterogeneity, which is characteristic of soil microbial communities. The quality of DNA and total RNA was assessed by NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Inc., Wilmington, DE, USA) and the ratio of absorbance at 260 nm and 280 nm was around 1.8–2.0, indicating good quality for downstream analysis. Total RNA of the paddy soils was subjected to an rRNA removal procedure using the Ribo-Zero rRNA removal bacteria kit (Illumina Inc., San Diego, CA, USA) following the manufacturer’s protocol.

Paired-end DNA and cDNA libraries (2 × 150 bp) were constructed using TruSeq™ DNA Sample prep kit (Illumina Inc., San Diego, CA, USA) and TruSeq™ RNA Sample prep kit (Illumina Inc., San Diego, CA, USA), respectively, following the manufacturer’s protocol. The six DNA and cDNA libraries were sequenced respectively on an Illumina HiSeq 3000 platform at Majorbio Bio-pharm Technology Company in Shanghai, China with HiSeq 3000/4000 PE Cluster kit and HiSeq 3000/4000 SBS kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s protocol. Metagenomics and metatranscriptomics datasets have been deposited in National Center for Biotechnology Information (NCBI)’s Short Read Archive (SRA) database, under the bioproject PRJNA616041. NCBI SRA numbers for all datasets were provided in Supplementary Table S1.

Taxonomic classification, functional annotation and population genome binning of soil metagenomes and metatranscriptomes

The raw paired-end reads of metagenomic and metatranscriptomics datasets were trimmed using SolexaQA [29], removing reads with Phred quality score (Q < 20) and minimum fragment read length of 50 bp after trimming for downstream analyses. Parallel-META version 2.4 [30] was used with default settings to recover 16S rRNA gene fragments from metagenomes and to identify and remove residual rRNA sequences after rRNA subtraction from metatranscriptomes. 16S rRNA gene fragments from metagenomes were then processed for Operational Taxonomic Unit (OTU) picking using close OUT picking, defined at the 97% nucleotide sequence identity threshold, and taxonomic identification with the SILVA database [31] as implemented in QIIME 1.9.1 [32] and classified at the phylum and class levels. The relative abundance of each taxon at each level was normalized by the total number of 16S rRNA gene sequence fragments recovered from the corresponding metagenome.

Protein-coding sequences from short-read metagenomes and metatranscriptomes were predicted using FragGeneScan [33] with default settings and searched against Swiss-Prot database (UniProt, downloaded in April 2019) using BLASTp (BLAST + version 2.2.28) [34]. Only the best match with amino acid identity ≥ 50% and query sequence length coverage ≥ 70% by the alignment was retained. Protein-coding sequences were also functionally annotated, using the same threshold, against the SEED database using the subsystem categories [35]. Differentially abundant taxa based on 16S rRNA gene sequence fragments and SEED subsystems of metagenomics and metatranscriptomics between high and low As levels were determined with the DESeq2 package using the negative binomial model, adjusted for false discovery rate [36]. Heatmap in pheatmap R package (R version 3.5.1) [37] was used to visualize the significantly differentially abundant taxonomic levels and SEED subsystems (adjusted p value < 0.05). Metagenome-assembled genomes (MAGs) were recovered using both MaxBin Version 2.1.1 [38] and MetaBAT Version 2.12.1 [39] as described in more detail in Supplementary methods.

Assessment of the abundance and transcriptional activity of arsenic metabolism genes

Building an As metabolism protein reference database

To identify and quantify As metabolism proteins containing reads, in-house databases of eight As resistance proteins [Ars i.e. ArsR (ArsR1, ArsR2, ArsR3 and ArsR4), ACR3, ArsB, ArsC (ArsC1, ArsC2 and ACR2), ArsM, ArsI, ArsP, ArsH] and As respiratory oxidation (AioA) and reduction (ArrA) proteins were constructed and manually curated as described before [5]. Briefly, verified arsenic metabolism protein sequences were downloaded from UniProt or NCBI’s NR database to build HMM profiles, and searched against 786 representative species which are sampled from an original dataset used for reconstruction of tree of life (species tree) [40], to identify homologs of As metabolism proteins. The resulting matching sequences were manually inspected for the presence of conserved domains based on the Pfam models (release 28.0) [41] and forming deep-branching clades in the phylogenetic tree of the (manually) verified reference sequences to further determine if the sequence in question was truly As metabolism protein. The tree was built using maximum likelihood as implemented in RAxML v8.4.1 with the PROTGAMMAAUTO model [42]. Sequences forming long-branches and/or lacking one or more of the conserved domains were removed from further analysis.

To build a more comprehensive As metabolism proteins database that better captures the diversity of As metabolism proteins, proteins on the contigs assembled using IDBA with default settings [43] from metagenomes and metatranscriptomes (for more details, see Supplementary methods) were also searched against the curated As metabolism proteins database using BLASTp (BLAST + version 2.2.28) with a cut-off for a match of amino acid identity ≥ 40% and reference length coverage ≥ 70% by the alignment. A total of 1040 Ars proteins, including 87 ACR3, 281 ArsC, 1 ArsH, 69 ArsI, 173 ArsM, 21 ArsP and 408 ArsR were identified using this approach. However, for the AioA and ArrA proteins that are longer than the Ars proteins (800–900 vs 200–400 amino acids), no match was identified, and it was likely attributable to the limiting assembly, consistent with the relatively low sequencing effort. Therefore, in order to obtain a comprehensive database for AioA and ArrA proteins, the AioA and ArrA references in curated As metabolism database were searched against UniRef90 (TrEMBL, downloaded in October 2018) database for additional homologs. A more stringent cut-off compared to the Ars proteins mentioned above, i.e. identity ≥ 40% and reference length coverage ≥ 90% by the alignment, was used to filter the AioA and ArrA proteins in UniRef90 database. Identified Ars protein sequences from metagenomes and metatranscriptomes, and AioA and ArrA proteins sequences from UniRef90 database were searched against Pfam database (release 32.0) [41] in Jun 2019 and inspected for the presence of known conserved domains. The resulting matching sequences were aligned using MAFFT version 7.407 [44] and visually inspected for forming deep-branching clades in a maximum likelihood phylogenetic tree as described above to determine if the sequences were truly As metabolism proteins. These phylogenetic trees were also used for placement of As metabolism proteins encoding reads as described below.

Estimation of As metabolism genes abundance and transcriptional activity

Short sequences from metagenomes and metatranscriptomes were searched against the comprehensive Ars, AioA and ArrA references databases using BLASTx search (BLAST + version 2.2.28) and filtered for best matches with amino acid identity ≥ 90%, alignment length ≥ 25 amino acid and e-value ≤ 1E-5. MicrobeCensus was used to assess the average genome size of each sample community [45] and normalize the total base-pairs assign to each As gene by gene length, i.e. RPKM = (the number of reads mapping on reference sequences/gene length in kb per million reads), and genome equivalent, i.e., RPKM/average genome size, providing an estimate RPKG [(reads mapped to gene)/(gene length in kb)/(genome equivalents)] of how many cells from those sampled encode the gene of interest. Relative abundance of As metabolism genes in metatranscriptomics was calculated by the statistic RPKM (reads per kb per millions of reads), i.e. RPKM = (the number of reads mapping on reference sequences)/(gene length in kb)/( the total number of reads after removal of rRNA reads). The As metabolism genes transcriptional activity was evaluated by the relative abundance of As metabolism genes in metatranscriptomes (RPKM) normalized by the relative abundance of the same gene variant or OTU (at > 90% amino acid level) identified in metagenomes (RPKM; we used RPKM to normalized the relative abundance of the same gene identified in metagenomes here instead of RPKG in order to keep the normalization criteria consistent between metagenomes and metatranscriptomes in calculating the gene’s transcriptional activity), and visualized in heatmap using the pheatmap package in R 3.5.1 [37]. One-way ANOVA analysis, performed with the vegan package [46] in R 3.5.1, was used to compare the difference of As metabolism genes abundances and transcriptional activities in high vs low As levels paddy soils.

Phylogenetic placement of As metabolism gene carrying reads

Arsenic metabolism protein-containing reads identified by BLASTx search were firstly translated to amino-acid sequences by FragGeneScan [33] with default settings. The amino acid fragments were then added to the As metabolism protein references alignments using MAFFT version 7.407 [44] with ‘addfragments’ option and were placed in the corresponding phylogenetic tree of different As metabolism proteins (described above) respectively, using RAxML v8.4.1 [42] with -f v option. The placement of the target As metabolism protein-containing reads was visualized in iTOL [47] after processing the resulting visualization jplace file generated by RAxML with JPlace.to_iTOL.rb from the Enveomics Collection [48].

Network analysis

To examine potential co-occurrence patterns and associations between As metabolism genes, the inferred correlations between the abundances of the As genes based on metatranscriptomic data were calculated using Sparse Correlations for Compositional data (Sparcc), which predicts interpretable microbial interaction networks through graphical model inference [49]. Only As genes detected in all triplicate samples of the same As level (i.e. low vs. high) were used in the correlation analysis. This preliminary filtering step removed poorly represented As metabolism genes that occurred in a limited number of samples and thus, reduced spurious associations and artifacts [50]. The nodes in the resulting Sparcc network represent the individual As metabolism genes, and the edges correspond to significant associations between the nodes. The network was explored and visualized with the interactive platform Gephi [51]. Associations examined and reported below were only those based on correlation coefficients ≥ 0.6 and statistical significance at p < 0.01.

Alpha and beta community diversity estimation

The abundance-weighted average coverage of the sampled microbial communities achieved by our sequencing effort was estimated by Nonpareil [52, 53] based on the level of redundancy of the sequence reads of metagenomes. Mash distances [54] were used to assess the sample-to-sample sequence composition similarity at the whole metagenome level, based on shared kmers. Because the coverage of these metagenomics varied among different studies, the datasets were subsample to 2Gbp to make them comparable. Principle coordinates analysis (PCoA) and/or non-metrical multidimensional scaling (NMDS) plot were used to visualize the site-to-site similarity matrices of microbial community compositions based on Mash distances, OTU tables, and annotation counts from SEED subsystems based on Bray-Curtis distance. NMDS plot of annotation counts of ars [i.e. arsR (arsR1, arsR2, arsR3, arsR4), acr3, arsB, arsC (arsC1, arsC2 and acr2), arsM, arsI, arsP, arsH], aioA, arrA genes from metagenomes and metatranscriptomes based on Bray-Curtis distances were used to visualize the site-to-site similarity of As-associated microbial community compositions.

The Adonis test was used to determine correlations between physicochemical variables measured at each site and similarity in microbial community composition (Mash distance) or ars genes abundances/transcriptional activities (Bray-Curtis distances) and tested for significance based on a permutation test with 999 interactions. The interactions between different physicochemical variables were assessed by PERMANOVA analysis while the correlation of physicochemical variables with aioA and arrA gene abundances and transcriptional activities in different sites were tested by Pearson correlation. These analyses were performed in R 3.5.1 with the vegan [46] and ggplot2 [55] packages.

Soil physicochemical characteristics

A total of 18 samples, three biological replicates from the same field for each site, were collected from six paddy soils in Hunan, China with different As concentrations (Fig. 1). Arsenic (As) concentrations in these paddy soils varied significantly (ANOVA, p < 0.01) between high (67.3–104.0 mg kg^− 1) and low (2.5–10.8 mg kg^− 1) As levels. The high As concentrations in the paddy soils likely resulted from the historic mine activities occurring in the same area. The total concentrations of nitrogen (TN) in high As paddy soils were also significant higher (ANOVA, p < 0.01) than those with low As (1.64–1.77 g kg^− 1 vs 1.17–1.24 g kg^− 1). No significant difference was found in the other soil physiochemical characteristics, including pH, Eh, total concentration of carbon (TC), phosphorus (TP), sulfur (TS), and concentrations of heavy metals including Pb, Cd, Cr, Sb and Fe. Details of the soil properties were summarized in supplementary Table S2.

Broad microbial community diversity indices

Six shotgun metagenomes and six shotgun metatranscriptomes were obtained after the triplicate DNA and RNA samples from each site were pooled together, respectively. For metagenomes (Supplementary Table S3), 12–15 Gbp (billion base pairs, 13 Gbp, on average) of shotgun metagenomic data and an average read length of ~ 115 bp were acquired after trimming for each pooled DNA sample. On average, 6.6 Gbp (ranged from 3.9 to 8.3 Gbp) of shotgun metatranscriptomic data and average read length of ~ 129 bp were obtained after trimming and removal of rRNA reads for each pooled RNA sample (Supplementary Table S4). The coverage achieved by sequencing based on the read redundancy value estimated by Nonpareil was around 22–27%, except for sample SKS_AsHig, which had a lower alpha-diversity (measured by Nonpareil) than the others (Nd 22.6 vs 23.6–23.8; note that Nd is in log scale) and thus, higher coverage (42%; Supplementary Table S3). No significant difference (p > 0.05) was found in the alpha-diversity between high and low As paddy soils (average of 23.3 vs 23.7). MicrobeCensus was used to evaluate the average genome size of microbial community in paddy soils and revealed no significant shifts due to the different As levels (5.7 Mbp vs 5.6 Mbp, on average, p > 0.05). Considering the difference of average genome size of microbial community between these paddy soils was no more than 2-fold, estimation of As metabolism genes abundance was essentially the same between the RPKG and RPKM metrics.

Paddy soil community taxonomic and functional composition

16S rRNA gene (16S) fragments extracted from the metagenomes were about 0.03% – 0.04% of the total reads (Supplementary Table S3) and around 65% – 70% of them were classifiable at the phylum level. Bacteria were predominant in these paddy soils, accounting for 61% – 69%, while Archaea accounted for only 0.5% – 2.4%; the remaining reads were unclassified against the SILVA database [31]. The dominant bacteria comprised of phyla Proteobacteria, Acidobacteria, Chloroflexi, Nitrospirae, Verrucomicrobia and Bacteroidetes, representing 50% – 59% of the total 16S-carrying reads (Fig. 2A). Archaea community was mostly represented by the Crenarchaeota phylum at 0.05% – 0.78% of the total, followed by Euryarchaeota at 0.36% – 0.77% (Fig. 2A). There were overall no significant relationships between the level of As and community composition, assessed by Bray-Curtis similarity in 16S rRNA gene-based OTUs (NMDS plot; Supplementary Figure S1). DESeq2 analysis revealed statistically significant lower abundance of members of Euryarchaeota, i.e. Methanomicrobia and Thermoplasmata (0.15% vs 0.85% and 0.02% vs 0.07%, respectively; p < 0.05) in high vs low As paddy soils (Fig. 2B). In contrast, the abundance of Chloroflexi, mostly comprised of Ellin6529 (1.44% vs 0.52), S085 (0.28% vs 0.09%), Thermomicrobia (0.06% vs 0.01%), and Chlorobi-OPB56 (0.16% vs 0.05%), GAL15 (0.21% vs 0.01%) and OP11-WCHB1-64 (0.07% vs 0.01%) were significantly (p < 0.05) higher in high vs low As paddy soils (Fig. 2B).

Protein-coding fragment sequences were predicted for about 90% – 92% and 44% – 59% of the trimmed reads in metagenomics and metatranscriptomics, respectively (Supplementary Table S3 and S4). The high As microbial communities clustered separately than their low As counterparts based on functional gene content (summarized as annotation counts from SEED subsystems) albeit they were also quite different from each other and thus, did not cluster tightly together (Fig. 3A). It was not unexpected to observe clustering at the functional but not taxonomic composition of the sampled communities considering that functional gene content has been shown to be more strongly and faster influenced by environmental conditions than phylogenetic beta-diversity [56]. DESeq2 analysis revealed no statistically significant changes in the relative abundance of genes predicted from metagenomics of different As levels. In contrast, genes associated with As resistance, phage tail fiber, phage photosynthesis, phage capsid, phage cyanophage and CBSS-316279 were all significantly (adjusted p < 0.05) more expressed based on metatranscriptomes by a log₂ fold change of 1–2.3 in high vs low As paddy soils (Fig. 3B).

Relative DNA abundance and expression of As resistance and respiratory genes

Overall, the relative abundance of ars genes, i.e., the cumulative abundance of arsR (arsR1, arsR2, arsR3, arsR4), acr3, arsB, arsC (arsC1, arsC2 and acr2), and arsM, arsI, arsP and arsH genes, were comparable between high and low As paddy soils as revealed by metagenomes (1.2 × 10^− 2 ± 9.7 × 10^− 3 vs 6.1 × 10^− 3 ± 5.3 × 10^− 4 RPKM; Fig. 4). An outlier to this result was sample SKS_AsHig that had three times higher relative abundance of ars genes compared to other samples (Supplementary Figure S2). To robustly assess transcriptome activity while accounting for DNA abundance, e.g., higher DNA abundance would likely result in more mRNA transcripts originating from these cells without necessarily an increase in the rate of gene transcription, the transcriptional level of ars genes were estimated by normalizing metatranscriptomic counts by the relative abundance of same ars genes identified in metagenomes (BLASTx search and filtered for best matches with minimum amino acid identity ≥ 90%, alignment length ≥ 25 amino acid and e-value ≤ 1E-5). More than 80% of the As metabolism genes detected in metatranscriptomes were identified in the metagenomes, representing the majority of the As metabolism genes in the samples. The normalized transcriptional expression levels of ars genes were significantly higher (p < 0.05) in high vs low As paddy soils (2074 ± 701 vs 710 ± 337 RPKM; Fig. 4). Overall, the 13 ars genes responsible for As transcriptional repressor (arsR), As(III) efflux (acr3 and arsB), As(V) reduction (arsC including arsC1, arsC2 and acr2), As(III) methylation (arsM), MMAs(III) demethylation (arsI), MMAs(III) oxidation (arsH) and MMAs(III) efflux (arsP) all consistently showed higher transcriptional activities in high vs low As paddy soils (Fig. 5A). However, for sample HP_AsLow, the arsH, arsR and arsP showed a slightly higher transcriptional activity than one of the samples from high As levels sites (CZ_AsHig or SKS_AsHig), which could be attributed to high soil heterogeneity or biases during library preparation and/or sequencing.

For aioA genes, which encode respiratory As(III) oxidase in the periplasm, a higher relative DNA abundance was revealed, although not statistically significant (p = 0.114), in high vs low As paddy soils (1.0 × 10^− 4 ± 6.2 × 10^− 5 vs 2.9 × 10^− 5 ± 3.2 × 10^− 7 PRKG, respectively; Fig. 4). The transcriptional activity of aioA genes was strongly significantly (p < 0.05) increased in high vs low As paddy soils (117 ± 39 vs 34 ± 31 PRKM, respectively; Fig. 4). Further, there was no significantly higher abundance (p = 0.598) or transcriptional activity of arrA genes (p = 0.33), which is responsible for respiratory As(V) reduction, between high and low As paddy soils (Fig. 4). The relative abundance of arrA genes in metagenomes were comparable between high and low As paddy soils (1.0 × 10^− 4 ± 5.8 × 10^− 5 vs 7.4 × 10^− 5 ± 5.9 × 10^− 5 RPKG, respectively), while their transcriptional activity was slightly higher, but not significant (p = 0.33), in paddy soils with high vs low As levels (39 ± 27 vs 19 ± 16 RPKM, respectively; Fig. 4).

Gene network patterns reveal the co-occurrence patterns of arsR, acr3 and arsM

Network analysis using Sparcc was employed to assess co-occurrence patterns of As-associated genes in the metatranscriptomics datasets (Fig. 5B). The resulting network was strongly structured (modularity 0.78) and consisted of 59 nodes (11 arsR, 23 arsM, 15 acr3, 5 arsI, 4 arrA, and 1 arsC) and 60 edges. 42 positive and 18 negative correlations were found between these As metabolism genes (see details in Supplementary Table S5). The majority of the positive (~ 88% of the total) and negative (~ 89% of the total) correlations are shown in Fig. 5C. As expected, the genes, i.e. arsR, acr3 and arsM, which are on the same ars operon were mostly positively correlated. In addition, the arsI genes were found positively associated with acr3 and arsM genes as well. Further, negative associations among these genes were also observed, but were fewer in number and weaker compared to the positive associations. These negative associations could be attributed to competition pathways of As-associated populations that carry different As gene types for the same As (or another, not As-associated) substrate. For example, the positive association of acr3 with arsM gene was detected between acr3_BAL26929.1 with arsM_AHG89793.1 (cor = 0.834, p = 0.01), while the negative association was detected between acr3_BAL26929.1 with another arsM gene, i.e., arsM_SKSAsHig_gene_id_443484 (cor = -0.902, p = 0.004). Corroboratively, the arsM_AHG89793.1 was also detected to negatively associate with arsM_SKSAsHig_gene_id_443484 (cor = -0.849, p = 0.022). The two arsM genes, i.e. arsM_AHG89793.1 and arsM_SKSAsHig_gene_id_443484 (amino acid identity of 50.7% and query coverage of 77%) could be different gene type that exist in different As-associated populations that sharing potential competing pathways. The limited number and quality of MAGs recovered did not allow us to test this hypothesis further.

Responses to high As level in paddy soils are community-wide rather than taxon-specific

The As-associated microbial populations in high As soils mostly clustered together and separate from their low As counterparts based on the similarities among the expression levels of the ars, aioA and arrA genes (Fig. 6A), except for sample SKS_AsHig, which could be attributed to high soil heterogeneity or the non-significant effect of As levels on arrA genes. In low As soils, the As-associated populations were quite different from each other and thus, did not cluster tightly together (Fig. 6A). In contrast, the As-associated populations did not show any obvious clustering in terms of high As vs low As soils based on metagenomic (DNA) abundances (Supplementary Figure S3). These results are overall consistent with the relative abundance and expression of As metabolism genes, indicating that high As select more for transcriptome response as opposed to gene abundance at DNA level.

To evaluate which clades expressed higher As genes in high As soils, and thus, assess whether the transcriptomic shifts of As metabolism gene described above were due to systematic community-wide response as opposed to a few taxa, the Ars, AioA and ArrA-carrying reads in metagenomes (DNA reads) and metatranscriptomes (RNA reads) were placed on the their representative phylogenetic clade (typically > 95% nucleotide identity within- vs < 95% identity between-clades) and the ratio of RNA/DNA reads assigned at the same reference clade (or gene-based OTU) was compared between samples (Fig. 6B and Supplementary Figure S4-S12). Overall, more than 60% (ranging between 63% and 98%) of the total clades had higher RNA/DNA ratio in high As soils for each As metabolism gene, except for arsB (Supplementary Table S6), suggesting that the high As levels mostly induced a community-wide response. For arsB, due to relative lower abundance of this gene compared to the others, by about 75–375 times in metatranscriptomes compared to acr3 that is also responsible for As(III) efflux, only 17% of the total clades recruited arsB-carrying reads from high As paddy soils and no clade recruited reads from the low As paddy soils (Supplementary Table S6). Moreover, around 54–78% of the total clades for Ars, AioA and ArrA proteins showed a higher ratio of RNA/DNA reads in high vs As paddy soils (Supplementary Table S6), revealing that the As metabolism genes were expressed more strongly in high than low As samples by many distinct clades (taxa). Further, the number of distinct gene-based OTUs (chao1 value of total expected OTUs) was similar or slightly higher in the high than low As paddy soils (1-1.5 fold, in general; Supplementary Table S6) but the transcriptional response of some of these OTUs was much higher (2–3 fold) in high vs low As higher in high than low As samples further corroborating the stronger response in the RNA (transcriptome) than the DNA level to high As concentrations.

Environmental factors affecting the composition of As-associated microbial communities

To investigate the differences between microbial communities in high vs low As soils, five previously characterized samples with high As (As concentration ranged from 34 to 821 mg kg^− 1) affected by a mine field in Hunan, China [57] and eight samples with low As (As concentration ranged from 2 to 16 mg kg^− 1) in Hunan, Zhejiang, Jiangxi and Guangdong, China [25] were included in comparisons with the samples determined by our study. Because these metagenomes had more than two-fold variation in sequencing effort applied, we subsampled the metagenomes to the same level (2 Gbp) in order to avoid false positive signal in detecting features as differentially abundant due to differences in coverage alone [58]. PCoA- or NMDA-based analysis of Mash distances of whole metagenomes showed that samples were mostly separated based on their geographic distances, e.g., the most distinct metagenomes were those that originated form the most distant sites (Fig. 7A and Supplementary Figure S13). The Adonis dissimilarity test corroborated that geographic distribution patterns imposed significant (p < 0.01) changes in microbial community compositions. Notably, As level did not cause any significant changes (p > 0.05) of the corresponding microbial communities.

Based on the samples determined by our study, the concentrations of total phosphorus (TP) and iron (Fe), in contrast to As, were significant factors in driving the variations observed between the microbial communities (Fig. 7B), explaining about 30% of the differences among the sampled total microbial communities. A significant (p < 0.01) interaction was also found between these two factors. Interestingly, the TP concentration also significantly correlated with the relative abundance of different ars genes in the metagenomic datasets (R² = 0.83, p < 0.05). The As concentration was found to significantly correlate with the variance in ars gene transcriptional activities assessed by Bray-Curtis distance between paddy soils (Adonis test, R² = 0.69, p < 0.05) and the transcriptional activities of aioA genes (Pearson’s correlation = 0.85, p < 0.05). In addition to As, the concentration of total nitrogen (TN) (R² = 0.73, p < 0.01) and chromium (R² = 0.72, p < 0.05) also significantly correlated with differences in the transcriptional activities of ars genes, and the concentration of Antimony (Sb) (Pearson’s correlation = 0.83, p < 0.05) significantly correlated with the transcriptional activities of aioA genes (Table 1). No significantly interaction was found between these environmental factors as revealed by PERMANOVA analysis. For the transcriptional activities of arrA genes, no environmental factor was found to affect these activities among the different paddy soils examined.

Table 1

**Correlation of As metabolism genes relative abundances and transcriptional activities with soil physicochemical properties.** For *ars* genes, the Adonis test was used. For *aioA* and *arrA* genes, the Pearson correlation was employed. Significant correlations are denoted by red fonts; * indicates p < 0.05 and ** p < 0.01 confidence levels.
Arsenic Metabolism Genes/Soil physiochemical characteristics		As	Pb	Cd	Cr	Sb	Fe	pH	Eh	C	N	P	S
Relative abundance	ars	0.31	0.02	0.02	0.20	0.92	0.05	0.05	0.03	0.12	0.12	0.83*	0.06
	aioA	0.69	-0.29	-0.29	0.61	0.89*	0.32	0.66	-0.60	0.06	0.58	0.67	0.42
	arrA	0.11	-0.59	-0.63	0.03	0.02	0.15	0.73	-0.77	0.22	0.11	-0.13	0.51
Transcriptional activity	ars	0.69*	0.14	0.54	0.72*	0.23	0.46	0.12	0.08	0.24	0.73**	0.03	0.06
	aioA	0.85*	-0.01	0.00	0.62	0.83*	0.25	0.35	-0.29	-0.06	0.80	0.67	-0.11
	arrA	0.62	-0.03	0.02	0.58	0.22	0.47	-0.31	0.39	0.17	0.64	0.05	-0.60

Arsenic metabolism genes were pervasively presented in the paddy soils either with high (67.3–105 mg kg^− 1) or low (2.5–10.8 mg kg^− 1) As levels, and their relative DNA abundance was mostly comparable between these soils (Fig. 4). Microbial capacity of As metabolism has been shown to be an ancient trait developed in response to the pervasive presence of As in early Earth [3]. Therefore, it is likely that the presence and relative abundances of As metabolism genes in the paddy soils are possibly a legacy effect of an early life characteristic, at least to some extent [4]. Further, the As resistant mechanisms, including the most well studied detoxification As(V) reduction system (ArsC) followed by the efflux of reduced As(III) (ACR3 and ArsB) [1], and the recently identified organoarsenical detoxification systems (ArsM, ArsI, ArsH and ArsP) [10, 59] have been detected in various microbes [60]. The respiratory reduction of As(V) or oxidation of As(III), which is thought to have been developed in the Archaean era [16, 61], have also been described in a broad diversity of microbes [62]. Consistently, the prevalent As metabolism genes have been previously reported in various wetlands with low As levels, such as paddy soils and estuarine sediments (< 15 mg kg^− 1) [21, 25, 26], corroborating the wide prevalence As metabolism genes in our paddy soils.

Notably, the As concentration of ~ 100 mg kg ^− 1 in high As paddy soils did not significantly correlated to a higher relative abundance of As metabolism genes in the resident As-associated microbial community. Previous study of As metabolism genes in mining impacted fields showed a significant (p < 0.05) linear correlation of arsC and aioA gene abundance with the increasing of As concentration from 34 to 821 mg kg ^− 1 [57]. We studied the correlation of ars genes abundances with As concentrations using our own and the previously determined samples by Luo and colleagues and found, consistent with the previous study, a significant correlation (p < 0.05) between As metabolism gene abundance and high As concentration (especially in the 410–812 mg kg ^− 1 range; Supplementary Figure S14). These findings indicated that the 12 times fold higher As concentration in high (average of 92 mg kg ^− 1) vs low (average of 8 mg kg ^− 1) paddy soils of our study was not strong enough to significantly differentiate As genes abundances, and that higher difference in As level and/or high absolute As concentrations (e.g., > 410 mg kg ^− 1) is probably necessary in order to elicit more clear gene content and abundance differences. The possibilities that i) some of the As we measured in our high As soils is not bioavailable (and thus, not selecting for more As-metabolism genes), especially considering that the water soluble or exchangeable (bioavailable) As only accounts for about 0.4% – 24% of the total As in paddy soils [63], or that ii) our measurements are not representative of the average As level in the soils due to soil heterogeneity (e.g., As levels can vary substantially at the millimeter scale) and relative low number of samples analyzed, cannot be excluded at present. To better understand, however, the selection pressure exercised by different As levels, more samples and As gradients ranging from low (i.e. < 15 kg mg^− 1) to higher than (at least) 400 kg mg^− 1 should be included in future studies.

In any case, the increased gradient of As in the paddy soils did upregulated the ars (arsR, acr3, arsB, arsC, arsM, arsI, arsP and arsH) and aioA genes expression level. The co-occurrence patterns of arsR, acr3, and arsM genes revealed by our network analysis (Fig. 5B) further corroborated that the acr3, and arsM genes are usually organized in a single ars operon, which is nearly always under the control of an ArsR transcriptional repressor [12, 13]. Both the ArsR and AioA proteins are induced by As(III), which is the dominant As species in the paddy soils during the flooded period [18, 19], and most likely contributed to the significantly increased transcriptional activities of ars and aioA genes in soils with higher As levels (Fig. 4). The positively association of arsI with acr3 and arsM genes (Fig. 5B and 5C) could be attributed to that the increased transcriptional activity of arsI introduced more As(III) in the soil and consequently, induced the As(III) efflux and As(III) methylation processes, considering that arsI gene is responsible for the demethylation of MMAs(III) to As(III), and acr3 is responsible for As(III) efflux while arsM is responsible for As(III) methylation to MMAs(III) (Fig. 5A). Nevertheless, it cannot be excluded that these gene respond to the same external stimuli (e.g. higher overall As concentrations or the concentration of specific species of As) as opposed to interacting among each other. Related to this, the total As concentration did not significantly correlate with the transcription activity of any As metabolism gene except one arsM gene (arsM_SKSAsHig_gene_id_443484), indicating that the Spracc network analysis might have revealed several true biological interactions among these As genes. Further, the overall transcriptional activity of arrA gene, which is responsible for dissimilarity As(V) reduction [62] was not significantly different between high vs low As paddy soils (Fig. 4). Because ArrA is induced by As(V), it is possible that the dominant As species, i.e. As(III) in paddy soils [18, 19] did not significantly upregulated its transcriptional activity directly. To the best of our knowledge, there are no studies of the transcriptional activity of As metabolism genes in the environment, with the possible exception of a metaproteomic study of an As rich ecosystem, which detected As metabolism proteins, consistent with the findings reported here [64]. Moreover, our study showed the increased transcriptional activity of As metabolism genes was due to several distinct clades (community-wide) in the high As paddy soils (Fig. 6 and Supplementary Figure S4-S12), and this result again corroborated the pervasive presence of As metabolism mechanisms in microbes of the paddy ecosystem.

It’s also noticeable that the relative abundance of aioA genes (DNA level) were about 3.5 times higher in high than low As paddy soils, although not statistically significant (p = 0.114). In contrast, the relative abundance of ars genes at the DNA level was comparable between high and low As levels (Fig. 4). Considering that the As respiratory As(III) oxidation catalyzed by Aio in microbes can couple As(III) oxidation to ATP production, while the Ars is only responsible for As detoxification [14, 60], our finding is consistent with the expectation that energy-production genes (e.g., aioA) will be at higher demand at both the DNA and transcriptional (RNA) levels with increased availability of the substrate, while for the detoxification metabolism genes (e.g., ars) transcriptional response might be adequate with increased substrate and the substrate might not be similarly toxic to all species in order to select for the corresponding genes.

In terms of whole microbial communities, geographic distances, rather than As levels, were found to shape microbial community differences among the paddy soil sites studied here (Fig. 7). The geographic patterns of soil microbial communities have been widely studied previously, and are mostly attributed to varied physicochemical and carbon sources in soil environments [65–67]. It is possible that the As levels in these paddy soils did not possess a strong-enough selective pressure on the whole microbial communities (discussed above). In contrast to As, the concentrations of TP (0.31–0.46 g kg^− 1) and Fe (14.8–27.0 mg kg^− 1), which were about 1.5–1.8 times different between high and low paddy soils, were significant in determining microbial community diversity among our samples (Fig. 7). Previous studies have shown that available organic soil P has an important influence on microbial community composition, shifting both the composition and function of soils microbes [68, 69]. Especially for rhizosphere-associated microbes such as the dominant Proteobacteria in paddy soils [70] their community structure can be altered by soil P availability [71, 72]. The Fe in paddy soils is likely to impact the microbially mediated redox processes [73], and its presence in the wetland environments has been demonstrated to correlate with microbial biogeographic patterns. It is also important to mention that the effect of Fe and TP on As microbes and genes may be interlinked due to the significant (p < 0.01) interaction of P and Fe since P is often the limiting nutrient in many terrestrial ecosystems due to strong sorption on Fe(III) hydroxides [74]. A significant positive correlation between P and Fe in porewaters as well as soil matrix has been previously observed in flooded soils with rice growth [75]. Moreover, their interaction could also affect As behavior in paddy soils [76] because iron oxyhydroxides could sequester As(V) [77] and reduce As bioavailability to microbes.

The dominant microbes in these paddy soils included Proteobacteria, Acidobacteria, Chloroflexi, Nitrospirae, Verrucomicrobia and Bacteroidetes, which have been all documented in paddy soils previously [23, 24]. It should be noted that this result was based on the classification of 16S rRNA gene fragments extracted from metagenomic datasets (Fig. 2). Because of the relatively low coverage as well as the high diversity of microbes in paddy soils, only five metagenomic assembled genomes (MAGs) of adequate quality were obtained from the high As paddy soils after population genome binning (completeness > 80% and contamination < 5%; Supplementary Table S7). These MAGs were assignable to the class of Deltaproteobacteria, Thermomicrobia, Clostridia, Acidobacteriaceae and Solibacteres based on MiGA [78]. Only a couple of the As metabolism genes including arsR, acr3, arsC, arsM and arsP were detected in these MAGs, indicating that these MAGs were capable of As metabolism. However, the remaining As metabolism genes, i.e., arsH, arsI, arsB, aioA and arrA, were not detected in any of the MAGs, possibly due to the incompleteness of these MAGs, and also the relative long length of these genes compared to the other As metabolism genes, especially for AioA and ArrA (800 vs 100–400 amino acids, respectively), which might have prevented their robust assembly into long enough contigs for binning. To better understand, however, higher sequencing efforts should be applied in order to robustly assemble more MAGs, especially those representing the As-transforming taxa, and understand the distributions of As metabolism pathways in different taxa.

Overall, our study revealed that the high As level of ~100 kg mg^-1 increased As metabolism genes activities rather than their abundance or diversity in paddy soils. These findings advance our understanding of how microbes respond to high As levels and the diversity of As metabolism genes in paddy soils. Therefore, the potential functional contributions of the ubiquitous As metabolism genes in the paddy ecosystem, and likely other environments, should be examined not only at the DNA level but also at the transcript and possibly protein levels. The As-associated microbial populations in soils with As could have significant contributions to the biogeochemical cycling of As and the MAGs reported here could facilitate future studies in this area. The work presented here could also be served as a guide for the number of samples to obtain, amount of sequencing to apply, and what bioinformatics analyses to perform for studying the dynamics of the As or other elements metabolism pathways.

As: arsenic; As(V): arsenate; As(III): arsenite; MMAs(III): methylarsenite; MMAs(V): methylarsenate; ars: As resistance gene; arsC: detoxification As(V) reduction gene; arsB and acr3: As(III) efflux gene; arsM: As(III) methylation gene; arsI: MMAs(III) demethylation gene; arsH: MMAs(III) oxidation gene; arsP: MMAs(III) efflux gene; arsR: ars operon transcriptional repressor; aioA: As respiratory oxidation gene; arrA: As respiratory reduction gene; TOC: total organic carbon; TN: total nitrogen; TP: total phosphorus; TS: total sulfur; Pb: lead; Cd: cadmium; Cr: chromium; Sb: antimony; Fe: iron; NCBI: National Center for Biotechnology Information; SRA: Short Read Archive; NMDS: non-metrical multidimensional scaling; PCoA: Principle coordinates analysis; MAG: metagenomic assembled genome;

Acknowledgements

We would thank the China Postdoctoral Science Foundation for their kind support this project. We also appreciate the support from the Brook Byers Institute for Sustainable Systems and the Hightower Chair at the Georgia Institute of Technology.

Author’s contributions

SZ conceived the study. XX collected the paddy soil samples and performed the soil properties analysis. SZ analyzed the metagenomic and metatranscriptomic datasets. SC built the in-house databases of As resistance proteins. SZ and KTK wrote the paper. YZ provided suggestions and revised the paper. All authors read and approved the final manuscript.

Funding

This research was supported, in part, by the China Postdoctoral Science Foundation (No. 212400241) and Major Project of the Ministry of Science and Technology of Jiangxi Province (No. CK201302055). KTK’s research was supported, in part, by the U.S. National Science Foundation (Award No. 1831582 and 1759831).

Availability of data and materials

All relevant data are within the paper and its additional files. Raw sequence data of metagenomics and metatranscriptomics datasets have been deposited in National Center for Biotechnology Information (NCBI)’s Short Read Archive (SRA) database, under the bioproject PRJNA616041.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declaim no conflicts of interest.

Author details

1. School of Ecological and Environmental Sciences, East China Normal University, Shanghai, China. 2. School of Civil and Environmental Engineering, Georgia Institute of Technology, Ford Environmental Science & Technology Building, GA, USA. 3. Department of Biological Sciences, Georgia Institute of Technology, Ford Environmental Science & Technology Building, GA, USA. 4. Key Laboratory of Jiangxi Province for Persistent Pollutants Control and Resources Recycle, Nanchang Hangkong University, Nanchang, Jiangxi, PR China. 5. Department of Isotope Biogeochemistry, Helmholtz Centre for Environmental Research – UFZ, Leipzig, Germany. 6. State Key Laboratory of Urban and Regional Ecology, Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing, PR China. 7. Key Laboratory of Urban Environment and Health, Institute of Urban Environment, Chinese Academy of Sciences, Xiamen, PR China.

Oremland RS, Stolz JF. The ecology of arsenic. Science. 2003;300(5621):939-44.
Zhu Y-G, Yoshinaga M, Zhao F-J, Rosen BP. Earth abides arsenic biotransformations. Annual Review of Earth and Planetary Sciences. 2014;42:443-67.
Saunders JK, Fuchsman CA, McKay C, Rocap G. Complete arsenic-based respiratory cycle in the marine microbial communities of pelagic oxygen-deficient zones. Proceedings of the National Academy of Sciences. 2019;116(20):9925-30.
Sforna MC, Philippot P, Somogyi A, Van Zuilen MA, Medjoubi K, Schoepp-Cothenet B, et al. Evidence for arsenic metabolism and cycling by microorganisms 2.7 billion years ago. Nature Geoscience. 2014;7(11):811-5.
Chen S-C, Sun G-X, Yan Y, Konstantinidis KT, Zhang S-Y, Deng Y, et al. The Great Oxidation Event expanded the genetic repertoire of arsenic metabolism and cycling. Proceedings of the National Academy of Sciences. 2020.
Mukhopadhyay R, Rosen BP. Arsenate reductases in prokaryotes and eukaryotes. Environmental health perspectives. 2002;110(suppl 5):745-8.
Rosen BP. Biochemistry of arsenic detoxification. FEBS letters. 2002;529(1):86-92.
Qin J, Rosen BP, Zhang Y, Wang G, Franke S, Rensing C. Arsenic detoxification and evolution of trimethylarsine gas by a microbial arsenite S-adenosylmethionine methyltransferase. Proceedings of the National Academy of Sciences. 2006;103(7):2075-80.
Chen J, Bhattacharjee H, Rosen BP. ArsH is an organoarsenical oxidase that confers resistance to trivalent forms of the herbicide monosodium methylarsenate and the poultry growth promoter roxarsone. Molecular microbiology. 2015;96(5):1042-52.
Chen J, Madegowda M, Bhattacharjee H, Rosen BP. ArsP: a methylarsenite efflux permease. Molecular microbiology. 2015;98(4):625-35.
Yoshinaga M, Rosen BP. AC⋅ As lyase for degradation of environmental organoarsenical herbicides and animal husbandry growth promoters. Proceedings of the National Academy of Sciences. 2014;111(21):7701-6.
Chen J, Nadar VS, Rosen BP. A novel MAs (III)‐selective ArsR transcriptional repressor. Molecular microbiology. 2017;106(3):469-78.
Qin J, Fu H-L, Ye J, Bencze KZ, Stemmler TL, Rawlings DE, et al. Convergent evolution of a new arsenic binding site in the ArsR/SmtB family of metalloregulators. Journal of Biological Chemistry. 2007;282(47):34346-55.
Slyemi D, Bonnefoy V. How prokaryotes deal with arsenic. Environmental microbiology reports. 2012;4(6):571-86.
Ahmann D, Roberts AL, Krumholz LR, Morel FM. Microbe grows by reducing arsenic. Nature. 1994;371(6500):750-.
Lebrun E, Brugna M, Baymann F, Muller D, Lievremont D, Lett M-C, et al. Arsenite oxidase, an ancient bioenergetic enzyme. Molecular biology and evolution. 2003;20(5):686-93.
Malasarn D, Saltikov C, Campbell K, Santini J, Hering J, Newman D. arrA is a reliable marker for As (V) respiration. Science. 2004;306(5695):455-.
Takahashi Y, Minamikawa R, Hattori KH, Kurishima K, Kihou N, Yuita K. Arsenic behavior in paddy fields during the cycle of flooded and non-flooded periods. Environmental Science & Technology. 2004;38(4):1038-44.
Yamaguchi N, Nakamura T, Dong D, Takahashi Y, Amachi S, Makino T. Arsenic release from flooded paddy soils is influenced by speciation, Eh, pH, and iron dissolution. Chemosphere. 2011;83(7):925-32.
Meharg AA, Zhao F-J. Biogeochemistry of arsenic in paddy environments. In: Arsenic & Rice. Springer; 2012. p. 71-101.
Zhang S-Y, Zhao F-J, Sun G-X, Su J-Q, Yang X-R, Li H, et al. Diversity and abundance of arsenic biotransformation genes in paddy soils from southern China. Environmental science & technology. 2015;49(7):4138-46.
Huang H, Jia Y, Sun G-X, Zhu Y-G. Arsenic speciation and volatilization from flooded paddy soils amended with different organic matters. Environmental science & technology. 2012;46(4):2163-8.
Jia Y, Huang H, Zhong M, Wang F-H, Zhang L-M, Zhu Y-G. Microbial arsenic methylation in soil and rice rhizosphere. Environmental science & technology. 2013;47(7):3141-8.
Jia Y, Huang H, Chen Z, Zhu Y-G. Arsenic uptake by rice is influenced by microbe-mediated arsenic redox changes in the rhizosphere. Environmental science & technology. 2014;48(2):1001-7.
Xiao K-Q, Li L-G, Ma L-P, Zhang S-Y, Bao P, Zhang T, et al. Metagenomic analysis revealed highly diverse microbial arsenic metabolism genes in paddy soils with low-arsenic contents. Environmental pollution. 2016;211:1-8.
Zhang SY, Su JQ, Sun GX, Yang Y, Zhao Y, Ding J, et al. Land scale biogeography of arsenic biotransformation genes in estuarine wetland. Environmental microbiology. 2017;19(6):2468-82.
Ragot SA, Huguenin-Elie O, Kertesz MA, Frossard E, Bünemann EK. Total and active microbial communities and phoD as affected by phosphate depletion and pH in soil. Plant and soil. 2016;408(1-2):15-30.
Thompson I, Bailey M, Ellis R, Maguire N, Meharg A. Response of soil microbial communities to single and multiple doses of an organic pollutant. Soil Biology and Biochemistry. 1998;31(1):95-105.
Cox MP, Peterson DA, Biggs PJ. SolexaQA: At-a-glance quality assessment of Illumina second-generation sequencing data. BMC bioinformatics. 2010;11(1):485.
Su X, Pan W, Song B, Xu J, Ning K. Parallel-META 2.0: enhanced metagenomic data analysis with functional annotation, high performance computing and advanced visualization. PloS one. 2014;9(3).
Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic acids research. 2007;35(21):7188-96.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nature methods. 2010;7(5):335.
Rho M, Tang H, Ye Y. FragGeneScan: predicting genes in short and error-prone reads. Nucleic acids research. 2010;38(20):e191-e.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC bioinformatics. 2009;10(1):421.
Overbeek R, Olson R, Pusch GD, Olsen GJ, Davis JJ, Disz T, et al. The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST). Nucleic acids research. 2014;42(D1):D206-D14.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology. 2014;15(12):550.
Kolde R, Kolde MR. Package ‘pheatmap’. R Package. 2015;1(7).
Wu Y-W, Simmons BA, Singer SW. MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets. Bioinformatics. 2016;32(4):605-7.
Kang DD, Li F, Kirton E, Thomas A, Egan R, An H, et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ. 2019;7:e7359.
Hug LA, Baker BJ, Anantharaman K, Brown CT, Probst AJ, Castelle CJ, et al. A new view of the tree of life. Nature microbiology. 2016;1(5):1-6.
Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, et al. The Pfam protein families database. Nucleic acids research. 2002;30(1):276-80.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30(9):1312-3.
Peng Y, Leung HC, Yiu S-M, Chin FY: IDBA–a practical iterative de Bruijn graph de novo assembler. In: Annual international conference on research in computational molecular biology. Springer; 2010: 426-40.
Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular biology and evolution. 2013;30(4):772-80.
Nayfach S, Pollard KS. Average genome size estimation improves comparative metagenomics and sheds light on the functional ecology of the human microbiome. Genome biology. 2015;16(1):51.
Wagner H. Vegan: community ecology package. R package. 2015.
Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic acids research. 2016;44(W1):W242-W5.
Rodriguez-R LM, Konstantinidis KT: The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes. In.: PeerJ Preprints; 2016.
Friedman J, Alm EJ. Inferring correlation networks from genomic survey data. PLoS Comput Biol. 2012;8(9):e1002687.
Li B, Yang Y, Ma L, Ju F, Guo F, Tiedje JM, et al. Metagenomic and network analysis reveal wide distribution and co-occurrence of environmental antibiotic resistance genes. ISME J. 2015;9(11):2490-502; doi: 10.1038/ismej.2015.59.
Bastian MH, Sebastien; Jacomy, Mathieu. Gephi: an open source software for exploring and manipulating net- works. Proceedings of the Third International ICWSM Conference 2009:361-2.
Rodriguez-r LM, Konstantinidis KT. Nonpareil: a redundancy-based approach to assess the level of coverage in metagenomic datasets. Bioinformatics. 2014;30(5):629-35.
Rodriguez-R LM, Gunturu S, Tiedje JM, Cole JR, Konstantinidis KT. Nonpareil 3: fast estimation of metagenomic coverage and sequence diversity. MSystems. 2018;3(3):e00039-18.
Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, et al. Mash: fast genome and metagenome distance estimation using MinHash. Genome biology. 2016;17(1):132.
Wickham H. ggplot2: elegant graphics for data analysis. Springer; 2016.
Kuang J, Huang L, He Z, Chen L, Hua Z, Jia P, et al. Predicting taxonomic and functional structure of microbial communities in acid mine drainage. The ISME journal. 2016;10(6):1527-39.
Luo J, Bai Y, Liang J, Qu J. Metagenomic approach reveals variation of microbes with arsenic and antimony metabolism genes from highly contaminated soil. PLoS One. 2014;9(10).
Rodriguez-r LM, Konstantinidis KT. Estimating coverage in metagenomic data sets and why it matters. The ISME journal. 2014;8(11):2349-51.
Yang H-C, Rosen BP. New mechanisms of bacterial arsenic resistance. Biomedical journal. 2016;39(1):5-13.
Bhattacharjee H, Rosen BP. Arsenic metabolism in prokaryotic and eukaryotic microbes. In: Molecular microbiology of heavy metals. Springer; 2007. p. 371-406.
Duval S, Ducluzeau A-L, Nitschke W, Schoepp-Cothenet B. Enzyme phylogenies as markers for the oxidation state of the environment: the case of respiratory arsenate reductase and related enzymes. BMC evolutionary biology. 2008;8(1):206.
Silver S, Phung LT. Genes and enzymes involved in bacterial oxidation and reduction of inorganic arsenic. Appl Environ Microbiol. 2005;71(2):599-608.
Bhattacharyya P, Tripathy S, Kim K, Kim S-H. Arsenic fractions and enzyme activities in arsenic-contaminated soils by groundwater irrigation in West Bengal. Ecotoxicology and environmental safety. 2008;71(1):149-56.
Bertin PN, Heinrich-Salmeron A, Pelletier E, Goulhen-Chollet F, Arsène-Ploetze F, Gallien S, et al. Metabolic diversity among main microorganisms inside an arsenic-rich ecosystem revealed by meta-and proteo-genomics. The ISME journal. 2011;5(11):1735-47.
Ma B, Wang H, Dsouza M, Lou J, He Y, Dai Z, et al. Geographic patterns of co-occurrence network topological features for soil microbiota at continental scale in eastern China. The ISME journal. 2016;10(8):1891-901.
Plassart P, Prévost-Bouré NC, Uroz S, Dequiedt S, Stone D, Creamer R, et al. Soil parameters, land use, and geographical distance drive soil bacterial communities along a European transect. Scientific reports. 2019;9(1):1-17.
Cao H, Chen R, Wang L, Jiang L, Yang F, Zheng S, et al. Soil pH, total phosphorus, climate and distance are the major factors influencing microbial activity at a regional spatial scale. Scientific reports. 2016;6(1):1-10.
DeForest JL, Scott LG. Available organic soil phosphorus has an important influence on microbial community composition. Soil Science Society of America Journal. 2010;74(6):2059-66.
Heuck C, Weig A, Spohn M. Soil microbial biomass C: N: P stoichiometry and microbial use of organic phosphorus. Soil Biology and Biochemistry. 2015;85:119-29.
Somenahally AC, Hollister EB, Loeppert RH, Yan W, Gentry TJ. Microbial communities in rice rhizosphere altered by intermittent and continuous flooding in fields with long-term arsenic application. Soil Biology and Biochemistry. 2011;43(6):1220-8.
Castrillo G, Teixeira PJPL, Paredes SH, Law TF, de Lorenzo L, Feltcher ME, et al. Root microbiota drive direct integration of phosphate stress and immunity. Nature. 2017;543(7646):513-8.
Jackson BP, Miller W. Effectiveness of phosphate and hydroxide for desorption of arsenic and selenium species from iron oxides. Soil Science Society of America Journal. 2000;64(5):1616-22.
Kögel-Knabner I, Amelung W, Cao Z, Fiedler S, Frenzel P, Jahn R, et al. Biogeochemistry of paddy soils. Geoderma. 2010;157(1-2):1-14.
Borch T, Fendorf S. Phosphate interactions with iron (hydr) oxides: Mineralization pathways and phosphorus retention upon bioreduction. Developments in Earth and Environmental Sciences. 2007;7:321-48.
Wang Y, Yuan J-H, Chen H, Zhao X, Wang D, Wang S-Q, et al. Small-scale interaction of iron and phosphorus in flooded soils with rice growth. Science of the Total Environment. 2019;669:911-9.
Ji Y, Luo W, Lu G, Fan C, Tao X, Ye H, et al. Effect of phosphate on amorphous iron mineral generation and arsenic behavior in paddy soils. Science of The Total Environment. 2019;657:644-56.
Liu W, Zhu Y, Hu Y, Williams P, Gault A, Meharg AA, et al. Arsenic sequestration in iron plaque, its accumulation and speciation in mature rice plants (Oryza sativa L.). Environmental Science & Technology. 2006;40(18):5730-6.
Rodriguez-R LM, Gunturu S, Harvey WT, Rosselló-Mora R, Tiedje JM, Cole JR, et al. The Microbial Genomes Atlas (MiGA) webserver: taxonomic and gene diversity analysis of Archaea and Bacteria at the whole genome level. Nucleic acids research. 2018;46(W1):W282-W8.

Download PDF

Journal Publication

published 28 Sep, 2021

Read the published version in Applied and Environmental Microbiology →

Version 1

posted 13 Jan, 2021

You are reading this older preprint version

Read the latest preprint version →

High Arsenic Levels Increase Activity Rather Than Diversity or Abundance of Arsenic Metabolism Genes in Paddy Soils

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Methods

Sites description and sample collection

Soil DNA and RNA isolation, library preparation and sequencing

Taxonomic classification, functional annotation and population genome binning of soil metagenomes and metatranscriptomes

Assessment of the abundance and transcriptional activity of arsenic metabolism genes

Building an As metabolism protein reference database

Estimation of As metabolism genes abundance and transcriptional activity

Phylogenetic placement of As metabolism gene carrying reads

Network analysis

Alpha and beta community diversity estimation

Results

Soil physicochemical characteristics

Broad microbial community diversity indices

Paddy soil community taxonomic and functional composition

Relative DNA abundance and expression of As resistance and respiratory genes

Responses to high As level in paddy soils are community-wide rather than taxon-specific

Environmental factors affecting the composition of As-associated microbial communities

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1