Cave microbial community composition in oceanic islands: disentangling the effect of different colored mats in diversity patterns of Azorean lava caves

Processes determining diversity and composition of bacterial communities in island volcanic caves are still poorly understood. Here, we characterized colored microbial mats in 14 volcanic caves from two oceanic islands of the Azores using 16S rRNA gene sequences. Factors determining community diversity (α) and composition (β) were explored, namely colored mats, caves and islands, as well as environmental and chemical characteristics of caves. Additive partitioning of diversity using OTU occurrence showed a greater influence of β-diversity between islands and caves that may relate to differences in rare OTUs (singletons and doubletons) across scales. In contrast, Shannon diversity partitioning revealed the importance of the lowest hierarchical level (α diversity, colored mat), suggesting a dominance of cosmopolitan OTUs (>1%) in most samples. Cosmopolitan OTUs included members involved in nitrogen cycling, supporting the importance of this process in Azorean caves. Environmental and chemical conditions in caves did not show any significant relationship to OTU diversity and composition. The absence of clear differences between mat colors and across scales may be explained by (1) the geological youth of the cave system (cave communities have not had enough time to diverge) or/and (2) community convergence, as the result of selection pressure in extreme environments.

bacterial diversity, additive partitioning, lava cave, sampling grain, levels of organization, Azores

INTRODUCTION

Microorganisms living in subsurface ecosystems frequently enhance their survival by assembling into biofilms (Northup et al. 2003, 2011; Hall-Stoodley, Costerton and Stoodley 2004; Barton and Jurado 2007; Dapkevicius 2013; Jones, Schaperdoth and Macalady 2014). The presence and relative abundance of bacterial phyla in caves all around the world have been reviewed (Lee et al. 2012), with Proteobacteria, Chlorobi/Bacteroidetes, Actinobacteria, Acidobacteria, Nitrospirae and Chloroflexi being the most abundant groups. Within this broad pattern, subsurface microbial communities have shown predominant, non-random co-occurrence patterns, with ecologically similar taxa coexisting to a greater degree than expected by chance, indicating that deterministic processes are important in structuring communities (Stegen et al. 2012). Furthermore, cave-dwelling oligotrophic biofilms (microbial mats) are complex microecosystems, encompassing several evolutionarily distant bacterial phylotypes. Moreover, the relative abundance of phyla differs depending on the type of caves, pointing to a correlation between the main biochemical processes occurring in caves and the microorganisms they harbor. Therefore, the analysis of the composition and evenness of the phylotypes in cave microbial communities may give insight into the ecological and evolutionary context within the caves (Northup et al. 2003; Barton et al. 2007; Chen et al. 2009; Engel et al. 2010).

Caves are sheltered from atmospheric disturbances, and represent ecosystems in which many environmental variables remain relatively constant. Connectivity to the surface occurs mainly through entrances, skylights and rock cracks, the latter acting like narrow channels. Some cracks, particularly in lava caves close to the surface, originate from tree roots that penetrate into the cave, and can introduce rhizosphere microorganisms into the cave habitat (Snider 2010). Microbial mats and organic oozes cover the lava cave walls and ceilings with a range of colorful deposits, with white, yellow and tan being the dominant colors among the microbial mats (Northup et al. 2011, 2012). In spite of large geomorphological differences among caves, it is possible to identify mat types that look indistinguishable amongst lava caves. Color has recently been suggested as a proxy for the environmental conditions and biogeochemical processes operating in caves, as the distribution of different colored biofilms apparently can be influenced by nutrient availability and microclimate variations (Cuezva et al. 2009). Microbial mat distribution also varies depending on pH and the mineral association pattern, suggesting that preferential metabolic processes could be predominant in different colored biofilms (Cuezva et al. 2009). The link between biomineralization activities and bacterial diversity of communities has been described for different mineral transformations occurring in cave environments (Northup et al. 2003; Barton et al. 2007; de los Ríos et al. 2011; Miller et al. 2014).

Differences in bacterial communities according to mat color have been reported (Portillo, Gonzalez and Saiz-Jimenez 2008; Portillo and Gonzalez 2009). A comparative analysis of yellow colored spots along distant carbonate caves in Europe revealed a common core of microorganisms, some of them being close relatives to microbes recovered from yellow mats in lava cave walls (Porca et al. 2012). Also, in a pyrosequencing study of seven lava caves in Lava Beds National Monument (KH Lavoie, pers. comm.), mat color showed a trend to differentiate among operational taxonomic units (OTUs), but was not significantly different. Nevertheless, in a recent paper, Hathaway et al. (2014a) found no clustering pattern for mats according to their pigmentation in Azorean and Hawaiian lava caves. However, geographical location, as well as the levels of some elements in the rock composition (e.g. nitrogen, organic carbon, copper), was correlated with differences in the bacterial community composition.

This study provides the first broad, multiscale assessment of microbial diversity of lava caves in an isolated oceanic volcanic archipelago. The biogeographic isolation of the Azores, a remote archipelago in the North Atlantic Ocean, the numerous accessible lava caves and the vast bacterial diversity dwelling inside (Northup et al. 2011, 2012; Hathaway et al. 2014a,b) make Azorean islands an ideal testing ground for investigating the determinants of microbial diversity in volcanic caves. Three hierarchical levels of diversity with different isolation levels and distances between them were assessed, i.e. island (90 km), cave (c. 11 and 23 km for Terceira and Pico islands, respectively) and mat color (between 0.005 and 120 m). Linear mixed models (LMM) and β additive diversity partitioning (Crist et al. 2003) were used to quantify the variation among the three levels of organization. Furthermore, we investigated putative geographical, environmental and chemical factors that could explain the differences in OTU diversity and composition between caves.

Azorean caves host unique geological structures and endemic troglobitic macroorganisms (Borges et al. 2012), and understanding how diversity is structured in these isolated ecosystems is necessary for conserving species diversity at both micro- and macroscales. Studying how different components of bacterial diversity vary among different cave communities within and between islands can help in the selection of priority areas for conservation (Gering, Crist and Veech 2003).

MATERIALS AND METHODS

Area of study

The Azores archipelago is situated in the North Atlantic Ocean, along the Mid-Atlantic Ridge (approximately 36°55^′ to 39°43^′ N and 24°33^′ to 31°17^′ W). It is made up of nine islands harboring a total of 250 lava caves and volcanic pits described so far, with Pico and Terceira islands having the highest number with 118 and 73 natural caves each, respectively (Borges et al. 2012). Geologically, the islands are young, with ages of 3.52 Ma years (Terceira) and 300 thousand years (Pico) (Forjaz 2004; see also Borges and Hortal 2009).

Sampling site description and sample collection

A total of 13 volcanic caves were sampled, 9 from Terceira and 4 from Pico Islands (Fig. S1, Supporting Information), representing a variety of geological and physical conditions. Due to the length of Gruta das Torres and because its entrance divides the cave into two areas, an upper and a lower part, both were analyzed as independent caves (Table S1, Supporting Information). Three different mat colors were sampled per cave, namely yellow, tan and white. The three mat colorations were found in three of the different zonations described for volcanic caves, entrance, twilight and deep zones (Howarth 1993; Northup et al. 2012). The spatial distribution of the mat colorations showed no discernible common patterns among the studied caves. For each sample, wall rock chips covered with microbial mats were collected aseptically and covered with sucrose lysis buffer to preserve the DNA (Giovannoni et al. 1990), transported to the laboratory and stored at −80°C until DNA extraction. One soil sample was collected inside of each lava cave from the floor near the entrance and one soil sample was collected deeper within the lava cave, close to the microbial mats sampling sites. Exceptions for the soil sampling scheme mentioned above include Algar do Carvão and Gruta do Natal, where two samples were collected inside the cave (entrances are located inside a tourist facility), and Gruta da Terra Mole where no soil was found in the interior (Table S1, Supporting Information). Percentages of nitrogen (%N) and carbon (%C) in the soil samples were determined by high-temperature combustion. The resulting gases were eluted on a gas chromatography column and detected by thermal conductivity and integrated to yield carbon and nitrogen content. Samples for organic carbon analysis were placed in crucibles inside a desiccator with 6N HCl for 24 h prior to analysis to remove carbonate (Harris, Horwath and Van Kessel 2001). Analyses were performed on a ThermoQuest CE Instruments NC2100 Elemental Analyzer (ThermoQuest Italia S.p.A., Rodano, Italy) (Pella 1990a,b). One infiltrating water sample was collected per cave, with the exception of Gruta dos Balcões, Gruta dos Montanheiros and Gruta da Ribeira do Fundo where two water samples were collected. Results from the analyses of duplicated water samples were averaged (Table S1, Supporting Information). Concentrations of chloride (Cl⁻), bromide (Br⁻), sulfate (SO₄²⁻), nitrate (NO₃⁻) and phosphate (PO₄³⁻) were measured using a Dionex Ion Chromatograph DX-100 (Dionex, Sunnyvale, CA, USA) as described in Pfaff, Hautman and Munch (1997). The amount of ammonium (NH₄⁺) and total nitrogen (TN) in the water samples were analyzed using a Technicon AutoAnalyzer II (Technicon, Tarrytown, NY, USA). The limited water sampling was due to the scarcity of dripping water inside the caves. Samples were kept at 4°C until analysis. No soil or water samples were collected in Gruta da Achada because at the time of the collection, work on the trail precluded entrance into the cave (Table S1, Supporting Information). In order to avoid collinearity of chemical variables, which can lead to a misestimate of model parameters in subsequent regression analysis (see below), a principal component analysis was carried out including all log-transformed chemical variables (%C, %N, C:N, Cl⁻, Br⁻, SO₄²⁻, PO₄³⁻, NO₃⁻, NH₄⁺ and TN), and the three first PC axes that explained 82% of the variance were retained. The variance percentages and variable contributions to the axes are given in Table S2 (Supporting Information). For each of the 14 caves, temperature and humidity readings (wet bulb/dry bulb) were taken with an IMC Digital Thermometer probe and averaged (Table 1). Relative humidity is generally close to saturation once the deep cave zone was reached (i.e. cave walls were constantly covered by condensation water) and was therefore not considered. Additionally, geological age, where the value at the maximum boundary of the range was considered to proceed in a conservative manner (Calvert et al. 2006; Nunes et al. 2014), length and altitude of the caves (‘Os Montanheiros’ Association, pers. comm.) were compiled (Table 1). Surface precipitation was also recorded (Center of Climate, Meteorology and Global Changes, University of the Azores) although it was not considered in subsequent analysis because it was highly correlated with altitude (r = 0.85).

Table 1.

Names, location and main characteristics of the 14 caves sampled in the study. (a) Altitude was measured at the main entrance of the caves. (b) One sample of each mat color was collected from each cave.

			Age max.	Altitudea)	Length	Temperature	Sequences	Sequences	Sequences
Caves	Code	Island	(Years)	(m)	(m)	(°C)	yellow matb)	tan matb)	white matb)
Gruta do Lemos	GL	Pico	40000	15	28	17.1	75	65	70
Gruta dos Montanheiros	GM	Pico	1500	785	1805	11.0	79	46	71
Gruta das Ribeira do Fundo	GRF	Pico	50000	180	200	15.2	80	24	73
Gruta das Torres (Section 1)	GT	Pico	1000	300	2800	14.6	79	38	58
Gruta das Torres (Section 2)	GT2	Pico	1000	300	1000	14.6	79	22	63
Algar do Carvão	AC	Terceira	1730	585	90	11.3	64	24	68
Gruta das Agulhas	GA	Terceira	50000	1	250	22.2	30	28	45
Gruta da Achada	GAS	Terceira	1730	330	170	14.9	89	27	49
Gruta dos Buracos	GB	Terceira	7130	475	450	15.7	77	41	25
Gruta dos Balcões	GBL	Terceira	7130	422	4421	15.4	187	41	98
Gruta da Branca Opala	GBO	Terceira	7130	280	99	14.9	124	54	72
Gruta da Madre de Deus	GMD	Terceira	3000	59	245	14.6	59	101	86
Gruta do Natal	GN	Terceira	12861	551	697	15.8	90	72	111
Gruta da Terra Mole	GTM	Terceira	7130	387	120	14.9	62	40	20

DNA extraction and bacterial 16S rRNA gene clone library construction

DNA was extracted with a MoBio PowerSoil DNA Isolation Kit (MoBio, Carlsbad, CA) using manufacture's protocols with the exception of the substitution of bead beating for 1.5 min (Biospec Products, Bartlesville, OK, USA) instead of vortexing for cell lysis and elution in 30 μl C6 buffer. Nucleic acids were then extracted and purified in triplicate reactions from 0.5 g aliquots of the microbial samples. Bacterial 16S rRNA genes were amplified by PCR using the combination of universal primer 46F (5^′-GCYTAAYACATGCAAGTCG-3^′) and 1409 reverse (5^′-GTGACGGGCRGTGTGTRCAA-3^′) (Northup et al. 2010). Reactions were carried out in a 25 μL volume with 1X PCR buffer with 1.5 mM Mg²⁺, 0.4 μM of each primer, 0.25 mM of each dNTPs, 5 μg of 50 mg/mL BSA (Ambion Austin, TX, USA) and 1U AmpliTaq LD (Applied Biosystems, Foster City, CA, USA). The PCR reaction was performed with an MJ thermocycler with a program that consisted of preheating at 94°C for 5 min, 30 cycles at 94°C for 30 s, 57.8°C for 30 s, 72°C for 1.5 min and a final extension at 72°C for 10 min. The amplified products were purified using the QIAQUICK PCR purification kit (Qiagen, Germantown, Maryland). Bacterial 16S rRNA gene amplicons (ca. 1365 bases) were cloned using the TOPO TA Cloning kit (Invitrogen, Carlsbad, CA). Clone sequencing was carried out at the Washington University Genome Sequencing Facility using the T3 and T7 primers.

Sequence processing and OTUs delimitation

Sequences were checked for quality, edited and assembled with Sequencher 4.8. (Gene Codes Corporation, Ann Arbor, MI, USA). Sequences shorter than 1000 base pairs (bp) were removed from the data set. Sequence orientation was checked with OrientationChecker (http://www.softsea.com/review/OrientationChecker.html). Chimeras were detected using the Mallard/Pintail software (Ashelford et al. 2006). Alignment was accomplished with the SILVA Incremental Aligner (SINA) (Pruesse, Peplies and Glöckner 2012; http://www.arb-silva.de/aligner) and inspected visually. OTUs were assigned based on the uncorrected pairwise distances between aligned 16S rRNA gene sequences using the average neighbor clustering (Schloss and Westcott 2011). OTU delimitation was performed with the software package mothur 1.29.1 (Schloss et al. 2009), as well as the taxonomic assignment of the sequences performed by the Silva-based alignment tool using the default parameters.

Sequences were submitted to the NCBI GenBank database (http://www.ncbi.nlm.nih.gov/genbank/). Accession numbers are given in Supplementary Information SI1.

Statistical analysis

Our sampling unit was mat color (N = 42 samples) nested in caves (N = 14) and caves nested in islands (N = 2). All statistical analyses were implemented within the R statistical programming environment (R Core Team 2013). Details of the packages and functions used in the following analysis are given in Supplementary Information SI2.

Diversity indices

Sampling efficacy was assessed for each sample by calculating the ratio of observed OTU richness to several non-parametric richness estimators, namely Jackknife, Chao and the abundance-based coverage ACE. We tested for differences in OTU richness and Shannon Index among mat color, caves and between islands separately using Analysis of Variance (ANOVA) with Tukey's tests to identify differences among samples. ANOVA was implemented using LMM (Pinheiro et al.2014) with random factors to address the non-independence of samples with respect to colors, caves and islands (e.g. samples of the same caves share historic and geological conditions). Therefore, the effect of color was tested using caves nested in islands as random effect, the effect of caves was tested using mat color and islands as random effect and the effect of islands was tested using mat color and caves as random effect.

Cluster analysis

To compare community composition, OTU abundances were log-transformed to give lower weight to rare and possibly ‘erratic’ OTUs, and pairwise dissimilarities among samples were calculated using the relative abundance-based distance Bray–Curtis (BC). Non-metric multidimensional scaling (NMDS) analysis was performed on the dissimilarity matrices to visualize patterns of community composition. To test differences in OTU composition, we performed LMM following by Tukey's tests with adequate random structures (i.e. as implemented for diversity indices) using the three NMDS axis independently as response variables. To relate the resulting NMDS axes back to the original OTU matrix, we used vector overlays of the correlation between axes and OTU abundances using Rank Spearman correlations.

Diversity partitioning

In order to investigate changes in community structure across our three levels of organization i.e. mat colors, caves and islands, we used additive partitioning of diversity (e.g. Crist et al. 2003; Legendre, Borcard and Peres-Neto 2005). With this analysis, the OTU diversity at regional scale (total diversity combining Terceira and Pico islands) can be partitioned into components representing within-community diversity (α diversity) and among-community diversity (β diversity) (Fig. 1a). Using additive partitioning, we analyzed which levels contributed the most in explaining total OTU diversity considering both richness and Shannon diversity. The total observed OTU diversity can be therefore partitioned as

\begin{equation*} {{\rm{\gamma }}_{{\rm{obs}}}} = {\rm{ }}{{\rm{\alpha }}_{{\rm{colors}}}} + {{\rm{\beta }}_{{\rm{colors}}}} + {{\rm{\beta }}_{{\rm{caves}}}} + {{\rm{\beta }}_{{\rm{islands}}}} \end{equation*}

where α_colors is the mean α diversity per mat color, β_colors is the mean between-mat color β diversity within caves, β_caves is the mean between-caves β diversity within islands and β_islands is the between-islands β diversity. The statistical significance of level-specific α and β estimates was tested through a randomization procedure. We used a complete randomization (Crist et al. 2003) to generate 1000 random distributions of OTUs among samples at all hierarchical levels to generate the null distribution of each α and β estimates for each hierarchical level. Each original level-specific estimate was then compared with the appropriate null distribution and used to test the null hypothesis that the observed α and β diversities were obtained by a random distribution of OTUs among samples at all hierarchical levels. Statistical significance was assessed by the proportion of null values that were greater than (or smaller than) the actual estimate (Manly 1997).

Figure 1.

(a) Schematic representation of the three levels of organization considered to investigate microbial diversity. (b) OTU composition of the each sample (N = 42). The circles associated with each sample show the 16S reads classified by best matching bacterial class, with rare groups in the ‘other’ category. Sizes of the circles are proportional of the log-transformed rarefied abundances (See methods).

Relationship between diversity, community composition and environmental and chemical variables

As one single measure of each environmental/chemical variable was obtained for each cave, we tested the relationship between environmental variables, diversity indices and community composition by first pooling OTUs per cave and secondly by rerunning the following analyses for each mat color independently. The cave Gruta da Achada (GAS) at Terceira Island was excluded from this analysis because no environmental/chemical data were recorded. Relationships between environmental/chemical variables, OTU richness and Shannon Index were assessed by using LMM with islands as random effect while the relationships with community composition was assessed by using distance-based redundancy analyses (dbRDA). For both analysis, we first included all variables (i.e. island age, cave altitude, cave length, cave temperature and the three PC axes) in a global model and then simplified this model by sequentially removing the least significant variable until only significant ones remained. Because dbRDA does not allow the inclusion of random effects as it is implemented in mixed models, we performed dbRDA with and without islands as a covariate to evaluate the importance of island of origin in the analysis.

Scanning electron microscopy

Samples were examined on a JEOL 5800 SEM equipped with an energy dispersive X-ray analyzer, at high vacuum with an accelerating voltage of 15 KeV with a beam current between 0.1 and 0.01 ηA. Rock chips with microbial mats were mounted directly on scanning electron microscopy (SEM) sample stubs in the field, air-dried and coated with Au-Pd metal for imaging in the laboratory.

RESULTS

Descriptive statistics

We identified 2706 sequences corresponding to 872 OTUs (clustering sequences at the 97% similarity level) belonging to 20 different phyla (Acidobacteria, Actinobacteria, Armatimonadetes, Bacteroidetes, candidate division OD1, candidate division TM7, Chlamydiae, Chlorobi, Chloroflexi, Cyanobacteria, Elusimicrobia, Firmicutes, Gemmatimonadetes, Nitrospirae, Planctomycetes, Proteobacteria, TA06, Thermotogae, TM6 and Verrucomicrobia). Sequences that could not be assigned to taxonomic affiliations at the 97% similarity level (i.e. 10.75% of all sequences) were labeled as ‘unclassified’ and corresponded to 12.6% of the OTUs. Proteobacteria (sequences = 35.6%; OTUs = 14.3%), Actinobacteria (23.6%; 13.9%) Acidobacteria (11.8%; 13.3%) and Nitrospirae (9.3%; 12.2%) represented the majority of the cave phyla (Fig. 1b). Proteobacteria sequences included Alphaproteobacteria (36.3%), Betaproteobacteria (13.2%), Gammaproteobacteria (34.4%) and Deltaproteobacteria (16.1%). Among the 872 OTUs, 5 of them (0.6%) represent 33% of the total number of sequences found in our study and were recorded in most of the samples. More precisely, the most predominant OTU (OTU 866) belonged to the Crossiella genus, with 503 sequences (18.6%) in 39 of the 42 samples. The second most common OTU (OTU 870) belonged to the genus Nitrospira, with 142 sequences (5.2%) in 34 samples. In contrast, 601 OTUs (69%) were singletons (OTUs with only one sequence) and 118 (13.5%) were doubletons (OTUs with two sequences). Taxa with percentages higher than 1% represented abundant/cosmopolitan community members, while common and rare fractions were defined by 0.1–1% and <0.1% ranges, respectively.

OTU diversity and community composition among colors, caves and islands

The completeness values for each sample varied for each estimator but were in general low and varied greatly among samples: from 25.00 to 51.07% (mean = 36.20%) for Jackknife 1, from 9.50 to 64.00% (mean = 34.90%) for Chao1 and from 5.8 to 61.2% (mean = 28.0%) for ACE (See Table S3, Supporting Information). Therefore, to avoid biases arising from variation in completeness among samples, rarefaction per sample was used in all subsequent diversity analyses (diversity and distance matrices were calculated as averages of 1000 replicate subsamples of 20 reads per sample, see Supplementary Information SI3).

LMMs did not detect any difference between colors, caves or islands. OTU composition varied substantially among mat color and caves, but not between islands. Analyses on the axes NMDS 1 showed that two caves, Gruta da Achada (GAS) in Terceira island and Gruta dos Montanheiros (GM) in Pico island differed in OTU composition from other caves (LMM: F_13,23 = 4.268, P = 0.001). Post-hoc tests showed that GAS cave differed from Algar do Carvão (AC) and Gruta da Madre de Deus (GMD) caves, while GM cave differed from AC, Gruta das Agulhas (GA), GMD, Gruta do Lemos (GL) and Gruta das Torres (GT and GT2) caves (post-hoc test: positive difference, P < 0.05). The difference was mainly reflected in higher abundance of OTU 868 (Alphaproteobacteria), OTU 679 (Gemmatimonadetes), OTU 832 (Betaproteobacteria), OTU 629 (Acidobacteria) and OTU 870 (Nitrospirae), and a lower abundance of OTU 869 (Gammaproteobacteria) and OTU 435 (Alphaproteobacteria) (Table S4, Supporting Information). However, GAS and GM caves, on Terceira and Pico islands, respectively, did not differ between each other (post-hoc test P = 0.99). Analysis on the NMDS 2 also showed significant differences among caves but only for observed data (LMM: F_13,23 = 1.806, P = 0.104). However, NMDS 2 significantly discriminated tan from white mat color (LMM: F_2,26 = 5.801, P = 0.008; post-hoc test: negative difference, P < 0.05) (Fig. 2), with samples from tan mats displaying a lower abundance of OTU 870 (Nitrospirae), OTU 843 (Nitrospirae), OTU 807 (Nitrospirae), OTU 851 (Acidobacteria) and OTU 701 (Alphaproteobacteria) compared to white colored mats (Table S5 for details, Supporting Information). Analysis on the NMDS 3 did not reveal any differences among mat color, caves or between islands.

Figure 2.

NMDS ordination based on BC dissimilarities of OTU composition (stress = 0.15). NMDS was performed with three dimensions and only the first two are represented. Each level was highlighted separately for clarity with (a), (b) and (c) for colors, caves and islands, respectively.

Diversity partitioning

Additive partitioning analysis revealed differences between OTU richness and Shannon Indices across levels of organization (Fig. 3 and Table 2). For richness data, β_islands and β_caves contributed the most to the total richness (42.3% and 46.9%, respectively), although none were significantly greater or smaller than expected by chance (Table 2). Only α_colors was significantly higher than excepted by chance. The Shannon Index, within mat color (α_colors), contributed the most to the total diversity (56.8%) and was significantly smaller than expected by chance. Additionally, both β_colors and β_caves, contributing respectively 15.1% and 23.3% of the total richness, were significantly greater than expected by chance. In contrast, β_islands explained only 4.7% of the total richness and was not significantly different from the random expectation.

Figure 3.

Additive partitioning for OTU richness and Shannon Index observed at all scales. Results are expressed as percentages of γ diversity. Results of the statistical tests are given in Table 2.

Table 2.

Additive partitioning of α and β OTU diversity considering both richness and Shannon index. Observed partitions (Obs) are compared with those expected (Exp) after 1000 randomizations. Observed values are significantly higher or lower than expected at P < 0.05 or P > 0.95, respectively; otherwise they are not significant (ns).

	Observed	Expected	P-values
Richness
αcolors	34.476	32.723	<0.001
βcolors	58.952	59.380	0.346 (ns)
βcaves	409.571	410.084	0.864 (ns)
βislands	369	369.812	0.864 (ns)
Shannon
αcolors	2.973	3.370	<0.001
βcolors	0.791	0.727	<0.001
βcaves	1.218	0.899	<0.001
βislands	0.244	0.231	0.258 (ns)

Effects of environmental and geographical variables on OTU diversity

OTU richness and Shannon index calculated for both pooled samples per cave (three samples) and per mat color were not significantly explained by any variables (LMM). The only exception was for OTU richness for pooled samples, which was significantly explained by cave temperature (positive relationship, likelihood ratio test P = 0.003) and altitude (positive relationship, likelihood ratio test P = 0.001).

SEM for mat color samples

SEM revealed that all three colors, white, tan and yellow-colored samples, contained some similar predominant morphologies. These included filaments that were mostly covered with hair-like and knobby extensions (Fig. 4B–I). Biofilm was also a common feature in the SEMs of all three colors, but the iron oxide biofilm (Fig. 4A and H) only showed up in a limited number of samples. In some samples, the colonies were observed to be erupting from the biofilm (e.g. Fig. 4C, D and F). Smooth filaments were more rare, but did occur in all three colors of samples (Fig. 4B–D and H). Beads-on-a-string morphologies have been observed in all three colors of samples (e.g. Fig. 4C and D). Smooth coccoid morphologies were mainly observed in yellow-colored samples. Coccoid morphologies with hair-like extensions were observed mainly in yellow-colored samples, while white samples were seen to contain coccoid morphologies with knobby extensions (e.g. Fig. 4B and Hathaway et al. 2014a). Rod-shaped morphologies were observed in tan and white-colored samples (e.g. Fig. 4E and Hathaway et al. 2014a).

Figure 4.

SEM of white (upper row), tan (middle row) and yellow (lower row) microbial mats from Azorean lava caves.

DISCUSSION

This study represents the most intense cave microbial sampling effort done in the archipelago of the Azores. The predominance of Proteobacteria is in accordance with the findings of earlier studies in lava caves of this archipelago (Hathaway et al. 2014a). Alpha- and Gammaproteobacteria account for more than 70% of the total sequences (36.3% and 34.3%, respectively), while Delta- and Betaproteobacteria represent 16.1% and 13.2%, respectively. Actinobacteria appear with higher percentages of sequence recovery than Acidobacteria, in agreement with the percentages found in Hawai'i lava caves, but in contrast with previously described results for Terceira island (Hathaway et al. 2014a). Proteobacteria, Actinobacteria, Acidobacteria, Nitrospirae and Chloroflexi were the only phyla recovered from all caves. These dominant phyla are similar to the ones described in other cave systems (see Barton et al. 2007; Portillo, Gonzalez and Saiz-Jimenez 2008; Portillo and Gonzalez 2009; Portillo, Saiz-Jimenez and Gonzalez 2009; Pašić et al. 2010; Northup et al. 2011; Porca et al. 2012; Barton et al. 2014; Hathaway et al. 2014a), excluding sulfur caves, which present substantially different phyla content (Engel et al. 2010; Jones, Schaperdoth and Macalady 2014).

The analysis of diversity patterns across our three levels of organization showed that OTU diversity between mat color, caves and islands did not vary (Table S3, Supporting Information) and that only β components were informative in distinguishing the different scale of organization considered (Fig. 3, Table 2). For OTU richness, β_caves and β_islands were found to contribute the most to the total diversity, although these did not differ significantly from null expectation. This pattern may be explained by the presence of singletons and doubletons in our samples, inflating the role of higher levels of organization in their contribution to the total diversity. Hathaway et al. (2014a) reported the influence of geographical isolation in speciation comparing diversity patterns between Azores and Hawai'i lava cave microbial communities. Furthermore, the influence of geographical isolation on bacterial diversity supports the hypothesis that dispersal limitation plays an important role in determining assemblage structure of microbial communities (Fontaneto and Hortal 2012).

In contrast, the additive partitioning performed with Shannon index revealed that the majority of OTU total diversity was contained at the smallest level, i.e. in mat color (α_colors; 60% of the total diversity). Unlike what was found for richness, rare OTUs have little effect on the Shannon Index; therefore, the Shannon Index should reflect the high abundance of cosmopolitan OTUs dwelling in most of the caves and mats analyzed. The presence of these abundant widespread OTUs across caves and mats may be a consequence of the age of the studied caves, too young to have hosted extensive divergence events, or at least for the new endemic OTUs to dominate the OTU assemblages. Furthermore, Azorean islands do not present pronounced environmental gradients, therefore providing fewer opportunities for adaptive radiation. Such arguments have been put forward to explain why the Azores hosts so few endemic plants and arthropods (e.g. Borges and Hortal 2009). Although we could not draw conclusions from larger organisms’ ecological theory, some relationships between macro- and microorganisms have been found to be comparable with variations in the rates of the processes (Bell et al. 2005; Green and Bohannan 2006; reviewed by Soininen 2012). Another possible explanation is given through the analysis of species associations. It has been noted that some bacterial taxa tend to aggregate, with a more frequent occurrence for cosmopolitan taxa and for phylogenetically related taxa (Pascual-Garcia, Tamames and Bastolla 2014). It is likely that both habitat filtering (abiotic factor) and ecological interactions (biotic factor) could have shaped the current diversity pattern.

Concerning OTU composition, few differences among caves and among mat colors were observed (Fig. 2). The geographic distance between GM and GAS caves (located in different islands) and the relatively different environmental variables recorded provided no easy ecological explanation for the similarity between these two caves. OTUs responsible for the observed shifting of white and tan mat samples along the second axes of the NMDS could relate to nitrogen metabolism, a shared trait between some of the members of the order Nitrospirales and Rhizobiales (Alphaproteobacteria), but further studies should be conducted to specifically explore this issue. Overall, the absence of strong divergence in OTU composition between mat color, caves and islands when abundance was considered might be reflecting a similar ‘core’ of community functioning. Crossiella and Nitrospira were the only two genera present in all studied caves. Crossiella, an aerobic, non-motile actinomycete, nitrate-reducing genus, was by far the most abundant and ubiquitous found in Azorean lava caves and it has already been described in subsurface environments (Barton et al. 2007; Stomeo et al. 2008; Portillo, Alloza and Gonzalez 2009). Nitrospira is a chemolithoautotrophic nitrifying genus described in many caves (Vlasceanu et al. 2000; Schabereiter-Gurtner et al. 2004; Zhou et al. 2007; Chen et al. 2009; Porca et al. 2012; Hathaway et al. 2014a,b). Other cosmopolitan taxa found in our study include members involved in the nitrogen cycle, supporting the importance of this cycle in Azorean lava caves (Hathaway et al. 2014a). The Beijerinckiaceae family presents chemoorganotrophic lifestyles, i.e. N₂-fixing and methanotrophic activities and the Nitrosococcus genus oxidizes ammonia compounds to nitrites. Both taxa have already been described in caves (see van de Kamp 2004 and Barton et al. 2014 for the Beijerinckiaceae family and Schabereiter-Gurtner et al. 2004 and Chen et al. 2009 for the Nitrosococcus genus). An excess in nitrogen concentration in Azorean soils has been reported, due mainly to the animal production industry in the archipelago (Fontes et al. 2011), which might indicate that substantial amounts of nitrogen could be percolating into lava caves. Because of their abundance, these cosmopolitan OTUs could be regarded as important taxa for community functioning (Zhang et al. 2013). Investigating their relative abundance is a potentially important research topic to identify key ecosystem processes in cave environments.

Caves are characterized by zonal environments according to the distance to entrances (Poulson and White 1969; Howarth, 1983, 1993), passage geometry and microenvironments that result from several types of reactions, including microbial processes that often involve redox reactions (Barton and Northup 2007). Based on the environmental data collected for this study, we conclude that the variables analyzed were not the main drivers of the diversity patterns observed. However, the lack of correlation between environmental variables and bacterial diversity may be the result of our not including the relevant environmental descriptors. By not being able to properly characterize the microenvironments for each color mat, due to the difficulties of performing such measurements, we may have missed the environmental factors that truly determine bacterial diversity in the study caves.

In Azorean lava caves, no bacterial taxa were significantly associated with mat color, as opposed to what has previously been suggested for microbial mats in carbonate caves (Portillo, Gonzalez and Saiz-Jimenez 2008; Portillo and Gonzalez 2009; Porca et al. 2012). For instance, microbial communities from yellow mats of three carbonate caves were suggested to share a core microbial community composed of phylotypes affiliated with the Pseudonocardiaceae family, Chromatiales order and the genus Nitrospira (Porca et al. 2012). However, in our study, representatives of the above-mentioned groups (OTU 866, 869 and 870, respectively) were recovered in all three colored mats. The Crossiella genus (representing the Pseudonocardiaceae family in 96.8% of the sequences) was recovered from almost every sample and in every colored mat (OTU 866). Nitrospira and Chromatiales (99.5% of the sequences belonged to Nitrosococcus genus) also appeared in all three mat colors (OTU 870 and OTU 869). Several members of the Pseudonocardiaceae family showed pigment diffusion in culture and their possible contribution to the observed mat coloration has been mentioned (Porca et al. 2012). Other phylotypes have been suggested as contributors to the color appearance of yellow mats because of carotenoid production, such as members of the Xanthomonadales order (Portillo, Gonzalez and Saiz-Jimenez 2008) and the Steroidobacter genus (Porca et al. 2012). In our survey, Xanthomonadales and Steroidobacter were better represented in white mats compared to yellow and tan mats. Because active members of communities may exert the most relevant influence in color (Portillo, Saiz-Jimenez and Gonzalez 2009), differences in metabolic activity have also been explored to explain differences in mat color. Also, the presence of a larger proportion of active Desulfovibrio-related sulfate-reducing bacteria cells in yellow biofilms compared to biofilms with other colorations has been shown (Portillo, Gonzalez and Saiz-Jimenez 2008). Most sulfate-reducing bacteria are phylogenetically placed within the Deltaproteobacteria class, including the genus Desulfovibrio. We did not observe differences in the presence of Deltaproteobacteria among the diverse mat colors, although our results did not distinguish between metabolically active microorganisms and dormant or nonviable ones. Also, our SEM results suggested that many similar looking bacteria were found among yellow, tan and white mats.

Our study provides insight into lava cave microbial communities; however, rarefaction analysis indicated that the volcanic cave environment was not sampled to saturation (Fig. S2, Supporting Information) and therefore a more comprehensive sampling is required to provide a more complete assessment of these microbial communities. A survey of archaea, protists and fungi, which could also be contributing or determining the color of the mats, should be performed to have a more complete view of the microbial communities in volcanic caves. The presence of fungi, yeasts and slime molds in caves has been recently reviewed for carbonate caves but unfortunately, no data for lava caves were included in the study (Vanderwolf et al. 2013). Whether microbial mat color is determined by phylotype occurrence or a combination of traits in a community should be considered. Convergent evolution and horizontal gene transfer can result in the distribution of traits across multiple phylogenetic groups (Snel, Bork and Huynen 2002), the latter occurring in a more efficient way in biofilms (Molin and Tolker-Nielsen 2003). Several recent studies in soil microbial communities are pointing to the fact that quantifying diversity at high taxonomical levels could provide critical information on shared specific ecological traits (Fierer, Bradford and Jackson 2007; Philippot et al. 2010; Lennon et al. 2012). Indeed, it has been shown that complex functional traits are shared among members of deep clades of microorganisms due to phylogenetic conservatism (Martiny, Treseder and Pusch 2013). Furthermore, our results are based on OTU clustering, with the widely accepted assumption that 16S rRNA gene sequence similarity is a good proxy for ecological similarity (Schloss and Westcott 2011, but see Tikhonov and Wingreen 2015). However, identical 16S rRNA gene sequences may not necessarily correspond to identical genomes or ecotypes (Jaspers and Overmann 2004). Differences in biomineralization processes could also play an important role in the color appearance of the microbial mats. Different CaCO₃ deposits have been reported to predominate in diverse mat colors in carbonate caves (Cuezva et al. 2009); whether this fact is also occurring in volcanic caves needs to be studied. Biologically driven processes involving silica have been described in volcanic caves (Kashima, Irie and Kinoshita 1987), among them the precipitation of silica moonmilk, formed by the weathering of basaltic rocks (Onac and Forti 2011).

In conclusion, we found that bacterial diversity did not vary significantly among mat color, caves and between islands. Our results suggest that the color of the mat is not dictated by the abundant phylotypes shared across mat colors. Neither singletons nor doubletons are responsible for the mat color, as they are unique for each mat. Our analysis to detect differences in cave bacterial communities showed that microbial communities at the three levels of organization investigated, i.e., microbial mats, caves and island, were similar for the most dominant OTUs, but not identical when rare microbes were taken into account.

The studied lava caves have limited connectivity to the surface and occur on islands. These two levels of isolation may be reflected in the relevant contribution of β_islands and β_caves in richness values, with high percentages of singletons. The maintenance of higher levels of rare taxa in oligotrophic systems (Jones and Lennon 2010) has been confirmed in our study of Azorean lava caves. Genetic diversity could be the key to overcoming severe environmental conditions and a reversible dormant state could have implications in the maintenance of biodiversity (Lennon and Jones 2011). At the same time, our data suggest a set of communities with a low number of ubiquitous OTUs, which may be due to (i) the geological youth of the study system, i.e. cave communities have not had enough time to diverge; or (ii) community convergence, as the result of selection pressures in extreme environments (Kunin et al. 2008). These clusters could be sharing basic functional activities for the fitness of the community. A finer taxonomic resolution approach combined with a metatranscriptome analysis would shed light on this question.

Our results show that bacterial communities inhabiting Azorean lava caves are composed of consortia of cosmopolitan and more niche-specific members. This structure fits with what could be a general case for bacterial and archaeal communities (Tamames et al. 2010). Although, the evolutionary mechanisms that generate and maintain diversity of microbes in extreme environments are complex, additive partitioning can contribute to disentangling the input of the diversity components to the observed biogeographic patterns.

SUPPLEMENTARY DATA

Supplementary Data.

The authors wish to thank the Montanheiros association (http://www.montanheiros.com), especially Fernando Pereira, Maria João Leal and Paulino Costa, and Airidas Dapkevicius, Kenneth Ingham, Pedro Cardoso, Ana Rita Varela, Clara Gaspar, Ali Ghadimi, John Craig and Guida Pires for help during field and lab work. Island maps and surface precipitation data were kindly provided by the Center of Climate, Meteorology and Global Changes of the University of the Azores/Projecto CLIMAAT_Açores e ESTRAMAR (MAC/3/C177) supported by the European Union through the MAC Transnational Program of Cooperation – Madeira-Azores-Canaries. We acknowledge support from the UNM Molecular Biology Facility, which is funded by NIH (grant no. P20GM103452). The authors would like to thank the referees for their valuable comments which helped to improve the manuscript.

FUNDING

This work was supported by Portuguese national funds from the Foundation for Science and Technology of the Portuguese Government, [Understanding Underground Biodiversity: Studies in Azorean Lava Tubes (reference PTDC/AMB/70801/2006]. CR and IRA were funded by the Regional Fund for Science and Technology and Pro-Emprego program of the Regional Government of the Azores, Portugal [M3.1.7/F/013/2011 and M3.1.7/F/030/2011, respectively]. FR and PB were funded by the Foundation for Science and Technology of the Portuguese Government [FCT-PTDC/BIA-BIC/119255/2010]. This work was also supported by the Cave Conservancy of the Virginias, the Graduate Research Allocation Committee at UNM Biology, UNM Biology Grove Scholarship, the Student Research Allocation Committee at UNM, the National Speleological Society, the New Mexico Space Grant Consortium, and Kenneth Ingham Consulting (USA).

Conflict of interest. None declared.