A global survey of arsenic-related genes in soil microbiomes

Background Environmental resistomes include transferable microbial genes. One important resistome component is resistance to arsenic, a ubiquitous and toxic metalloid that can have negative and chronic consequences for human and animal health. The distribution of arsenic resistance and metabolism genes in the environment is not well understood. However, microbial communities and their resistomes mediate key transformations of arsenic that are expected to impact both biogeochemistry and local toxicity. Results We examined the phylogenetic diversity, genomic location (chromosome or plasmid), and biogeography of arsenic resistance and metabolism genes in 922 soil genomes and 38 metagenomes. To do so, we developed a bioinformatic toolkit that includes BLAST databases, hidden Markov models and resources for gene-targeted assembly of nine arsenic resistance and metabolism genes: acr3, aioA, arsB, arsC (grx), arsC (trx), arsD, arsM, arrA, and arxA. Though arsenic-related genes were common, they were not universally detected, contradicting the common conjecture that all organisms have them. From major clades of arsenic-related genes, we inferred their potential for horizontal and vertical transfer. Different types and proportions of genes were detected across soils, suggesting microbial community composition will, in part, determine local arsenic toxicity and biogeochemistry. While arsenic-related genes were globally distributed, particular sequence variants were highly endemic (e.g., acr3), suggesting dispersal limitation. The gene encoding arsenic methylase arsM was unexpectedly abundant in soil metagenomes (median 48%), suggesting that it plays a prominent role in global arsenic biogeochemistry. Conclusions Our analysis advances understanding of arsenic resistance, metabolism, and biogeochemistry, and our approach provides a roadmap for the ecological investigation of environmental resistomes. Electronic supplementary material The online version of this article (10.1186/s12915-019-0661-5) contains supplementary material, which is available to authorized users.


Background
Microbial communities drive global biogeochemical cycles through diverse functions. The biogeography of functional genes can help to predict and manage the influence of microbial communities on biogeochemical cycling [1]. These trait-based analyses require that the functional genes are well-characterized from both evolutionary and genetic perspectives [2]. The arsenic resistance and metabolism genes exemplify a suite of well-characterized functional genes that have consequences for biogeochemistry. Arsenic is a toxic metalloid that, upon exposure, can have negative effects for all life, including humans, livestock, and microorganisms. The toxicity and mobility of arsenic depends, in part, on its oxidation state: the trivalent arsenite is more mobile and more toxic than the pentavalent arsenate [3]. The toxicity of methylated arsenic species varies with oxidation state and number of methyl groups (monomethyl, dimethyl, trimethyl). Pentavalent methylarsenicals are progressively less toxic than inorganic arsenate, while trivalent methylarsenicals are progressively more toxic than inorganic arsenite with the exception of trimethylarsine which is the least toxic arsenic species [4,5]. Additionally, volatilization of arsenic can occur through methylation [6], which has varied impacts. Methylated forms of arsenic can be released to new areas through air [7], captured during bioremediation [8], or accumulate in crops such as rice [9]. Microbial transformations of arsenic can have consequences for arsenic speciation and methylation; therefore, they impact arsenic ecotoxicity and the fate of arsenic in the environment.
Arsenic biogeochemical cycling by microbial communities is both an ancient [10,11] and a contemporary [3,12] phenomenon. Changes to the methylation or oxidation state of arsenic alter biogeochemical cycling of arsenic, and microbes have evolved a variety of mechanisms to carry out these functions. Arsenic-related genes are generally separated into two categories: resistance and metabolism [13]. Arsenic resistance, or detoxification, is encoded by the ars operon [14]. The ars operon protects the cell from arsenic but does not detoxify arsenic itself in the environment. This operon includes arsenite efflux (ArsB, Acr3) which is potentially precluded by cytoplasmic arsenate reduction with either glutaredoxin (ArsC (grx)) or thioredoxin (ArsC (trx)) [14]. Arsenic metabolisms include methylation (ArsM), oxidation (AioAB, ArxAB), and dissimilatory reduction (ArrAB) [13]. While these genetic determinants of arsenic detoxification and metabolism are well-characterized, the full scope of arsenic detoxification and metabolism gene distribution, diversity, and interspecies transfer is unknown [15][16][17].
Microbial arsenic resistance is reportedly widespread in the environment. Arsenic-resistant organisms have been found in sites with low arsenic concentrations (< 7 ppm) [18,19], and it has been speculated that nearly all organisms have arsenic resistance genes [20]. While the number of identified microorganisms with arsenic resistance genes continues to grow [13], the number of microorganisms without arsenic resistance genes is unclear. Furthermore, though the complete arsenic biogeochemical cycle has been detected in the environment [10], the relative contributions of genes encoding detoxification and metabolism remain unknown [11]. A global, biogeographic perspective of environmental arsenic-related genes would improve understanding of their ecology. This information would expand foundational knowledge of arsenic detoxification and metabolism, including local and global abundances, gene diversity, dispersal across different environments, and representations over the microbial tree of life.
Knowledge gaps concerning the diversity of microbial arsenic-related genes are driven, in part, by numerous inconsistencies in nomenclature and detection methods. Though public microbial metagenome and genome data continue to surge, there are several practical hurdles to achieving a robust, global assessment of microbial arsenic-related genes from this wealth of data. First, tools to detect these genes rely on imperfect annotation [15] and widely vary in nomenclature [21]. Next, the use of different reference databases [12,[22][23][24][25] and normalization techniques [25,26] complicates comparisons between studies. To overcome these hurdles, we developed an open-access toolkit to examine arsenic resistance and metabolism genes in microbial sequence datasets. This toolkit allowed us probe genomic and metagenomic datasets simultaneously to investigate arsenic-related genes in soil microbiomes. We first asked whether arsenic-related genes are universal in soil-associated microorganisms. Next, we tested the hypothesis that genes encoding arsenic detoxification are more abundant than those encoding arsenic metabolism. We also tested the hypothesis that arsenic resistance genes with redundant function (i.e., acr3 and arsB; arsC (grx) and arsC (trx)) would have complementary environmental abundances. Third, we asked whether estimations of arsenic-related gene abundance are biased by cultivation efforts, as cultivation is often a research emphasis because cultivable, arsenic-resistant microorganisms can be used in bioremediation [17]. Finally, we tested the hypothesis that sequence variants of arsenic-related genes are endemic, not cosmopolitan.

A bioinformatic toolkit for detecting and quantifying arsenic-related genes
We developed a toolkit to improve investigations of microbial arsenic-related genes (Fig. 1a, b) [14,[31][32][33][34][35]. We selected these nine genes because they are markers of arsenic detoxification and metabolism [21,25] and because their genetic underpinnings are well established. Seed sequences (high-quality and full-length sequences) for each gene of interest were collected and used to construct BLAST databases [30], functional gene (FunGene) databases [27], hidden Markov models (HMMs [36]), and gene resources for gene-targeted assembly (Xander [28]) (Fig. 1a). Altogether, this toolkit relies on consistent references and nomenclature and can search both amino acid and nucleotide sequence data.
To demonstrate the utility of our toolkit, we performed an analysis of arsenic-related genes in soil-associated genomes and metagenomes. We used HMMs for marker genes for arsenic detoxification and metabolism to search RefSoil+ genomes, a set of complete chromosomes and plasmids from cultivable soil microorganisms [37]. Additionally, we used a gene-targeted assembler [28] to test 38 public soil metagenomes from Brazil, Canada, Malaysia, Russia, and the USA for arsenic resistance and metabolism genes (Additional file 1). Ultimately, these data serve as a broad baseline of arsenic detoxification and metabolism genes in soil.

Phylogenetic distributions and genomic locations of arsenic-related genes
We asked whether arsenic resistance and metabolism genes were universal in RefSoil+ organisms [37]. Of the 922 RefSoil+ genomes spanning 25 phyla ( Fig. 2b; Additional file 2), 14.3% (132 genomes) did not contain any tested arsenic-related genes. Of the 25 phyla in RefSoil+, two phyla (Chlamydiae and Crenarchaeota) did not have any of these genes. These phyla, however, had few RefSoil+ representatives (three and nine, respectively), so other members of these phyla may have arsenic detoxification and metabolism genes. Supporting this hypothesis, a Crenarchaeota isolate was previously reported to oxidize arsenic [38]. Nonetheless, these data suggest that arsenic-related genes are widespread, but not universal, even among cultivable soil organisms (Fig. 2).

Phylogenetic diversity of arsenic-related genes: insights into vertical and horizontal transfer Arsenite efflux pumps
We examined the phylogenetic diversity of distinct genes encoding arsenite efflux pumps, acr3 and arsB, for soil-associated microorganisms (Fig. 3, Additional files 3 and 4). Gene acr3 is separated into two clades: acr3(1) and acr3(2) [40]. Clade acr3(1) is typically composed of Proteobacterial sequences while acr3(2) is typically composed of Firmicutes and Actinobacterial sequences [21,40,41]. Though RefSoil+ genomes were mostly composed of acr3(2) sequences from Proteobacteria ( Fig. 3a; Additional file 3), we observed greater taxonomic diversity observed than previously reported for this clade [21,40,41]. Surprisingly, there were deep branches in acr3(2) that belonged to Bacteroidetes, Euryarchaeota, Firmicutes, Fusobacteria, and Verrucomicrobia. Similarly, acr3(1) contained closely related acr3 sequences present in a diverse array of phyla (10 out of 25). Both clades had sequences present on plasmids (6.1%). Plasmid-borne arsB sequences were only present in Proteobacteria and Deinococcus-Thermus strains ( Fig.  3b; Additional file 4). Sequences from Actinobacteria, Proteobacteria, and Firmicutes were each present in two distinct phylogenetic groups, and previous studies also observed separation of arsB sequences based on phylum [40,41]. Interestingly, our genome-centric analysis revealed that microorganisms with multiple copies of arsB did not harbor identical copies. For example, seven Bacillus subtilis subsp. subtilis strains had two copies of arsB, with one from each of the two clades (Additional file 4). Arsenic resistance and metabolism gene toolkit schematic. a Seed sequences for nine arsenic resistance genes were used to construct an arsenic resistance gene database with existing tools [27][28][29][30]. Lines indicate interdependence between modules. b Table of arsenic resistance and metabolism genes included in the toolkit. The toolkit is freely available on GitHub: https://github.com/ShadeLab/PAPER_Dunivin_meta_arsenic A B Fig. 2 Arsenic resistance and metabolism genes in RefSoil+ organisms. a Maximum likelihood tree of 16S rRNA genes in RefSoil+ organisms. Bootstrap support > 50 is shown with black circles. Tree branches and the first ring are colored by organism taxonomy. Each node is annotated with arsenic resistance genotype where color indicates the gene. Filled boxes indicate gene presence on chromosome, and open boxes indicate gene presence on plasmid. b Proportion of RefSoil+ organisms and organisms containing arsenic resistance genes are colored by the taxonomy of the organism containing the gene. "None" refers to the number of genomes that do not test positive for any of the nine arsenic resistance genes analyzed. Note the difference between y-axes

Cytoplasmic arsenate reductases
Cytoplasmic arsenate reductase (ArsC (trx)) was phylogenetically widespread in RefSoil+ microorganisms ( Fig. 4a; Additional file 5). While some arsC (trx) sequences were plasmid-borne, the majority were chromosomally encoded. Similarly, plasmid-encoded arsC (grx) made up 4.6% of RefSoil+ hits ( Fig. 4b; Additional file 6). Notably, several Proteobacteria strains have multiple copies of arsC (grx) with distinct sequences. It is possible that this is the result of an early gene duplication event or HGT of a second arsC (grx).
Arsenic metabolism genes aioA, arrA, and arxA were phylogenetically conserved (Fig. 6). Genes encoding arsenite oxidases aioA and arxA were restricted to Proteobacteria. aioA sequences clustered into two clades based on class-level taxonomy: all Alphaproteobacteria sequences cluster separately from Gamma-and Betaproteobacteria sequences. The gene encoding dissimilatory arsenate reduction arrA was also phylogenetically conserved in RefSoil+ strains, with strains from Proteobacteria clustering separate from Firmicutes (Fig. 6).

Cultivation bias and environmental distributions of arsenic-related genes
To gain a cultivation-dependent perspective of the abundances of arsenic-related genes in soils, we used inferred environmental abundances of RefSoil microorganisms [42,43]. The environmental abundance of RefSoil microorganisms, which are cultivable, soil-associated microorganisms, was previously estimated by comparing 16S rRNA gene sequences in RefSoil with those in soil metagenomes [42]. We used this estimated abundance of cultivable microorganisms along with arsenic-related gene information from this study ( Fig. 2) to estimate the environmental abundances of arsenic-related genes from the cultivated bacteria. Arsenic metabolism genes (aioA, arrA, arsM, arxA) were predicted to be less common in the environment compared with arsenic detoxification genes (acr3, arsB, arsC (grx), arsC (trx), and arsD) ( Fig. 7a; Mann-Whitney U test p < 0.01). Despite similar distributions of acr3 and arsB in RefSoil+ (Fig. 2b), acr3 was more abundant in most soil orders ( Fig. 7a; Mann-Whitney U test p < 0.05). For genes encoding cytoplasmic arsenate reductases, arsC (grx) was more abundant than arsC (trx) (Mann-Whitney U test p < 0.01).
To gain a cultivation-independent perspective of the abundances of arsenic-related genes, we examined their normalized abundance from soil metagenomes (Fig. 7b). An undetected gene does not confirm absence, so we present a conservative estimate that only includes metagenomes testing positive for a gene. Arsenic detoxification genes (acr3, arsB, arsC (grx), arsC (trx), and arsD) were more abundant than arsenic metabolism genes (aioA, arrA, arsM, and arxA) (Mann-Whitney U test p < 0.01; Fig. 7b). Genes encoding arsenite efflux pumps differed in their abundance with acr3 being more abundant than arsB (Mann-Whitney U test p < 0.01). We also observed differences in cytoplasmic arsenate reductases: arsC (grx) was more abundant than arsC (trx) (Mann-Whitney U test p < 0.01).
We explored cultivation bias of arsenic-related genes with a case study comparing cultivation-dependent (lawn growth on the standard medium TSA50) and cultivation-independent communities from the same soil. Genes in the ars operon (acr3, arsB, arsD, and arsC (trx)) were elevated in the cultivation-dependent metagenome (Fig. 7c). Additionally, arsenic metabolism genes were not detected (aioA, arrA, arxA) or in low abundance (arsM) in the cultivation-dependent sample; however, all four of these arsenic metabolism genes were detected in the cultivation-independent sample. Though this is a single-case study of cultivation-dependent and cultivation-independent methods, these results recapitulate the general discrepancies between RefSoil+ genomes and soil metagenomes (Fig. 7b). This bias has important implications for studies focusing on arsenic bioremediation because cultivation-dependent studies could misestimate the potential of microbiomes for arsenic detoxification and metabolism in situ.

Arsenic-related gene endemism
Arsenic-related genes are globally distributed, but their biogeography is poorly understood. Broadly, arsenic-related genes had comparable abundance among different soils (Fig. 7a, b). The relative distributions of distinct arsenic detoxification and metabolism mechanisms in one site, however, are relevant for predicting the impact of microbial communities on the fate of arsenic. To understand site-specific distributions, we explored soil metagenomes from Brazil, Canada, Malaysia, Russia, and the USA (Additional file 1). These 16 sites had differences in community membership (Additional file 9) and arsenic-related gene content (Fig. 8a). Geographic location was not predictive of arsenic-related gene content (Mantel's r = 0.03493; p > 0.05). Soils had different distributions of arsenic-related genes and therefore differed in their potential impact on the biogeochemical cycling of arsenic. While arsC (grx) and arsM dominated most samples, their relative proportions varied greatly (Fig. 8a). RefSoil+ data suggests that arsM can be found in Verrucomicrobia (100%, n = 2), which is of particular importance for soil metagenomes since Verrucomicrobia are often underestimated with cultivation-dependent methods [44]. The mangrove sample had the most even proportions of arsenic-related genes (Fig. 8a). This distribution was driven by a high abundance of arsC (trx) and arrA. We further examined the arsenic resistance gene abundance at individual sites. We did not include arr and arx in this analysis due to limited available data. For each gene, the abundance varied greatly, but replicates within one site had similar abundances (Fig. 8b). The majority of arsenic-related gene sequences (99.3%) were endemic and only found in one to two sites, but 24 sequences were detected in three or more sites (Fig. 8c; Additional file 10). The majority (70.8%) of cosmopolitan sequences belonged to arsC (grx). This analysis suggests that arsenic-related genes acr3, arsB, arsC (trx), arsD, arsM, and aioA are generally endemic.

Phylogenetic diversity and distribution of arsenic-related genes
It has been conjectured that nearly all organisms have arsenic resistance genes [20], and though this assumption has propagated in the literature, it had never been explicitly quantified. Our data suggest that arsenic detoxification and metabolism genes are ubiquitous, but not universal in RefSoil+ microorganisms (Fig. 2). It is possible for these 132 organisms to have untested or novel arsenic-related genes; nonetheless, these nine well-characterized genes were not universally detected. Additionally, phylogeny was predictive of the presence of acr3, arsB, arsC (grc), arsC (trx), and arsM. This correlation suggests that taxonomy is predictive of arsenic genotype despite documented potential for HGT [19,40,48,52,53]. This result could be explained by ancient rather than contemporary HGT, as seen with arsM [53] and arsC (grx) [48]. Therefore, we next assessed evidence for HGT by examining the phylogenetic congruence and genomic location (e.g., chromosome or plasmid) of arsenic-related gene sequences.
Horizontal transfer of arsenic-related genes has been well documented [19,40,48,[52][53][54][55] and is an important consideration for understanding the propagation and taxonomic identity of arsenic-related genes. We examined the phylogenetic diversity of arsenic-related genes in RefSoil+ microorganisms, including plasmids and chromosomes, and compared them with the 16S rRNA gene taxonomy.

Efflux pumps
While known acr3 sequences separate into two clades [21,40,41], plasmid-borne acr3 sequences were present across clades, suggesting a potential for transfer across unrelated taxa. Therefore, studies assigning taxonomy to acr3 in the absence of host information should consider the clade precisely and proceed with caution. Despite their functional redundancy as arsenite efflux pumps, acr3 and arsB have very distinctive diversity. As compared with acr3, arsB was less diverse and more phylogenetically conserved ( Fig. 3b; Additional file 4). This observation is in agreement with previous reports comparing the diversity of arsB to acr3 [40,41]. Multiple, phylogenetically distinct copies of arsB were present in some RefSoil+ organisms, which could be due to an early gene duplication and subsequent diversification or to an early transfer event. Therefore, despite relatively lower sequence variation, this arsB phylogeny suggests an interesting evolutionary history that could be investigated further.
(See figure on previous page.) Fig. 7 Comparison of arsenic resistance and metabolism gene abundance between cultivation-dependent and cultivation-independent methods. a Mean normalized abundance of arsenic-related genes based on RefSoil microorganisms abundance estimated from corresponding 16S rRNA gene abundance in Earth Microbiome Project datasets. Points are colored by soil order. b Normalized abundance of arsenic resistance genes in RefSoil+ and 38 metagenomes. Metagenome abundance was normalized to rplB, and RefSoil+ normalized abundance was calculated using the number of RefSoil+ genomes. Only metagenomes with an arsenic resistance gene detected are shown, and the total number of datasets (including RefSoil+) is shown in parentheses. c rplB-normalized abundance of arsenic resistance genes in cultivation-dependent and cultivationindependent metagenomes from the same soil sample Cytoplasmic arsenate reductases arsC (trx) was predominantly found on RefSoil+ chromosomes, not plasmids, suggesting vertical transfer of arsC (trx) is common. arsC (trx) was present in both Bacteria and Archaea, and sequences from the two domains formed two distinct clades. arsC (trx) sequences that cluster separately from Bacterial-arsC (trx) sequences have been documented in Thermococci, Archaeoglobi, Thermoplasmata, and Halobacteria [56]. Together, this distribution supports an early evolutionary origin for arsC (trx). Thus, arsC (trx) appears to be an evolutionarily old enzyme that is phylogenetically conserved despite its presence on plasmids and potential for HGT. Plasmid-encoded arsC (grx) were also observed in RefSoil+ microorganisms, highlighting a contemporary potential for HGT that has been documented in soil [48]. Thus, both genes encoding cytoplasmic arsenate reductases were more common on chromosomes.

Arsenic metabolisms
The evolutionary history of the gene encoding arsenite S-adenosylmethionine methyltransferase, arsM, was recently investigated [52,53]. Both studies independently determined that arsM evolved billions of years ago and was subject to HGT [52,53]. In this work, arsM sequences from Euryarchaeota were dispersed throughout the arsM phylogeny, supporting the potential for inter-kingdom transfer events that were recently suggested [52,53]. Very few RefSoil+ organisms had arsenic metabolism genes aioA, arrA, or arxA, which limits phylogenetic analysis. Nonetheless, they were mostly found in Proteobacteria, which is in agreement with previous work [13].

Cultivation bias and environmental distributions of arsenic-related genes
Cultivation-based assessments of arsenic-related gene content are important since cultivable strains are often favored for bioremediation [57]. We estimated distributions of arsenic-related genes in cultivable microorganisms from soils and found a greater abundance of arsenic detoxification genes acr3, arsB, and arsC (trx) (Fig. 7a). A previous study also reported an abundance of acr3 over arsB in cultivable microoganisms from forest soils and attributed this to the greater phylogenetic distribution of acr3 compared with arsB [41].
Additionally, they found that arsC (grx) was more abundant than arsC (trx) in cultivated microorganisms from these soils. It has been posited in cultivation-independent studies that arsC (trx) is more efficient than arsC (grx) and that high local arsenic concentrations result in a relatively greater abundance of arsC (trx) [21,58]. Our cultivation-dependent abundances suggest that acr3 and arsC (grx), rather than arsB and arsC (trx), predominantly comprise the arsenic detoxification pathway in soils.
To assess arsenic-related gene content without cultivation bias, we examined arsenic-related genes in soil metagenomes. As predicted by cultivable organisms, arsenic metabolism genes (aioA, arrA, arxA) were generally in low abundance while acr3 and arsC (grx) were in high abundance. Estimates of genes encoding arsenic detoxification (acr3, arsB, arsD, arsC (grx), arsC (trx)) were considerably lower in these cultivation-independent samples. This result could be due, in part, to the large number of RefSoil+ microorganisms with multiple copies of these genes (Additional file 8). Cultivation-independent genomes (e.g., single-cell-amplified genomes and metagenome-assembled genomes) could provide greater context about the environmental distributions of copy numbers of arsenic-related genes.
Notably, arsM was abundant in soil (median 48%), which greatly exceeds cultivation-dependent estimations, and in a case study of cultivation-dependent and cultivation-independent techniques, arsM was more abundant in the cultivation-independent sample (Fig. 7c). Due to the early phylogenetic origins of arsM and its independent functionality [53], this abundance of arsM in soil metagenomes is not unexpected. arsM is typically studied in paddy soils [6,59,60], but metagenomes in this study suggest it is an important component of the arsenic biogeochemical cycle in a variety of soils.

Arsenic-related gene endemism
We examined the relative abundance of arsenic-related genes in soil metagenomes and observed differences in genetic potential for arsenic transformation that could impact biogeochemical cycling (Fig. 8a). Notably, the mangrove sample had the most even proportions of arsenic-related genes. While the arsenic concentrations in this sample are unknown, mangroves are considered sources and sinks for arsenic [61][62][63]. This could explain (See figure on previous page.) Fig. 8 Arsenic resistance and metabolism gene biogeography. a Relative abundance of arsenic resistance genes in soil metagenomes. b Rank rplB-normalized abundance of arsenic-related genes in soil metagenomes. Sites are ordered by rank mean abundance. Note the differences in yaxes. c Abundance-occurrence plots of arsenic-related gene sequences clustered at 90% amino acid identity. Number of samples included are as follows: Brazilian forest n = 3, California grassland n = 2, Centralia active n = 7, Centralia recovered n = 5, Centralia reference n = 1, Disney preserve n = 2, Illinois soybean n = 2, Illinois switchgrass n = 1, Iowa agricultural n = 2, Iowa corn n = 2, Iowa prairie n = 3, Mangrove n = 2, Minnesota grassland n = 2, Permafrost Canada n = 2, Permafrost Russia n = 1, and Wyoming soil n = 1 the greater abundance of arsC (trx), which is hypothesized to be more abundant in high arsenic sites [21,58]. Additionally, arrA encodes a dissimilatory arsenate reductase that functions in an anaerobic environment [34], so its greater abundance in sediment is expected. Soil geochemical data was not available for all metagenomes examined in this work, so direct comparisons of arsenic-related gene content and soil geochemistry were not possible. This highlights the importance and utility of depositing geochemical data with DNA sequences. Future work, however, could further examine relationships between arsenic resistance genes and soil geochemical data, including arsenic concentration and redox potential.
We also measured whether arsenic-related gene sequence variants were endemic or cosmopolitan in soil metagenomes (Fig. 8c). We found that genes acr3, arsB, arsC (trx), arsD, arsM, and aioA were generally endemic, suggesting regional dispersal limitation. Only one aioA and three acr3 sequences were detected in multiple sites. This supports a previous finding that acr3 and aioA from the acid mine drainage in Carnoulès were endemic [64]. Conversely, arsC (grx) was cosmopolitan which could suggest genetic migration via HGT or vertical transfer and a limited gene diversification. Both are plausible since arsC (grx) was common in RefSoil+ plasmids and had low phylogenetic diversity ( Fig. 4b; Additional file 6).

Conclusions
We developed a bioinformatic toolkit for detecting arsenic detoxification and metabolism genes in microbial sequence data and applied it to analyze the genomes and metagenomes from soil microorganisms. This toolkit informs hypotheses about the evolutionary histories of these genes (including potential for vertical and horizontal transfers) and how community ecology in situ may influence their prevalence and distribution. This study reports the phylogenetic diversity, genomic locations, and biogeography of arsenic-related genes in soils, integrating information from different 'omics datasets and resources to provide a broad synthesis. The toolkit and the synthesis presented here can catalyze future work to understand the ecology and evolution of microbial arsenic biogeochemistry. Furthermore, the toolkit acts as a framework for similar studies of other functional genes of interest.

Gene selection and functional gene (FunGene) database construction
Marker genes can be used to estimate their potential to influence the arsenic biogeochemical cycle [21,25], so we selected nine well-characterized genes: acr3, aioA, arsB, arsC (grx), arsC (trx), arsD, arsM, arrA, and arxA. FunGene databases [27] were constructed for the following arsenic-related genes: arsB, arsC (grx), arsC (trx), acr3, aioA, arrA, and arxA. The arxA database was constructed with seed sequences from [12]. For all other genes, UniProt [65] was used to obtain full-length, reviewed sequences when possible. NCBI clusters of orthologous groups (COG) [66] for each gene were examined for evidence of function in the literature. All COG and UniProt sequences were aligned using MUSCLE [67]. Aligned sequences were included in a maximum likelihood tree with 50 bootstrap replications made with MEGA (v7.0, [68]). Sequences that did not cluster with known sequences and had no evidence of function were removed. A final FASTA file for each gene was submitted to the Ribosomal Database Project (RDP) to construct a FunGene database [27]. All arsenic-related gene databases are freely available on FunGene (http:// fungene.cme.msu.edu/).

Arsenic-related genes in cultivable soil microorganisms
The RefSoil+ database [37] was used to obtain high-quality genomes (chromosomes and plasmids) from soil microorganisms in the Genomes OnLine (GOLD) database [69]. RefSoil+ chromosomes and plasmids were searched with hmmsearch [29] using HMMs from Fun-Gene with an e-value cutoff of 10 − 10 . The top hits were analyzed in R [70]. For each gene, scores and percent alignments were plotted to determine quality cutoffs. Stringent percent alignment scores were included since this search was against completed genome sequences: only hits with scores > 100 and percent alignment > 90% were included. Hits with the lowest scores were manually examined to test for false positives. Due to false positives, hits against aioA, arrA, and arxA were further quality filtered to have scores > 1000. When one open reading frame (ORF) contained multiple hits, the hit with a lower score was removed. Taxonomy was assigned using the RefSoil database [42], and the relative abundance of arsenic-related genes within phyla was examined. A 16S rRNA gene maximum likelihood tree of RefSoil+ bacterial strains was constructed with RAxML (v.8.0.6 [71]) based on the Whelan and Goldman (WAG) model with 100 bootstrap replicates ("-m PROTGAM-MAWAG -p 12345 -f a -k -x 12345 -# 100"). Based on accession numbers, gene hits were extracted from RefSoil+ sequences and used to construct maximum likelihood trees for each gene.

Reference database construction
Reference gene databases of diverse, near full-length sequences were constructed using limited sequences from FunGene databases [27] for the following genes: acr3, aioA, arrA, arsB, arsC (grx), arsC (trx), arsD, arsM, and arxA. Seed sequences and hidden Markov models (HMMs) for each gene were downloaded from FunGene, and diverse protein and corresponding nucleotide sequences were selected with gene-specific search parameters (Additional file 11). Briefly, minimum amino acid length was set to 70% of the HMM length; minimum HMM coverage was set to 80% as is recommended by Xander software for targeted gene assembly; and a score cutoff was manually selected based on a dropoff point. Sequences were de-replicated before being used in subsequent analysis, and final sequence counts are included in Additional file 11. Reference databases were converted to publicly available BLAST databases using BLAST+ [30]. Reference and BLAST databases are publicly available on GitHub (https://github.com/ShadeLab/PAPER_ Dunivin_meta_arsenic)

Sample collection and preparation
A soil surface core (20 cm depth and 5.1 cm diameter) was collected in October 2014 from Centralia, PA (GPS coordinates: 40 48.070, 076 20.574). For cultivation-dependent work, a soil slurry was made by vortexing 5 g soil with 25 mL phosphate-buffered saline (PBS) for 1 min. Remaining soil was stored at − 80°C until DNA extractions. The soil slurry was allowed to settle for 2 min. One hundred microliters of the slurry was then removed and serial diluted using PBS to a 10 − 2 dilution. One hundred microliters of the solution was added to 50% trypticase soy agar (TSA50) with 200 μg/ mL cycloheximide to prevent fungal growth. Plates were incubated at 60°C for 72 h. Lawns of growth were extracted by adding 600 μL trypticase soy broth with 25% glycerol to plates. The plate scrapings were stored at − 80°C until DNA extraction.

DNA extraction and metagenome sequencing
DNA for cultivation-independent analysis was manually extracted from soil using a phenol chloroform extraction [72]

Public soil metagenome acquisition
In total, 38 soil metagenomes were obtained for this work (Additional file 1). Datasets from Centralia, PA, were generated in our research group. All other metagenome datasets were obtained from MG-RAST (http:// metagenomics.anl.gov/). The MG-RAST database was searched on May 15, 2017, with the following criteria: material = soil, sequence type = shotgun, public = true. The resulting list of metagenome datasets was ordered by the number of base pairs (bp). Metagenomic datasets with the most bp were only included if they were sequenced using Illumina to standardize sequencing errors, had an available FASTQ file for internal quality control, and contained < 30% low quality as determined by MG-RAST. Within high-quality Illumina samples, priority for inclusion was given to projects with multiple samples so that comparisons could be made both within and between soil sites. When a project had multiple samples, datasets with the greatest bp were selected. While we prioritized samples with multiple datasets, several replicate samples were omitted early on due to > 30% of data removed during quality filtering, and samples Illinois soil, Russian permafrost, and Wyoming soil have just one sample. This search ultimately yielded 26 datasets from 12 locations and 5 countries (Additional file 2).

Soil metagenome processing and gene targeted assembly
Sequences from MG-RAST datasets as well as Centralia sample Cen13 were quality controlled using the FASTX toolkit (fastq_quality_filter, "-Q33 -q 30 -p 50"). Twelve datasets from Centralia, PA, were obtained from the Joint Genome Institute and quality filtered as described previously [73]. Quality-filtered sequences were used in all downstream analyses. For each dataset, a gene targeted metagenome assembler [28] was used to assemble each gene of interest. For each gene of interest, seed sequences, HMMs, and reference gene databases described above were included. For rplB, reference gene database, seed sequences, and HMMs from the Xander package were used. In most instances, default assembly parameters were used except to incorporate differences in protein length (i.e., protein is shorter than default 150 amino acids) or to improve quality (i.e., maximum length is increased to improve specificity) (Additional file 11). While the assembler includes chimera removal, additional quality control steps were added. Final assembled sequences (operational taxonomic units, OTUs) were searched against the reference gene database as well as the non-redundant database (nr) from NCBI (August 28, 2017) using BLAST The funders had no role in the design of the study and collection, analysis, and interpretation of data.