Characterization of the bacterial communities on recent Icelandic volcanic deposits of different ages

Site parameters and lava flow composition

The three different sampling locations represented deposits from different eruptions of the 1975–84 Krafla Fires (32-, 35- and 39-year old, respectively) and showed clear progression from poorly to fully vegetated sites. At each flow, three different points were chosen, resulting in three sites for the three different lava deposits. Two sites were in proximity to each other (10 m distance; designated as year and year*) while a third one (D) was at a distance of 100 m from the two other sites but still within the same lava deposit. Based on the meteoblue climate diagrams (www.meteoblue.com), the mean daily maximum and minimum temperature ranges between − 3° and 12 °C, and − 9 °C and 5 °C, respectively. The mean monthly precipitation ranges between 30 and 84 mm (as rain- and snowfall). The main wind direction is south. The 32-year flow is still exposed to the release of volcanic SO₂ and H₂S gasses and higher temperatures. No visible vegetation was observed on the 32-year flow. The 35-year flow showed lichen formation while the 39-year old flow already showed the formation of grasses and mosses (Fig. 1). The total number of cells determined by flow cytometry showed that the three sampling sites within the three flows contained between 6.80 × 10⁷ and 7.54 × 10⁸ cells/g (Fig. 2).

Basalt composition

Taxonomic distribution within each flow and comparison between the different flows and sites

Subsampling was performed based on the lowest amount of reads obtained over the nine different samples, i.e., 32-year sample with a coverage of 3060 reads (Fig. 3). Such a subsampling procedure allowed us to make a fair comparison between the diversity indices for all samples, e.g. avoiding a potential higher diversity due to a higher sequencing coverage.

Taxonomic profiling indicated that a large fraction of the operational taxonomic units (OTUs) could not be classified (Fig. 4; Additional file 1: Table S3). At the level of taxonomic families, sites 32, 32* and 32^D contained about 35%, 52% and 53% unclassified sequences, respectively. For the 35-year old flow, 35 and 35* contained 41 and 46% unclassified sequences while 35^D contained 66% unclassified sequences. The three sites within the 39-year old flow (39, 39* and 39^D) contained 40%, 53% and 47% unclassified sequences, respectively. An additional overview of the classified and unclassified number of raw reads, OTUs, single-read OTUs as well as OTUs contributing more than 1% to the taxonomic distribution is given in Table 3. Nevertheless, none of the OTUs identified in the samples were observed in the negative control, which was comprised of sterile water that was subjected to the same processing pipeline as the basalt samples. The major OTUs (relative abundance above 1%) constituted 1 to 18% of all OTUs, the single-read OTUs constituted 0 to 56% of all OTUs. The OTUs shared by the different sampling sites within each flow clearly differed (Fig. 5).

	32	32*	32D	35	35*	35D	39	39*	39D
Classified
# reads	1792	1725	1221	1600	1436	851	1695	881	1539
%	59.03	56.82	40.22	52.70	47.30	28.03	55.83	29.02	50.69
# OTUs	332	323	156	336	249	87	52	103	48
# single-read OTUs	170	148	29	149	103	23	0	30	3
# OTUs > 1%	8	11	3	6	7	5	10	5	9
Unclassified
# reads	1268	1335	1839	1460	1624	2209	1365	2179	1521
%	41.77	43.97	60.57	48.09	53.49	72.76	44.96	71.77	50.10
# OTUs	457	475	158	564	294	70	15	135	20
# single-read OTUs	263	295	28	298	124	9	0	43	1
# OTUs > 1%	2	4	7	1	8	18	2	9	3

In the 32-year old flow, the most abundant classified phylum across the three sites was Proteobacteria, although site 32 and 32^D were predominated by Acidobacteria (mainly Gp4) and Actinobacteria, respectively (Fig. 6). The most abundant OTU belonged to Gp4 (Acidobacteria; OTU9) for site 32 (6.7%) and 32* (4.3%), and was unclassified OTU2 (24.6%) for the more distant site (Fig. 7).

In the 35-year old flow, the most abundant classified phylum across the three sites also was Proteobacteria, although site 35 was predominated by Acidobacteria (Gp4 family) and site 35^D by Actinobacteria (Fig. 6). The most abundant OTU belonged to Gp4 (Acidobacteria; OTU9) for site 35 (7.3%), was unclassified OTU2 for site 35* (11.7%) and unclassified OTU5 for site 35^D (9.4%) (Fig. 7).

In the 39-year old flow, the most abundant classified phylum across the three sites was Firmicutes (mainly Clostridiaceae), although site 39* was predominated by Actinobacteria (mainly Geodermatophilaceae) (Fig. 6). The most abundant OTU in site 39 and 39^D was OTU5 (37.7% and 44.2%, respectively), and belonged to Geodermatophilaceae for site 39* (OTU10; 16.8%) (Fig. 7).

One hundred and five OTUs were shared between at least one sampling site of the 32-, 35- and 39-year old flow (Additional file 1: Table S3). Three were common to all samples and classified as Clostridium (OTU16), Propionibacterium (OTU22) and Bradyrhizobium (OTU48). Next to the OTUs shared by all three flows, additional OTUs were shared when flows were pairwise compared. The 32- and 35-year old flows additionally shared seven OTUs classified as Acidobacteria (OTU9), Actinobacteria (OTU65, OTU146 and OTU170), Proteobacteria (OTU67 and OTU79) and unclassified (OTU68). The 32- and 39-year old flows additionally shared one unclassified OTU (OTU302). The 35- and 39-year old flows additionally shared six OTUs classified as Bacteroidetes (OTU195), Firmicutes (OTU176), Proteobacteria (OTU113, OTU262 and OTU405) and unclassified (OTU5). The three most abundant OTUs in each sample and their abundance in the other samples indicated the high abundance of a few unclassified OTUs, especially in the 39-year old flow (Fig. 7).

Diversity and NMDS analysis between sites and comparison

Beta-diversity analysis (Yue & Clayton) indicated that the proximal 32-year old sampling sites as well as the 39 and 39^D sampling sites were highly similar (Fig. 8). Likewise, nonmetric multidimensional scaling (NMDS) analysis indicated that the proximal 32-year old sampling sites were similar as well as the 39 and 39^D sites, since these ordinated closer together (Fig. 9). The 35-year old sites do not cluster together, possibly indicating that differences for the 35-year flow are more pronounced than for the other flows. Twenty-eight OTUs significantly impacted the downward movement along NMDS axis 2 (with 12 being unclassified). Seventeen OTUs significantly impacted the leftward movement along NMDS axis 2 (with 8 being classified as Proteobacteria). In addition, age significantly (p = 0.026) influenced the leftward movement along NMDS axis 1.

Sample sites and sampling

Samples were taken from three different locations on the Leirhnjúkur fissures formed around Krafla, Iceland from the 1975–84 Krafla Fires. During this period, there was significant periodic rifting and faulting along the plate boundary, which was confined to a single system around Krafla. Both rifting episodes included a series of nine small (< 0.2 km³) effusive eruptions that occurred on fissures within the Krafla caldera volcano as well as on the nearby sectors of the fissure swarms [13]. On 27 July 2016, 9 sample sites were defined, three on each of the three lava flows (32-, 35- and 39-year old) originating from the fissure swarms (Fig. 1). The three flows were clearly different in vegetation ranging from almost none (youngest) to fully (oldest) vegetated and were relatively accessible for sampling. Within each flow, three sites were selected. Two sites were in proximity to each other (10 m distance; designated as year and year*) while a third one (^D) was at a distance of 100 m from the two other sites but still within the same lava deposit (Fig. 1). Coordinates of the 32-, 35- and 39-year old lava flows are 65°43’39.828’ ’N-16°47’23.3372’ ’W, 65°42’50.5332’ ’N-16°47’26.7576’ ’W and 65°42’56.8872’ ’N-16°47’31.3908’ ’W, respectively.

About 20 g of smaller sized rocks and soil material was sampled at 5 cm below the upper part of the deposits with a sterilized spoon and collected in a 50 ml sterile falcon tube (CELLSTAR® Centrifuge Tubes, Polypropylene, Sterile, Greiner Bio-One, Belgium). These tubes were sealed in three layers of plastic bags and transported at ambient temperature to the laboratory for analysis (see below). The time between collection and analysis was 7 days.

Composition analysis

Five grams of basalt rock material, with a mesh size below 200, were prepared and used to analyze the composition of the different basalt rocks, with lithium metaborate/tetraborate fusion and inductively coupled plasma optical emission spectrometry (ICP-OES)/inductively coupled plasma mass spectrometry (ICP-MS) analysis, and titration for FeO (Actlabs, Canada) (Tables 1 and 2).

Flow cytometry

Five grams of sample was added to 50 ml Butterfield’s Phosphate Buffer (0.25 M KH₂PO₄ in distilled water, adjusted to pH 7 with 1 M NaOH) and incubated at 30 °C in static conditions. After overnight incubation, part of the buffer solution was analysed with flow cytometry to analyse the total cell numbers (SYBR Green; SG) [38]. Samples were stained according to optimized procedures [39, 40]. Samples were prepared by diluting the sample (to obtain an event rate between 200 and 2000 events/μl) in 0.2 mm filtered mineral water. Next, SG dye (Sigma Aldrich; final concentration of 1×) was added and cell suspensions were incubated at 35 °C for 13 min. Stained bacterial suspensions were analysed on an Accuri C6 (BD, Erembodegem) with a blue (488 nm, 20 mW) and red (640 nm, 14.7 mW) laser, which was calibrated according to the manufacturer’s recommendation. Standard optical filters were used and SG was detected in FL-1 (530/30 nm) with the blue laser. A quality control with 6- and 8-peak fluorescent beads (by manufacturer BD, Erembodegem) and a cleaning cycle were performed prior to experiments to assess both the accuracy (bead count and position) and the cleanliness. Samples were analysed using the Accuri C6 software (version 1.0.264.21).

DNA extraction

To limit DNA extraction bias, two different protocols were used. DNA was extracted from 1 g basalt using the bead-beating cetyl trimethylammonium bromide (CTAB) phenol extraction protocol [41] and the FastPrep protocol [42] without the last Sepharose 4B based purification step [43]. Afterwards, extracted DNA was dissolved in 1xTE buffer and stored at -20 °C.

Illumina sequencing

The DNA concentration was quantified by the Quantifluor dsDNA sample kit on a multi-detection system (Promega, Leiden, the Netherlands). High-throughput amplicon sequencing of the V3-V4 hypervariable region [44] was performed with the Illumina MiSeq platform according to the manufacturer’s guidelines at LGC Genomics GmbH (Berlin, Germany). Sequences were preprocessed using the OCToPUS pipeline combining various preprocessing algorithms [45]. First, both forward and reverse reads are quality checked via looking at k-mer frequency to identify potential false k-mers using the Hammer algorithm implemented in the SPAdes tool. Next, the contigs were created by merging the paired-end reads using the heuristic based on the difference in Phred quality scores of both reads using an updated version of the mothur command “make.contigs” as described in [46]. These contigs were aligned to the SILVA database [47], followed by removing those contigs having an ambiguous base, homopolymer longer than 8 nt, or a length below 400, and those incompliant with the targeted region within the 16S rRNA gene using the mothur commands “align.seqs” and “screen.seqs”, respectively. Next, the aligned sequences were filtered and dereplicated using the mothur commands “filter.seqs” and “unique.seqs”. Secondly, sequencing errors were removed using the IPED algorithm - dedicated to denoise MiSeq amplicon sequencing data [45]. Chimera were detected using the de novo mode of the CATCh algorithm [46]. Creation of the operational taxonomic units (OTUs) was performed using UPARSE ([48] with default parameters (v7.0.1001_i86linux32 – commands sortbysize, cluster_otus, and usearch_global) [42, 45, 46, 49–52]. The datasets generated and analyzed during the current study are available in the NCBI Sequence Read Archive (SRA) repository (SRP106138).

Statistics

The datasets generated by the two DNA extraction methods were combined for further analyses. Subsampling was performed using the the mothur command “sub.sample” using the amount of reads of the smallest sample. Alpha diversity indices (Chao, Shannon, and inverse Simpson) were calculated using the mothur command “summary.single”. The beta-diversity between the samples was assessed using the Yue & Clayton measure of dissimilarity (mothur command “dist.shared”) and non-metric multidimensional scaling (NMDS) plots (mothur command “nmds”). For the NMDS plot, statistical significant correlation with the coordinates was indicated (p < 0.05) (mothur corr.axes command). The structure of the microbial populations was assessed using the analysis of molecular variance (mothur command “amova”) and homogeneity of molecular variance (mothur command “homova”). Statistical comparison of indices was performed using a one-way ANOVA analysis followed by a post-hoc Tukey test (including multiple testing).