Salt flat microbial diversity and dynamics across salinity gradient

Sabkha soil analysis

Soil biochemical analysis shows higher concentration of salt in the middle (average of 200 ppt) which decreases to 20 ppt outside of the sabkha (Fig. 1B). Soil conductivity and total dissolved solids (TDS) were positively correlated with salinity (r = 0.8, p < 0.05). Cl⁻; Na₂O was only abundant in the middle, whereas Cu, Fe₂O₃, P₂O₅, SiO₂ and Al₂O₃ were abundant in the inner and edge of the sabkhas compared to the outside. In contrast, SO₃ was more pronounced on the outside (Fig. 1C, Supplementary Tables 1, 2).

Dynamic shift in metagenomic profile across sabkha salinity gradient

Species diversity within the samples (alpha-diversity) were measured by Shannon diversity index method. The highest diversity was identified in the high salt regions Site I (diversity: 4–4.8) and the edge of sabkha (E) samples (diversity: 4.1–4.6) (Fig. 1D). The lowest microbial diversity (2.5–3.8) was observed in the samples collected from the combined sand and overflow sabkha region (outside). Alpha diversity measures for upper, lower and middle across the eighteen samples show no significant difference among layers (ANOVA F-value = 0.86572, p = 0.53) (Supplementary Fig. 1). The difference in the diversity is only observed between inside and outside the sabkha, which correlates with the salt gradient (ANOVA F-value = 18.962, p-value = 7.8 × 10^–5). Diversity between the samples (Beta diversity) analysis of the sabkha microbial composition along the three different sites (I, E, O) identified similar microbial populations between inside and edge sabkha samples, clustered together in a PCoA analysis. Samples collected from outside sabkha region formed a separate cluster in the PCoA analysis (Fig. 1E,F) (PERMANOVA F-value = 45.062, R² = 0.857, p-value < 0.001).

Taxonomic identity shows that there is shift in abundance between inside and outside the sabkha (Fig. 2A,B) and it is correlated with the salt gradient (R = 0.85, p < 0.05). At phylum level, Bacteroidetes (~ 22 to 34%), Euryarchaeota (~ 18 to ~ 30%), unclassified bacteria (~ 24 to ~ 35%), Actinobacteria (~ 0.01 to ~ 11%) and Cyanobacteria (less than 1%) are predominantly found in both inside and edge sabkha regions, whereas Proteobacteria (~ 92 to ~ 97%) and Parcubacteria (~ 1 to ~ 2%) are predominant in outer sabkha, Cyanobacterial classes were found to be over-represented inside the sabkha compared to the outside. Several archaeal taxa, as well as some unclassified bacteria, were rich in the inner and edge sabkha samples, while halophilic bacteria were overrepresented on the inside hypersaline region of the sabkha (Fig. 2A). Halobacteriaceae, Actinobacteria, Planctomycetaceae, Bacteroidetes and Planctomycetes are found more abundant in inside sabkha. Gammaproteobacteria, Rhodobacteraceae, Alphaproteobacteria, Rhizobiales and Oceanospirillales are found more abundant in the hyposaline outer sabkha. We see a shift in the Proteobacteria structure (Fig. 2B, Supplementary Fig. 2), where Delta/Epsilon Proteobacteria were more abundant in the hypersaline region inside the sabkha, while Gammaproteobacteria, and Alphaproteobacteria were more abundant in the hyposaline outer sabkha (Fig. 2B). Interestingly, we see similar taxonomic groups in the middle and lower layer of inside sabkha and less groups in the upper salt crust of each section (Fig. 2B), but not on the outside the sabkha.

For the sabkha whole metagenome assembly, the 4 sites sampled on inside, edge (s1, s3, s8, s12) and the 4 sites outside (s13, s14, s17, s18) (Fig. 1B) were used for assembly to resolve the extent of genome completeness. In total ~ 830 million PE reads (150 bp) were generated (Supplementary Table 3) for the metagenome assembly, resulting in 664,455 contigs (final genome assembly length ~ 2.3 Gb) with the N50 value of 3.945 Kb (Supplementary Table 4). All the assembled contigs were taxonomically classified (at phylum level) based on the sample location. The contigs related to the phyla Euryarchaeota, Bacteroidetes, Actinobacteria, Planctomycetes, Cyanobacteria, and Firmicutes are predominantly found in the inside sabkha samples, whereas phyla Proteobacteria is predominant in outside sabkha region (Fig. 3A,B). In addition, contigs were binned into 225 genomes across the sabkha (Fig. 4A). Assembled genome bacterial size ranges from ~ 0.5 to ~ 9.3 Mb, with the N50 of the assembled genome from ~ 2 Kb to ~ 1 Mb (Supplementary Table 5). From the assembled bins, we annotated CDS, tRNA, rRNA and tmRNA sequences (Supplementary Table 5). In total, 651,373 protein coding (CDS) and 6698 noncoding genes were annotated from the assembled genome. The phylogenetic tree generated using assembled genomes is shown in (Fig. 4B).

Taxonomically resolved genomes show high nucleotide diversity (π) levels on the inside compared to the outside of the sabkha (Fig. 4C). We observed comparable taxonomic microbial profiles using whole metagenome and 16S rDNA analysis. Interestingly, Halobacteria have the highest diversity on the inside followed by Rhodotermales, whereas Rhizobiales and Pseudomonadales are highest from the outside (Fig. 4C).

The pairwise genome-wide ANI (Average Nucleotide Identity) and alignment coverage of 125 dereplicated but quality-filtered genomes across the sabkha showed a decrease of pairwise ANI at ~ 97% (Supplementary Fig. 3), which is comparable to other reported thresholds in bacterial species⁵². The clustering of genomes using dRep into groups of species-like populations using ANI = 97% (Supplementary Fig. 4) was done and genomes with > 74% completeness and less than 10% contamination kept for further analysis. The focus on abundant populations/sites with more than 1% abundance resulted in 129 dereplicated genomes with a least 3 genomes per population summarized in Fig. 4D. Using metagenome-assembled genomes, the measure of Irep index show that the inside and outside is fully replicating (Fig. 4D). There is statistically significant difference (t test p = 2.43 × 10⁻⁵) between Irep for overall taxonomic order group from hypersaline inside (1.744 ± 0.274) compared to hyposaline outside (1.527 ± 0.145) sabkha (Supplementary Table 6). On the inside, we observe Halobacteriale with an average of 1.754 ± 0.184, while on the outside, Pseudomonadales has 1.554 ± 0.124, which is significantly lower (t test p = 0.02) than inside sabkha. Rhizobiales on the outside is replicating more than the overall outside with Irep around 1.632 ± 0.218 (Fig. 4E). We found no correlation between Irep and Pi for the different taxonomic order (Pearson R = − 0.12, p = 0.53) (Supplementary Fig. 5, Supplementary Table ,6), which suggest involvement of replication of these bacteria in nutrient cycling of sabkha (Supplementary Fig. 7).

Nucleotide diversity in the sabkha

The overall combined microbial population nucleotide diversity (π) was highly significant for all inner microbial communities compared to outside the sabkha (Fig. 5A) as is evident for each sampling site on the inner (s1, s3, s8, s12) compared to outer (s13, s14, s17, s18) (t test Bonferroni p < 0.05) (Supplementary Fig. 8). Clearly, Halobacteriales and Pseudomonadales populations have high diversity from inside the sabkha, while Rhizobiales and Pseudomonadales have highest diversity from outside the sabkha (Supplementary Fig. 8). The comparable diversity of combined as well as per site samples suggests that diversity persists within soils at even finer scale. We also found comparable ranges of per gene nucleotide diversity values (Pearson R = 0.9, p < 0.05) (Supplementary Fig. 9) for different populations at 96% and 98% ANI (Supplementary Fig. 10).

We found in almost all sabkha populations that ribosomal genes have lower nucleotide diversities that the average for genes in the genome (Fig. 5A), consistent with purifying selection⁵³. An exception is for the inner sabkha Halobacteriales population, and outer sabkha populations such as Propionibacteriales, Pseudomonadales, Rhizobiales, and Rhodobacteriales, which have higher ribosomal nucleotide diversities, suggesting that they play an adaptive role in diverse ecological niches and microenvironment. Moreover, biosynthetic genes annotated with antiSMASH have significantly higher diversity than the genome average (pairwise t test Bonferroni p < 0.05), except for Halobacteriales, Thermoplasmatota, Pseudomonadales, and Rhizobiales.

When we looked at the diversity of different classes of biosynthetic genes from inside and outside the sabkha, we find unique patterns of biosynthetic gene class counts (Supplementary Fig. 11, Fig. 5B). Overall shared classes of inside and outside the sabkha include terpene, bacteriocin, betalactone, hserlactone, T3PKS, T1PKS, siderophores and ectoine (Supplementary Fig. 11, Fig. 5B). Inside the sabkha, however, we observe for Halobacteriales LAP (linear azol (in) e-containing peptides) that are synthetized in ribosomes, lanthipeptide and butyrolactones. On the other hand, we observed in outside of the sabkha, we see biosynthetic genes from different unique populations classes, such arylpolyene, CDPS, furan, lassopeptide, nucleoside, NRPS (NonRibosomal Peptide-Synthetase), phenazine, phosphonate, resorcinol and TfuA-related (Supplementary Fig. 11, Fig. 5B).

Horizontal gene transfer, homologous recombination and sabkha microbial diversity

Our analysis at the order level show that there is extensive horizontal gene transfer on the inside and outside the sabkha (Fig. 6A). On the inside, there is an exchange of genetic information between different levels especially Halobacteriales–Rhodotermales, reaching close to 42 HGT events (Fig. 6A, Supplementary Fig. 12). On the other hand, our analysis shows 30–47 HGT events between Rhizobiales–Rhodobacteriales, Kiloniellales–Rhizobiales and Pseudomonadales–Nevskiales outside the sabkha (Fig. 6A, Supplementary Fig. 12).

Using PFAM domains, we observed enriched transfer of abiotic stress genes, detoxification, biosynthetic of small molecules, transport, homeostasis and transport genes associated with HGT events (hypergeometric test, p < 0.05) (Fig. 6B). Two examples of HGT are highlighted in Fig. 6C. Interestingly, we see a gene for social cooperative development inside the sabkha. Cooperation among bacteria is a common strategy of interaction via cell–cell signaling, termed ‘quorum sensing’ (QS), which is an evolutionary stable strategy under extreme abiotic or biotic stress⁵⁴.

Measuring linkage disequilibrium (LD) of SNPs spanned by paired end reads^48,55 in different bacterial population can tell us on the impact of homologous recombination and its effect on genetic diversity. Our results are consistent with the expectation that natural microbial population experience a large amount of recombination, and this is observed with linkage disequilibrium (r²) that decay as the distance increases between two SNPs, both inside and outside the sabkha (Fig. 7A, Supplementary Fig. 13). The estimate of the ratio of the neutral rate of recombination to mutation in synonymous third position (Fig. 7B) shows a range for different species on the inside and outside of the sabkha. High rates of recombination is observed in some Halobacteriales, Balneolales, Enterobacteriales, Pseudomonadales, Rhizobiales and Rhodobacteriales species. On other hand, we observe species like Bacillales on the inside and some Pseudomonadales on the outside with very low rates of recombination (Fig. 7B). We also note, the ratio of recombination of nonsynonymous to synonymous (r2N/r2S) increases with higher diversity population (Fig. 7C; Supplementary Fig. 14).

We also measured the LD measure, namely, D′, which is less than 1 suggesting all possible combinations of a pair of biallelic sites, an serves as indication that recombination is present. Our results for inside the sabkha shows that D′ is correlated with mean r² (R² = 0.78, p = 0.005) (Fig. 7D); this correlation is similarly found outside the sabkha (Supplementary Fig. 15). D′ was less than 1 in 9–21% of all site pairs in inside sabkha populations, while D′ < 1 in 2–12% outside the sabkha (Supplementary Fig. 16; Supplementary Table 7). Given that D′ less than one is observed in a significant fraction of loci, as well as the decay of linkage with genomic distance suggests that homologous recombination in the sabkha microbiome.

Population dynamics and gene-specific selective sweeps across salinity gradient of sabkha

Changes in factors such as soil geochemistry, particle structure, salinity, and many other cues make soil ecosystems heterogeneous, as these factors can change over millimeter scales⁵⁶. Given the complexity in examining the effect of changing abiotic environmental factors over small spatial scales, we estimated the pairwise fixation index Fst for each gene between inside sites (s1, s3, s8, s12) and outside sites (s13, s14, s17, s18); these display a gradual decrease in salinity between sabkha sites (Fig. 8A). Our results show that Fst per populations between sites increases as salinity increase among sites (Fig. 8A), consistent with geographic population structure at most loci across the genome (Supplementary Fig. 17). In contrast, a fraction of Fst values were low (< 5%), consistent with dispersal of alleles among sites.

Population-specific selective pressures acting on specific loci have a signature of high Fst compared to the average genome background⁵⁷. We identified loci with high Fst using a sliding window of 5 genes and tested if these regions had Fst values greater than 1.5 standard deviations above the mean genome average (Supplementary Table 8). The length of genomic region with high Fst was defined when the mean Fst dropped below this threshold. We identified 72 loci with high Fst, in different microbial genomes (e.g., Halobacteriales 166 s1 vs. s3) despite having low average Fst values (Fig. 8B, Supplementary Fig. 18).

However, the low nonsynonymous to synonymous ratios in these genes that suggests a signal of purifying selection (Supplementary Fig. 19), as well as the reduction in sequencing coverage at some genomic regions, are both correlated with nucleotide diversity (Supplementary Fig. 20). This suggests that we cannot rely solely on identifying selective sweeps by a significant reduction in nucleotide diversity (Fig. 8C), therefore we looked for a significant increase in LD at these loci compared to the genomic background. Using this criterion, we identified 18 loci with an observed significant reduction in nucleotide diversity and increased LD for different population (Halobacteriales, Pseudomonadales, Nitriliruptorales, Bipolaricaulia) within inner and outside sabkha sites (Fig. 8C, Supplementary Fig. 21) suggesting selective sweep events at these genes at different sites in sabkhas. Several of these genes are related to uptake, movement and transport while others have functions related to abiotic stress sensing and response regulation (Supplementary Table 9). Selective sweeps are also evident between sabkha edge and outside sabkha sites (Fig. 8B) (e.g., Pseudomonadales 194), which together suggests that selection is driving some of the differences in microbial population structure across this ecosystem.

Diverse metabolic potential of microbial communities of sabkha

We investigated the ability of microbial communities of sabkha in harboring specific metabolic activity, and nutrient cycling (e.g.; carbon cycle, Nitrogen cycle, sulfur cycle) (Supplementary Fig. 22). We found, Perchlorate reduction and Chlorite reduction which are enriched on the inside of sabkha (Supplementary Table 10). Interestingly, when we look at the number of genes for combined metabolic pathway; on the inside extreme Halophilic microbes (e.g., Halobacteriales, Cyanobacteriales) have more than 50 percent of genes (t test p < 0.05), while outside, moderately halophilic microbes (e.g., Pseudomonadales, Rhizobiales, Nitriliruptorales) have them (t test p < 0.05) (Supplementary Fig. 6, Supplementary Tables 10, 11).

We found evidence for multiple carbon fixation pathways on the inside sabkha, including the Calvin–Benson–Bassham (CBB) cycle in Cyanobacteriales, Halobacteriales, Kiloniellales and Thiohalorhabdales, the 3-hydroxypropionate bicycle (3HP bicycle) in Enterobacterales and Pseudomonadales and the reductive citric acid cycle (rTCA) only found in Nanosalinales. In addition to the latter, carbon fixation pathways that are found in most of the outside metagenomes, the outside of sabkha harbors 3HP/4HB cycle, which is only found in Acidimicrobiales, Burkholderiales, Nitriliruptorales and Polyangiales (Supplementary Tables 10, 11, Supplementary Fig. 7). The key enzyme methyl-coenzyme M reductase (mcr) for methanogenesis was absent from sabkha (Supplementary Tables 10, 11, Supplementary Fig. 7). However, methanotrophy was present through soluble methane monooxygenase (mmoB) gene in Halobacteriales and Pseudomonadales from inside and in Pseudomonadales, Mycobacteriales, and Solirubrobacterales from outside. We see three abundant genes encoding aerobic monoxide dehydrogenase (coxSML) on the inside of sabkha and found only in our Halobacteriales, Kiloniellales and Bipolaricaulia, whereas it is only in Rhizobiales, Rhodobacteriales and Nitriliruptorales from outside (Supplementary Tables 10, 11).

In term of Nitrogen cycling, we did not observe any metagenomic evidence for nitrogen fixation for inside or outside sabkha. On the other hand, we found evidence for the first step of nitrification and denitrification process, using nitrite oxidoreductase (nxrAB), which is involved in nitrite to nitrate and nitrate reductase (narGH), for the reverse reaction (Supplementary Tables 10, 11; Supplementary Fig. 7). The second step of denitrification involves 2 genes (nirSK, norBC) and nirK is much abundant on the inside in Halobacteriales, Kiloniellales and Bipolaricaulia than nirS. nirK genes are found only in 6 metagenomes from the outside, such as in Rhodobacteriales, Pseudomonadales. As for the last step of dissimilatory nitrate reduction to ammonia mediated by nitrite reductase, we see abundance of genes encoding the large and small subunit of the NADH-dependent for of the enzyme (nirB, nirD) on the outside of sabkha and less abundant on the inside and this shift is the same genes encoding the cytochrome c for of the enzyme (nrfAH) (Supplementary Tables 10, 11, Supplementary Fig. 7).

In terms of Sulfur cycling, genes encoding sulfide-quinone reductase (sqr) are not abundant in Sabkha, with 3 genomes from the inside (Cyanobacteriales, Enterobacteriales and Pseudomonadales) and 3 genomes on the outside (Nitrosococcales, Gammaproteobacteria, Enterobacteriales). Instead, others genes such as Sulfur dioxygenase (sdo) and Sulfate adenylyltransferase (sat) are equally abundant on both sides of sabkha. Genes associated with thiosulfate oxidation (soxBY) were most abundant in the outer sabkha and included in metagenomes such as Burkholderiales, Rhizobiales, Rhodobacteriales and Rhodospirillales. The inside has only soxY and soxB genes in Thiohalorhabdales and Kiloniellales respectively (Supplementary Tables 10, 11; Supplementary Fig. 7).

We see evidence for aromatic compounds, complex carbon and fatty acid biodegradation pathways for bacteria and archaea metagenomes. Other pathways, such Arsenic cycling, Urea utilization, oxidative phosphorylation, fermentation and Amino acid utilization are metabolic functions that sabkha ecosystem is harboring (Supplementary Tables 10, 11).

Source: Ecology - nature.com