The selectivity of 16S primers influences the results of the diversity of cyanobacteria from small water bodies

1. Introduction

Advanced sequencing technologies and metagenomic-based techniques, which appeared during the last decade have created great opportunities for the exploration of the taxonomic and functional diversity of microorganisms on Earth. However, many authors point out that primer specificity and sensitivity but universality on the other hand is key to the effective amplification of targeted Bacteria and Archaea [1,2,3,4,5,6,7,8,9]. Our studies of benthic [10,11] and soil-inhabiting [12,13,14,15] as well as endolithic microbial communities [10,11,16,17] in mountain desert ecosystems, indicated that when using universal primers designed by Klindworth and co-authors [2] targeting the V3-V4 region of 16S rRNA primers some information concerning representative genera of cyanobacteria occurring in these communities can be lost. To our knowledge, up to now most of the primers targeting the 16S rRNA gene and the hypervariable regions have been tested to investigate cyanobacteria from pelagic habitats selecting phytoplanktonic, often bloom-forming, and potentially toxic cyanobacteria [10,11,12,13,17,18,19,20,21]. Besides studies of planktonic cyanobacteria, however, investigations of benthic cyanobacteria become more relevant due to the presence of potentially toxic cyanobacteria [22,23].

Microbial mats are often dominated by cyanobacteria that thrive in water in a benthic community, at the border between aquatic and terrestrial environments, and even in arid places. Cyanobacteria prevailing in such communities may produce cyanotoxins, thus posing a health risk for humans and animals using that water [14,18,22,23,24,25]. Advanced molecular techniques based on good quality long reads are still costly, and the analyses take a longer time, while amplicon-based identification is much faster belonging to low-cost methods. Apart from the advantages, there are some challenges associated with fast amplicon-based analysis. One of them is the presence of the chemical substances (i.e., humic compounds and exopolymers produced by microorganisms) in the environmental samples that can be considered the inhibitors influencing DNA extraction. The presence of the inhibitors can also influence the genes’ amplification, thus complicating metabarcoding-based studies. Microbial mats, soil-inhabiting, and endolithic communities, which are characterized by a much higher density of the environment than planktonic samples, may potentially accumulate these compounds which may, even more, hinder the molecular analyses. For this reason, special reagents and DNA extraction kits are suggested for such analyses [26]. However, these communities are also largely inhabited by different sets of organisms than planktonic environments and it is important to find primers that will amplify products in such difficult circumstances. As has been reported, microbial mats from Eastern Pamir are characterized by vast morphological and taxonomic diversities [11]; thus, the present study was performed to test and optimize amplicon-based identification of cyanobacteria in exemplary natural microbial mats’ communities growing under heterogeneous conditions in the extreme environment such as the cold mountain desert of Eastern Pamir. Additionally, a mock community composed of various benthic cyanobacteria isolated from the studied environment as well as some planktonic, well-described cyanobacteria were chosen for the study.

We evaluated the amplification of three primer pairs; primers designed by Klindworth and co-authors (2013) targeting the V3-V4 hypervariable region of the 16S rRNA gene, primers designed by Lee et al. (2017) [27], targeting the V6 region of 16S rRNA gene and primers designed in the present study. Klindworth and co-authors [2] described primer pairs, which are specific to the Bacteria domain and are widely used to investigate bacterial phyla present in different environments [11,27]. The primers proposed by Lee et al. (2017) [27] were characterized by the authors as cyanobacteria-specific and were used for water quality monitoring in lakes and ponds [27] focusing mainly on representatives of heterocystous cyanobacteria (Nostocales). The third primer pair was designed for this study and was supposed to target a wide range of cyanobacteria taxa, including specific genera and species occurring in microbial mats in small water bodies in highly elevated areas with limited water resources. Thus, the aims of the present study were to search for optimal primers to study cyanobacteria from microbial mats that can be a source of toxins in small water bodies. To do so we aimed to: 1) design a new primer pair (C_773F - C_1083R) and 2) test in silico and in vitro three primer sets: newly designed primers, the primers that were described in the literature as cyanobacteria-specific (V6 - 1328F - 1664R) as well as those universal bacterial primers used in our previous studies (V3-V4- S-D-Bact-0341-b-S-17- S-D-Bact-0785-a-A-21).

2. Materials and Methods

2.1. Study Design

The presented experiment consisted of six consecutive steps from sampling and morphological identification of taxa, through the preparation of mock community (MIX) and isolation of DNA, primers testing in silico, and amplification of given 16S rDNA regions in studied microbial mats, followed by bioinformatic and statistical analyses (Figure 1)

The in silico and in vitro experiments were performed to test the specificity of three pairs of primers: S-D-Bact-0341-b-S-17 - S-D-Bact-0785-a-A-21 (V3-V4) [2], 1328F - 1664R (V6) [27] and C_773F - C_1083R (V4-V6) (this paper) in order to explore the diversity of benthic microbial communities mainly represented by cyanobacteria according to our previous study [10]. Primer names, target regions, sequences, and references are given in Table 1. To make it easier to follow, in the rest of the text, we use the targeted regions instead of long primer names. Each primer pair was applied to investigate the diversity and relative abundance of bacterial communities, focusing mainly on Cyanobacteria.

2.2. Sampling and Samples

For the analysis, we used a mock community as well as three environmental samples of microbial mats collected in August 2017 in Eastern Pamir from three water bodies with different types of mats. They were: (i) an unlayered type of mat from a hot spring used for bathing - Cyx 15, (ii) a multilayer soft mat from a thermokarst pond potentially used as a source of drinking water for wild and domestic animals - Cyx8 and (iii) shallow pool with an unlayered mat beneath the soil - Cyx9 (Table 2). The types of mats are described in Khomutovska et al. [11] and Sandzewicz et al. [12].

Mock community design. The mock community (further addressed also as “MIX”) consists of the DNA isolated from fifteen cyanobacterial strains including those from culture collection (one SAG), two planktonic culture collections of Institute of Oceanography (University of Gdansk), as well as strains isolated from environmental samples collected in Eastern Pamir (Table S1). The cultures from Pamir included representative cyanobacterial genera that were common in studied microbial mats.

2.2. Environmental Samples

The set of three environmental samples was selected based on the macroscopic features, and the morphology of microbial mat and taxonomic composition of communities was investigated based on optical microscopy analysis Nikon ECLIPSE Ni with digital camera Nikon DS-Fi1c. As mentioned above they represented three types of mats: multilayer soft, unlayered beneath the soil and unlayered [12]. The samples contained different representatives of Nostocales, Oscillatoriales, Synechococcales, and Chroococcales. The genera identified in the studied environmental samples are listed below. Investigated microbial mats differed in morphology, the content of EPS (extracellular polymeric substances), and secondary metabolites (cyanotoxins) [12,28]. While in Cyx9 and Cyx15 no cyanotoxins were detected, the Cyx8 mat was characterized by the presence of dmMC-LR (m/z 981.5).

2.3. DNA Extraction

The cultures for the mock community were grown in liquid WC [29] and BG11[30] culture medium. The biomass of the cultures was centrifuged, next the pellet was discarded; about 200 mg of wet biomass was homogenized and frozen. From 150 to 250 mg of dried microbial mats were homogenized in sterile tubes. For each sample, the DNA extraction was performed at least 3 times, and then the subsamples were mixed into one sample. For the DNA extraction of both cultures and environmental samples, the Soil DNA Purification Kit (GeneMATRIX, EURx Ltd., Gdansk, Poland) was used. The genomic DNA was isolated according to the manufacturers’ protocol of EURx (Gdansk, Poland). The concentration and quality of the DNA were checked using a Hybrid Multi-Mode Reader (Synergy H1, BioTek, Bad Friedrichshall, Germany).

2.4. Primer Design and Characterization

The analysis included the whole set of cyanobacterial 16S rDNA sequences (352) publicly available at the National Center for Biotechnology Information (NCBI) website and the Silva database. The sequences from these two databases were compared using FSA [31] and the most conserved fragment DNA for cyanobacteria was identified (Figure 2). Cyanobacteria-specific primers were designed using the Prime Design DNASTAR Lasergene v. 15, and AliView [32] tools. In silico analysis was performed using SILVAngs TestPrimer1.0 tool/pipeline.

2.5. Sequencing and Bioinformatics

The next-generation sequencing was performed by Genomed Joint-Stock company (http://www.genomed.pl/index.php/en) using the Illumina MiSeq platform and V2 MiSeq Reagent Kit. The demultiplexed paired-end reads were analyzed in QIIME 2 (2020.11) using the following plugins: “dada2” for sequencing filtering and clustering of the representative sequences, “demux” for the visualization of the quality of reads, “classify-sklearn” to classify representative sequences with SILVA 138 classifier (Table S2).

The Sanger sequencing 16S rDNA reads from the mock community were joined into contigs using SeqMan software (Lasergene v. 14.00, DNASTAR).

2.6. Statistical and Phylogenetic Analysis

The statistical analysis was performed using R packages 4.0 such as “vegan”, “tidyverse”, “pheatmap”, “devtools” and “factoextra” [R core]. The non-metric multidimensional analysis was carried out based on the Bray-Curtis dissimilarity matrix calculated for studied communities. The heatmaps were constructed based on the Spearman rank correlation matrix and Euclidean distance phylogenetic analysis and visualization.

Phylogenetic placement of cyanobacterial ASVs was conducted to verify the taxonomic position of the most abundant taxa of Cyanobacteria. The query sequences were aligned to reference sequences using the Cydrasil database [33]. The sequences were mapped on Cydrasil-ML-tree-bs 1000, and the Algaebase-based nomenclature was used [33,34].

3. Results

3.1. In Silico Testing of Designed and Commonly Used Primers

The result of the in silico, TestPrimer-based analysis indicated that the primers targeting V4-V6 region amplified mostly cyanobacterial sequences with the bacterial amplicons accounting only for about 25% of all sequences (Figure 2). The amplified sequences covered 72% of cyanobacteria and about 0.2% of other bacteria deposited in the Silva database. Amplifying the V3-V4 region, resulted in a Cyanobacteria share of about 1,1 % in the whole obtained sequences which covered 53% of cyanobacterial taxa in the Silva database. Additionally, the V4-V6 primer pair detected only 3 other phyla of the Bacteria domain (i.e., Bdellovibrionota, Firmicutes, Fusobacteriota) while the pair of V3-V4 primers detected representatives of 86 bacterial phyla, among which are Sericitochromatia, Vampirivibrionia, Firmicutes, Chloroflexi as well as some Archaea and eukaryotic algae (based on chloroplast DNA). The third pair of the tested primers, targeting the V6 region [24] identified only 10% of cyanobacterial sequences present in the Silva database although the second bacterial phylum, Bdellovibrionota accounted only for about 3% of all amplified sequences (Figure 2).

2.2. In Vitro Analysis

The QIIME2-based analysis of the reads revealed significant differences in the numbers of reads obtained for the samples (Table S1). Those differences were observed in the reads results for divergent primer pairs within the same sample as well as between the samples. The highest number of reads within all sample types were obtained for V4-V6 primer pairs (designed in this study), while the lowest number of reads in most cases were for V3-V4 primer pairs. Also, a high number of reads was obtained for primers amplifying the V6 region that was designed by Lee [27]. The tested primer pairs showed different results of reads that passed quality filters (Table S2.) which varied for the particular communities. The result of V6 amplicons showed the highest percentage of reads that passed quality filters for studied communities, from 90.3% (84 433) to 91.6% (83 944) in various samples, although the total number of reads that passed quality filters was, in most cases, the highest for the primers designed within the present study. The results of amplification of the V3-V4 region demonstrated the highest differences in the numbers of reads that passed quality filters. For this last primer pair, we obtained the lowest number of reads for particular samples.

Concerning the cyanobacterial reads, the results of V4-V6 and V6 primer pairs show that amplified communities were mostly represented by Cyanobacteria. Cyanobacterial reads accounted for 100% and 97% of all amplified ASVs for the V4-V6 and V6 regions respectively. The results of V3-V4 amplification revealed that the lowest number of cyanobacterial reads accounted for 4.33%, while the maximum reached only 67% (Table S3).

The results of the mock community amplification showed that the highest number of cyanobacterial reads was obtained for the V3-V4 and the V4-V6 amplicons. The V6 regions that were described as cyanobacteria-specific gave the lowest number of reads. The general picture for all studied primers could be described as follows: the largest numbers of cyanobacterial ASVs come from V4-V6 primers amplification analysis and a little less from the V6 amplification (Figure 3., Table S2).

Next, we analyzed the query (amplicon) cyanobacterial sequences using references of full-length 16S rRNA sequences and the reference phylogenetic tree. This analysis helped to verify the cyanobacterial sequences obtained for studied primers and samples and to identify them up to the genus (names prefixed on heatmap). Comparing results for the mock community (MIX) demonstrated that the primer pairs V4-V6 and V3-V4 identified identical genera i.e., Tychonema, Calothrix, and Geitlerinema, while V6 primers prevailed in the determination of distinct representatives of Nostoc, Calothrix and Tychonema (Figure 3). Moreover, in detailed analysis with mapping obtained sequences on reference cyanobacterial phylogenetic tree (Cydrasil tree), the V4-V6 demonstrated more precise identification, because within these amplicons it was possible to identify the highest numbers of reads up to genus or family level (Table S4). Additionally, although the number of identified cyanobacterial sequences in the case of V6 (about 72 000) was slightly higher than in the case of V4-V6 (about 60 000), half of the V6 reads were assigned to single sequence of Nostoc (over 36 000 reads).

Mapping the results of mock community on the phylogenetic Cydrasil tree allowed also for identification of the most strains used for the mock community after amplification with V4-V6 primers. Our primers showed a greater number of ASVs in specific cyanobacteria and gave no count of other bacterial readings in contrary to the other two primers (V3-V4 and V6 regions). The mapping of query sequences performed for the V6 region revealed that it was hard to compare this analysis to the 16S sequencing results of the L08 strain and unequivocally identify the presence of this strain in the MIX sample.

The structure of bacterial communities at the phylum level demonstrated distinct differences in the sequences amplified with V3-V4 primer pairs and the other two primer pairs, targeting cyanobacteria. The sample Cyx15a and Cyx9 grouped together demonstrating the prevalence of Pseudomonadota (Proteobacteria) Chloroflexi, Bacteroidota, Desulfobacteriota, and Actinomycetota (Actinobacteriota). The V3-V4 amplicons of the Cxy8 and MIX sample formed another clade characterized by the dominance of Cyanobacteria, followed by Pseudomonadota (Proteobacteria) and Bacteroidota. The sequences amplified with V4-V6 and V6 primer pairs grouped together with Cyanobacteria phylum determining the structure of the amplified communities and no visible variability within other phyla (Figure 4.).

The structure of the cyanobacterial community at the family level was grouped according to the mat type, not a primer (Figure 5), though in most cases the highest similarity was between the amplicons of V3-V4 and V4-V6 regions. In the Cyx15 sample, Leptolyngbyaceae and Limnotrichaceae families dominated, in Cyx8 Leptolyngbyaceae, Cyx9 Nodosilineaceae and in the mock community the structure was the most diverse, though with sequences identified as Phormidiaceae dominating the results. The results of all samples amplification with the V6 primers were characterized by lower specificity with a higher share of sequences identified only to the class or order level (Cyanobacteriia and Cyanobacteriales).

The ASVs results of the amplification with three different primers’ pairs were investigated further using the principal coordinate analysis (PCA) analysis (Figure 6 and Figure 7). The results showed variable clustering for the whole bacterial communities. The amplicons of V3-V4 of all samples clustered together in the middle of the PCA, exhibiting very low contribution to the DIM1 and DIM2. The amplicons of V4-V6 primer pairs of all samples had a higher contribution than the corresponding amplicons of V6 primer pairs. In the case of the MIX sample the results of all three primers were close to each other (Figure 6). The PCA analysis of Cyanobacteria phylum (at the ASV level) demonstrated that the majority of samples were grouped according to the type of environmental samples (Figure 7). Similarly, as in the results of the whole bacterial community we observed that in the case of cyanobacterial ASVs the highest contribution to the placement of various samples mostly had the amplicons of the V4-V6 primers. The one exception was the Cyx9 sample when the amplicon V6 had a similar contribution as the V4-V6. We also noticed that the MIX sample, in contrast to the analysis of the whole bacteria community, had the lowest contribution to the placement of the results. Again, there was one exception – in the Cyx8 sample the results of V3-V4 and V6 amplification were grouped together, with the MIX displaying very low contribution to the placement of sequenced communities. However, the sequences amplified with the V4-V6 primer pair were placed in the opposite quarter with the highest contribution.

The highest parameters of the diversity metrics were obtained using V3-V4 primers (Table 3). The sample Cyx9 was characterized by the highest Shannon index (7.83) for primers V3-V4, while the other two primer pairs scored much lower values. However, the V4–V6 amplicons exhibited a higher Shannon index than V6 amplification. The same pattern was present in other analyzed samples. Verifying the diversity for cyanobacteria phylum amplification with the V3-V4 primer pair again resulted in the highest value of the Shannon diversity index, with the V4-V6 having the second highest value. The differences between the V3-V4 and the other two primers were less pronounced. The number of cyanobacterial families identified with the V6 primer pair was generally lower (with one exception) than in the case of the other two primer pairs. The results concerning observed features of cyanobacteria at the ASV level were much more variable and did not repeat previous patterns.

4. Discussion

The results of the present study revealed that the primers designed within the present study amplifying fragment V4-V6 of 16S allowed for more precise identification of cyanobacteria in benthic microbial communities, compared to previously described primers [2,27] which were focused on planktonic communities. Additionally, the primers were tested experimentally with microbial mat samples from high altitude, arid environments such as small water bodies, that are less studied than the commonly investigated planktonic communities. Such types of samples are often claimed to be sources of novel microorganisms (including species of Cyanobacteria). Thus, using amplicon-based methods it is important to investigate the taxonomic composition of the microorganism at a higher level as it is often not possible to identify the species. Also, other researchers confirmed recently that there is a need to use different fragments of hypervariable 16S rRNA to study bacteria and archaea from less known environments [35]. The present study presents that the V4-V6 region of 16S rRNA seems to be a good solution to investigate diverse microbial benthic communities giving great resolution for Cyanobacteria, although not providing reliable results of the rest of the bacterial community. In turn the V3-V4 primer pair is more universal for all phyla of Bacteria including also some Archaea and eukaryotic algae (chloroplast DNA) but still provides considerably broad results for cyanobacteria. However, both results of the in silico as well as in vitro analyses showed that the V4-V6 primer pair allows for the detection of more cyanobacterial taxa (including representatives of orders Cyanobacteriales, Gloeobacterales, Leptolyngbyales, Phormidesmiales, Synechococcales, Thermosynechococcales) compared to the other studied primer pairs, such as V6, and V3-V4. The analysis of the MIX sample using mapping on the phylogenetic tree also demonstrated that the V4-V6 primers resulted in a much higher number of identified various ASVs than the amplicons of the other two primer pairs. The high number of V6 amplicons was actually misleading as Nostoc accounted for half of the reads. Interestingly also we were not able in the V6 amplicon to identify the strain L08, Hillbrichtia pamiria gen. nov. sp. nov. The reason for this can be the novelty of the L08 strain and the specificity of the V6 region concerning the identification of cyanobacterial orders/families. The PCA analysis of Cyanobacteria phylum (at the ASV level).

On the other hand, for the studies of the whole bacterial communities (including Cyanobacteria) the use of the V3-V4 primers would be advisable (Table 3.). Described comparison was made for the MIX sample, in which we used “unialgal” cultures of Cyanobacteria and although they were not axenic, the contribution of Bacteria was generally lower than in natural mat samples. Using primers designed within the present study (V4-V6) we obtained much higher number of reads of cyanobacterial amplicons, that were assigned to various taxonomic levels on the reference phylogenetic tree, without having amplicons of other bacteria as in the case of V3-V4 (Table S4). Such specificity of V4-V6 primers facilitates bioinformatics analyses. The primer pair amplifying the V6 region of the 16S rRNA gene gives fewer bacterial sequences than V3-V4, but less varied results of Cyanobacteria than V4-V6. This can be explained by their higher specificity than in the case of primers designed within this study. Our results revealed that microbial communities in which the dominated Nostocaceae genera can be explored more precisely while applying V6 primers as was noted by previous authors [27]. However, this comes at the expense of data on other cyanobacterial taxa, which have fewer reads, as more than half of the reads of the V6 amplicons we assigned to Nostoc. Also using V6 primers resulted with higher percentage of sequences identified to a class or possibly an order level only. We are aware the differences in the number of reads obtained for the individual taxa of microorganisms during the DNA amplification can be also related to the presence of different inhibitory substances in the analyzed environmental samples, which was also pointed out by other authors [9]. The environmental DNA can become fragmented during the isolation resulting in short sequences that in some cases can be identified up to order level only [36,37,38,39,40]. All these factors may negatively affect the results of the analyses and make their comparison difficult.

The PCA analysis of the whole Bacteria domain demonstrated that from among the three primer pairs used, the highest contribution to the PCA results of given samples had amplicons of the V4-V6 primers. This was further confirmed in the analysis of cyanobacterial reads only, where results showed an even more pronounced contribution of V4-V6. An exception was observed in samples of Cyx9 mat, in which V6 amplicons influenced the placement in a similar degree to V4-V6. The differences in the placement between various samples in PCA of Bacteria and Cyanobacteria amplicons are of course connected with the composition and structure of analyzed communities and the differences in the contribution of other bacterial phyla than Cyanobacteria. For example, in the case of the mock community, the results of V3-V4 amplification were close to the other two primer pairs. That’s possible because cyanobacterial strains used for the mock community (MIX), although not axenic, did not contain as much bacteria as the environmental samples. It was confirmed in the PCA for cyanobacteria only, in which results of all three primer pairs amplifications of MIX are close together. Regarding the Cyx15 sample, it grouped along the Y axis on the PCA of Cyanobacteria with distinctly the highest contribution of V4-V6 and more similar results in the case of V3-V4 and V6. The results may reflect the presence and predominance of Hilbrichtia pamiria gen. nov. sp. nov. in this mat [13]. It seems that the amplicons of both primer pairs, V3-V4 – the less specific as well V6 - the more specific towards Nostocales can be difficult to identify and may be classified to various taxonomic groups.

The primer pairs designed within the present study proved to amplify well and allowed for the identification, compared with other primers, of the most cyanobacterial taxa at the order level. However, our results also revealed that metabarcoding is still a limited method to study the composition and structure of cyanobacterial communities. Although the pair of primers designed in this study scored high, we are currently unable to determine unequivocally that it provides the optimal results for environmental analyses. Both the V4-V6 and V3-V4 primer pairs produced a large, though different, diversity, with the latter yielding the highest results. That may be due to chimeras and unspecific amplification of V3-V4 primer pair. Additionally, the V4-V6 primer pairs gave much fewer chimeras, and the higher number of cyanobacterial reads identified to family or even species. The results of the V4-V6-based amplification also gave higher diversity parameters than the results of V6 amplification, though the values were more comparable.

5. Conclusions

In our research, we have shown that choosing a suitable primer can be crucial for expected results. Primers amplifying the V3-V4 region of 16S rRNA can be applied to investigate the structure of the entire bacterial community. However, if a study aims to obtain as much information as possible about Cyanobacteria, including potentially toxigenic ones and from specific environment such as benthic communities, we suggest using V4-V6 primers, which are more specific for cyanobacteria than V3-V4 and less specific to a single taxonomic group, as is in the case of V6.