Sequencing yield and pre-processing
To investigate the microbiome associated with native plants that grow in phosphate mines in Morocco and its comparison with the phosphate bulk soil samples and the rhizosphere of wheat crop, a total of 15 samples were subjected to 16S rRNA metagenomics analysis using Ion Torrent PGMTM sequencing platform. Total obtained reads accounted for 2 782 897. they were filtered for chimeras, quality score and copy number. Reads were classified according to 16S rRNA variable regions for each sample and only reads with high number of copies were mapped against 16S rRNA databases (Supplementary materials Table S1 and Table S2).
Taxonomic resolution of the primers targeting variable regions of the 16S rRNA
Variable regions of the 16S rRNA were amplified in the DNA derived from the studied samples and the ZymoBIOMICS community standard using two separate primer sets with the Ion 16STM metagenomics kit. Read Classification to taxonomic profiles at four levels e.g., phylum, family, genus and species was compared between the consensus and each region according to the number of taxa identified. In this section, results obtained for the standard ZymoBIOMICS DNA were also presented.
Figure 1 showed that at the phylum level, V2, V3, V4, V67 and V8 performed well for the identification of the phyla present in the ZymoBIOMICS standard. However, V9 region identified only Proteobacteria phylum. Taxonomic resolution at more lower levels became challenging and V9 region gave less resolution up to family, genus and species levels. V67 and V8 regions gave better resolution than V9, nevertheless, new families and genera were found with these regions that are not present in the standard. V2, V3 and V4 regions gave the best resolution up to both family and genus levels (Figure 1).
The same analysis was performed for the collected soil samples ( Figure 2). However, since the bacterial composition of these samples is unknown compared to the standard ZymoBIOMICS, the maximum number of each taxonomic level present in each sample was taken as the target resolution. For each region, resolution percentages were calculated accordingly and presented as heatmap. The overall analysis of the figure 2 indicated variability in the taxonomic resolution of different 16S rRNA regions according to the taxonomic level and the sample types. However, some regions gave the lesser taxonomic resolution e.g., V9 and V2 whereas V3, V4, V67 and V8 performed better in the taxonomic resolution mainly at the phylum level. At higher taxonomic levels e.g., Genus and Species, the number of taxa recovered by every hypervariable region decreased for almost all the regions except the V3 region.
Bacterial diversity analysis by consensus sequences in the studied samples:
16S Metagenomics analysis included samples from bulk phosphate mine, samples from the rhizosphere of three native wild plants in phosphate mine and samples from the wheat rhizosphere. Figure 3 presented the distribution and abundance of bacterial phyla, families and genera in the studied samples using the consensus sequences from the all variable regions of the 16s rRNA gene. Results of the bacterial composition clearly showed differences according to sample types and the plant species. Bacterial composition of the samples obtained from bulk phosphate mine in the vicinity of the wild plants was dominated by Actinobacteria with average relative abundance of 97,73% (figure 3a). Proteobacteria was Less abundant in these samples compared to Actinobacteria. At the family and genus levels, an average of 94,54% and 94,29% of OTUs were assigned to Nocardioidaceae and Nocardioides, respectively in these samples. This result indicated the selective presence and dominance of Nocardioides members in this particular environment among the bacterial diversity. As showed in figure 3, the bacterial composition was completely different in the rhizosphere samples of the three native plants. Proteobacteria phylum was almost enriched than other phyla with relative abundance of 62,24%, 71,15% and 65,61%, in the rhizosphere Alternanthera caracasana (NP1); Digitaria sanguinalis (NP2) and Dittrichia Viscosa (NP3), respectively. Within this phylum, main representative classes were Betaproteobacteria (7,1%; 40,38%; 23,15% in NP1, NP2 and NP3, respectively), Alphaproteobacteria (26,84%; 19,65% ; 22,85% in NP1, NP2 and NP3, respectively) and Gammaproteobacteria (28,05%; 7,85%; 15,68% in NP1, NP2 and NP3, respectively). Deltaproteobacteria were less abundant. High relative abundances of Actinobacteria (22,53%, 15,24%, 22,30% in NP1, NP2 and NP3, respectively) were also found in these native plant samples, however, at lower level, Actinobacteria were mainly represented by Actinomycetales order (21,81% ; 12,58; 19,17% in NP1, NP2 and NP3, respectively) with wide taxonomic diversity at family and genus level. It is noteworthy that the relative abundance of Nocardioidaceae family was significantly reduced in these samples in comparison with the bulk phosphate mine samples (Figure 3). Sequences assigned to Actinobacteria were further analyzed in this study. Other major phyla that were found in these samples were Bacteroidetes (7,57% ; 4,23% ; 7,63% in NP1, NP2 and NP3, respectively), and Firmicutes (5,82%; 1,17% ; 2,83% in NP1, NP2 and NP3, respectively). Broad taxonomic diversity of minor phyla was also found in the native plants samples and included Acidobacteria, Armatimonadetes, Chloroflexi, Cyanobacteria, Gemmatimonadetes and Verrucomicrobia (Figure 3). The wheat rhizosphere samples showed a relatively similar pattern as the bacterial composition of the rhizosphere of native plants. In these samples, the assignment of taxonomic affiliation revealed high relative abundances of Proteobacteria (57,42%), Actinobacteria (35,24%), Firmicutes (5,52%) and Nitrospirae (2,32%). Broad taxonomic diversity of minor phyla was also found in the wheat rhizosphere samples. Actinobacteria were mainly represented by Actinomycetales and Gaiellales orders with 15,53% and 9,36%, respectively.
Alpha diversity measures (Observed OTUs, Chao, ACE and Shannon) and rarefaction curves are presented in Figure 4 and supplementary material Figure S1. Rarefaction analysis enabled us to gauge the extent to which total diversity had been recovered at each sample and indicated that the extent of bacterial diversity was captured by the sequencing approach. Figure 4 showed significantly lower diversity and OTUs number in the bulk mine samples. This finding was expected, because of the extreme conditions met in such habitat. Interestingly, the native plants in the phosphate mine showed higher observed OTUs, species richness and species diversity and evenness in comparison with the wheat crop rhizosphere samples, even if the observed differences were statistically not significant (Figure 4). Higher alpha diversity measures were recorded for the rhizosphere samples of the phosphate mine native plant Dittrichia Viscosa (NP3), this increase was statistically significant with other native plants NP1 and NP2 samples (Figure 4). In contrast to other samples, the wheat rhizosphere samples were characterized by variable intra-diversity, which makes the box-plot range larger than in other samples. Some replicates were highly diverse and rich in terms of species, while others were characterized by their lower richness in bacterial communities (Figure 4).
Beta diversity measures included Bray Curtis distance unweighted and weighted UniFrac distances and jaccard dissimilatory. PCoA of the latter was presented in figure 5 and clearly indicated differences in bacterial diversity among at least three kind of samples : The bulk phosphate mine, the rhizosphere of the native wild plants of the phosphate mine and the rhizosphere of the wheat crop. Statistical analyses using PERMANOVA revealed significant differences between the three groups of samples (p=0,021). This finding indicated that the plant genotypes had limited influence in the shaping the bacterial diversity in this particular environment.
The distribution of shared and unique OTUs among bacterial communities in the studied samples was presented in figure 6. High number of OTUs was shared between the samples at the phylum level such as Actinobacteria, Proteobacteria,Firmicutes and Bacteroidetes (Figure 6a). The shared OTUs between the rhizosphere samples (native plants and wheat crops) included Armatimonadetes, Planctomycetes and Verrucomicrobia, whereas only Nitrospirae phylum was shared between the bulk mine samples and the rhizosphere of native plants. Similarly, at lower taxonomic levels e.g., family (Figure 6b) and genus (Figure 6c), high number of OTUs was shared between the samples (28 OTUs among 86 and 33 among 317 at the family and genus levels, respectively). Interestingly, high number of OTUs at the family level (27) were shared between the rhizosphere samples of native plants and wheat crops which included some plant-associated bacteria such as Rhizobiaceae and Phyllobacteriaceae, and indicated the effect of the plant rhizosphere in the interactions with bacterial diversity. The same statement is true for the shared OTUs at the genus level with high number among the plant samples (figure 6c)
Actinobacteria phylogenetic analysis :
According to taxonomy and representative sequences tables of Qiime2, the number of sequences assigned to Actinobacteria varied between samples. Bulk phosphate mine samples accounted for 168 sequences assigned to Actinobacteria, whereas 336 sequences were assigned to Actinobacteria in the rhizosphere samples of the wild plants in phosphate mine. Figure 7 showed the phylogenetic relatedness and the inferred evolutionary relationships between the Actinobacteria sequences in the studied samples and revealed that Actinobacteria in Bulk RP were phylogenetically different and separated from those in native plant samples. Actinobacteria related sequences in the bulk phosphate mine samples yielded three different clades whereas, those in native plants samples were more diverse with many branches. Actinobacteria related Sequences in NP1 and NP3 were the most closely related among all samples with a common node of 72% (Figure 7). Actinobacteria in the Bulk phosphate samples were less evolutive than the sequences in the native plants samples which could be maybe related to adaptation and co-evolution process between the bacteria and the plants All together these results indicated the selection of the Actinobacteria by plants in their rhizosphere and the absence of possible migration of some Actinobacteria from the bulk mine samples to the plants rhizosphere in their surrounding environment.
Prediction of Bacterial functions enriched in the samples :
The Bioinformatic tool Tax4fun was used in this study to predict bacterial functions involved in nitrogen and phosphorus metabolism and stress across the studied samples (Figure 8). Figure 8a showed the enrichment of some nitrogen fixation functions, mainly in the bacteria of the rhizosphere of the phosphate mine native plants Alternanthera caracasana (NP1) and Dittrichia Viscosa (NP3). The enriched bacterial functions in these samples included nitrogen fixation proteins (mainly NifQ, NifB, NifN, NifT, NifX, NifZ, NifE) and nitrogenase cofactors. The nitrogen fixation protein NifU was mainly enriched in bacteria assessed in wheat rhizosphere and the phosphate mine native plant Digitaria sanguinalis (NP2). Nitrogen regulatory proteins PII 1 and PII 2 KEGG were also enriched in bulk phosphate samples and some native plants and are linked to the regulation of carbon metabolism under conditions of low nitrogen availability. Heatmap of the bacterial functions involved in phosphate metabolism (Figure 8b) showed the enrichment of 11 predicted KEGG functions involved in mineralization, solubilization and transport (Phosphatase, acid phosphatases and inorganic phosphate transporter) among bacterial communities in the bulk phosphate mine samples. Interestingly, the assigned bacteria to the rhizosphere of native wild plants showed similar patterns, with the enrichment of many functions involved in the mineralization, solubilization and transport of phosphate. In the wheat rhizosphere samples, phosphate transport functions were mainly enriched with some phosphatase involved the phosphate mineralization processes. Figure 8C showed the enrichment profiles of bacterial functions related to stresses, it indicated that the Na+/H+ antiporters NhaC and NhaB, GABA permease, anion: H+ symporter DASS family, 1-aminocyclopropane-1-acyteltransferase, lipopolysaccharide O-acyteltransferase, polysaccharides transporter PST family, lipopolysaccharides export system permease protein, lipopolysaccharides biosynthesis protein WzzE, acyl homoserine lactone synthase and acyl homoserine lactone efflux proteins, were highly enriched in phosphate rock samples compared to rhizospheric soil of native plants in phosphate mine and the rhizospheric soil of wheat crop. On the other hand, all of these bacterial functions are higher in the rhizospheric soil of native plants compared to the rhizospheric soil of wheat crops. On the other hand, the NhaC family and phosphate:Na+ symporter were enriched in the rhizospheric soil of native plants than in bulk phosphate samples.