Comparative metagenomic analysis of soil microbial communities across three hexachlorocyclohexane contamination levels

Abstract

This paper presents the characterization of the microbial community responsible for the in-situ bioremediation of hexachlorocyclohexane (HCH). Microbial community structure and function was analyzed using 16S rRNA amplicon and shotgun metagenomic sequencing methods for three sets of soil samples. The three samples were collected from a HCH-dumpsite (450 mg HCH/g soil) and comprised of a HCH/soil ratio of 0.45, 0.0007, and 0.00003, respectively. Certain bacterial; (Chromohalobacter, Marinimicrobium, Idiomarina, Salinosphaera, Halomonas, Sphingopyxis, Novosphingobium, Sphingomonas and Pseudomonas), archaeal; (Halobacterium, Haloarcula and Halorhabdus) and fungal (Fusarium) genera were found to be more abundant in the soil sample from the HCH-dumpsite. Consistent with the phylogenetic shift, the dumpsite also exhibited a relatively higher abundance of genes coding for chemotaxis/motility, chloroaromatic and HCH degradation (lin genes). Reassembly of a draft pangenome of Chromohalobacter salaxigenes sp. (∼8X coverage) and 3 plasmids (pISP3, pISP4 and pLB1; 13X coverage) containing lin genes/clusters also provides an evidence for the horizontal transfer of HCH catabolism genes.

Figure 1. Phylogenetic analysis of the microbiomes.

(A) Dual dendrogram of top 50 bacterial genera across three metagenomes obtained after TEFAP analysis using four bacterial primer sets. Genera and sample categories were clustered using Manhattan distance metric, top 50 genera with standard deviation >0.4 and having at least 0.8% of the total abundance were selected. Colour scale is representing the relative abundance of sequence reads (normalized by sample-mean). (B) Phylogenetic correlation of microbial communities across increasing HCH contamination, a subset of 1000 randomly selected OTUs from each metagenome was used to construct an elucidan distance matrix. Matrices were pair-wise compared using Mantel-test (1000 permutation, 0.05 as standard P -value) and Pearson correlation values were calculated. Asterisks indicate the statistical significance P<0.001(mean±sm). (C) Relative percentage of reads assigned to different archeal (I) and fungal (II) genera in TEFAP analysis.

Figure 2. Functional traits of the studied metagenomes.

(A) Cellular processes enriched over increasing HCH contamination. Metagenomic reads were compared against the COG database and relative percentage (y-axis) for each category (x-axis) was calculated. (B) Heat map showing the relative abundance of top 50 subsystems enriched over increasing HCH concentrations (percentage cut-off = 0.8%, standard deviation cut-off = 0.4%). (C) Rarefaction analysis performed on unique protein families (Pfam) sampled across three HCH gradients. (D) Comparison of functional categories similarity between metagenome gradient pairs. KEGG enzyme profile of each metagenome was compared. Asterisks indicate significant differences (Two sided Fishers exact test with Bonferroni multiple test correction, P<0.01). ABBREVATIONS: (1) DS = dumpsite gradient, 1 km = 1 km gradient and 5 km = 5 km gradient.

Figure 3. Enrichment of lindane degradation (aerobic) pathway.

Schematic representation for the enrichment of aerobic degradation pathway of lindane. Numerical values (on color gradient) at each enzyme represent the diversity (genera) of the corresponding gene present at each metagenome estimated using Transpipe analysis.

Figure 4. DNA-seq analysis of the community potential for HCH degradation.

DNA-seq analysis of metagenome sequences against reference lin genes using Array star, x axis represents the relative abundance of lin genes from different genera (Table S8) present at studied metagenomes.

Figure 5. Quantifying the enrichment of environmental genomes/plasmids.

Metagenomic recruitment plots of genomes/plasmids constructed using all three studied metagenome sequences. Reads were mapped with coverage parameter. (A) Assembled contigs of Chromohalobacter salexigens (5189 contigs) from mtegenomic reads, shaded region represents the location of 16SrRNA gene sequence. (B) pISP4, (C) Sphingobium japonicum UT26 chromosome 1, (D) pISP3 and (E) pLB1. Localization of lin genes on respective genomes is marked along with representation symbols for IS-elements.