(16.4 %, linear regression: P<0.001; Figure B.3), with diversity increasing from acidic to alkaline pH. Further multiple regression analysis showed NO3-, turbidity, ORP, dissolved oxygen, NO2-, Si, and Cd also had meaningful contributions (Table B.3). Cumulatively, along with pH, these factors accounted for 26.6 % of the observed variation in Shannon diversity.
Correlation of pH with Shannon index (Pearson’s coefficient: |r|=0.41, P<0.001) and significance testing between samples binned by pH increments (Kruskal-Wallis: H=179.4, P<0.001; Figure B.1) further confirmed pH as a major driver of variation in alpha diversity.
This finding is consistent with reports of pH as the primary environmental predictor of microbial diversity in several ecosystems, both in Aotearoa-New Zealand and globally (e.g.
soil91,227, water156,228, alpine68,229).
It has been previously hypothesised that pH has significant influence on microbial community composition because changes in proton gradients will drastically alter nutrient availability, metal solubility, or organic carbon characteristics91. Similarly, acidic pH will also reduce the number of taxa observed due to the decreased number that can physiologically tolerate these conditions191 compared to non-acidic habitats. Here, we demonstrate that pH had the most significant effect on diversity across all springs measured, but due to our high sampling frequency, we see this influence diminishes at temperatures >70
°C (Figure 2.2). Inversely, the effect of temperature on diversity was lessened in springs where pH was <4 (Figure B.5). There is some evidence that suggests thermophily predates acid tolerance191,230, thus it is possible the added stress of an extreme proton gradient across cell membranes has constrained the diversification of the thermophilic chemolithoautotrophic organisms common to these areas231. Indeed, a recent investigation of thermoacidophily in Archaea suggests hyperacidophily (growth pH <3.0) may have only arisen as little as ~0.8 Ga191, thereby limiting the opportunity for microbial diversification; an observation highlighted by the paucity of these microorganisms in extremely acidic geothermal ecosystems139,191. It is also important to note that salinity has previously been suggested as an important driver of microbial community diversity20,90. The quantitative data in this study showed only minimal influence of salinity (proxy as conductivity) on diversity (linear regression: R2=0.001, P=0.272; Table B.2), bearing in mind that the majority of the geothermal spring samples in this study had salinities substantially less than that of seawater.
Figure 2.3 - Constrained correspondence analysis (CCA) of beta diversity with significant physicochemistry. (a) A scatter plot of spring community dissimilarities (n=923), with letters corresponding to centroids from the model for geothermal fields (A-Q; White Island, Taheke, Tikitere, Rotorua, Waimangu, Waikite, Waiotapu, Te Kopia, Reporoa, Orakei Korako, Whangairorohea, Ohaaki, Ngatamariki, Rotokawa, Wairakei-Tauhara, Tokaanu, Misc). Coloured communities are from fields represented in the subpanel. Constraining variables are plotted as arrows (COND: conductivity, TURB: turbidity), with length and direction indicating scale and area of influence each variable had on the model. (b) The top panel represents a subset of the full CCA model, with select geothermal fields shown in colour (including 95 % confidence intervals). The bottom panel shows their respective geochemical signatures as a ratio of chloride (Cl-), sulfate (SO42-), and bicarbonate (HCO3-
The relationship between temperature and alpha diversity reported in this research starkly contrasts a previous intercontinental study comparing microbial community diversity in soil/sediments from 165 geothermal springs92, which showed a strong relationship (R2=0.40–
0.44) existed. In contrast, our data across the entire suite of samples, revealed temperature had no significant influence on observed community diversity (R2=0.002, P=0.201; Figure B.3, Table B.2). This result increased marginally for archaeal-only diversity (R2=0.013, P=0.0005), suggesting temperature has a more profound effect on this domain than it does Bacteria. However, the primers used in this study are known to be unfavourable towards some archaeal clades232, therefore it is likely extensive archaeal diversity remains undetected in this study. The lack of influence of temperature on whole community diversity was further substantiated via multiple linear modelling (Table B.3), and significance and correlation testing (Kruskal-Wallis: H=16.2, P=0.039; Pearson’s coefficient: |r|=0.04, P=0.201). When samples were split into pH increments, like Sharp et al. (2014)92, we observed increasing temperature only significantly constrained diversity above moderately acidic conditions (pH
>4; Figure B.5). However, the magnitude of this effect was, in general, far less than previously reported and is likely a consequence of the sample type (e.g., soil/sediments versus aqueous) and density of samples processed172. Many samples from geothermal environments are recalcitrant to traditional DNA extraction protocols and research in these areas has therefore focused on those with greater biomass abundance92,159 (i.e., soils, sediments, streamers or biomats). Whereas aqueous samples typically exhibit a more homogeneous chemistry and community structure, the heterogeneity of terrestrial samples is known to affect microbial populations (e.g., particle size, depth, nutrient composition)57. Our deliberate use of aqueous samples extends the results of previous small-scale work93,187 and also permits the robust identification of genuine taxa-geochemical relationships in these environments.
2.5.3 Microbial communities are influenced by pH, temperature and source fluid Throughout the TVZ, beta diversity correlated more strongly with pH (Mantel: ρ=0.54, P<0.001) than with temperature (Mantel: ρ=0.19, P<0.001; Figure 2.2, Table B.4). This trend was consistent in pH- and temperature-binned samples (Figure B.6; ANOSIM: |R|=0.46 and 0.18 respectively, P<0.001); further confirming pH, more so than temperature, accounted for observed variations in beta diversity. Congruent with our finding that pH influences alpha diversity at lower temperatures (<70 °C), the effect of temperature reducing beta diversity had greater significance above 80 °C (Wilcox: P<0.001; Figure B.6). The extent of measured
physicochemical properties across 925 individual habitats, however, allowed us to explore the environmental impact on community structures beyond just pH and temperature.
Permutational multivariate analysis of variance in spring community assemblages showed that pH (12.4 %) and temperature (3.9 %) had the greatest contribution towards beta diversity, followed by ORP (1.4 %), SO42- (0.8 %), turbidity (0.8 %), and As (0.7 %, P<0.001; Table B.5). Interestingly, constrained correspondence analysis of the 15 most significant, non-collinear, and variable parameters (Table B.5 & Table B.6; pH, temperature, turbidity, ORP, SO42-, NO3-, As, NH4+, HCO3-, H2S, conductivity, Li, Al, Si, and PO43-), along with geothermal field locations, only explained 10 % of variation in beta diversity (Figure 2.3), indicating physicochemistry, or at least the 46 parameters measured were not the sole drivers of community composition.
Geothermal fields are known to express chemical signatures characteristic of their respective source fluids233, implying autocorrelation could occur between location and geochemistry.
We therefore investigated whether typical geochemical conditions exist for springs within the same geothermal field and whether specific microbial community assemblages could be predicted. Springs are usually classified according to these source fluids; alkaline-chloride or acid-sulfate. High-chloride features are typically sourced from magmatic waters and have little interaction with groundwater aquifers. At depth, water-rock interactions can result in elevated bicarbonate concentrations and, consequently, neutral to alkaline pH in surface features. Acid-sulfate springs (pH <3), in contrast, form as steam-heated groundwater couples with the eventual oxidation of hydrogen sulfide into sulfate (and protons). Rarely, a combination of the two processes can occur; leading to intermediate pH values234. It is unknown, however, whether these source fluid characteristics are predictive of their associated microbial ecosystems. Bray-Curtis dissimilarities confirmed that, like alpha diversity (Kruskal-Wallis: H=240.7, P<0.001; Figure 2.4), community structures were significantly different between geothermal fields (ANOSIM: |R|=0.26, P<0.001; Figure B.7).
Gradient analysis comparing significant geochemical variables and geography further identified meaningful intra-geothermal field clustering of microbial communities (95 % CI;
Figure 2.3 & Figure B.8). Further, characteristic geochemical signatures from these fields were identified and analysis suggests they could be predictive of community composition.
For example, the Rotokawa and Waikite geothermal fields (approx. 35 km apart; Figure 2.3 position N & F) display opposing ratios of HCO3-, SO42-, and Cl-, with corresponding microbial communities for these sites clustering independently in ordination space. Despite
this association, intra-field variation in both alpha and beta diversity also occurred at other geothermal sites where geochemical signatures were not uniform across local springs (e.g., Rotorua, Figure 2.3 position D), demonstrating that correlation does not necessarily always occur between locational proximity and physicochemistry.
Table 2.1 - Average relative abundance and prevalence of phyla and genera. Only taxa above a 1 % average compositional threshold are shown. Maximum abundance of each taxon within individual features and standard deviation (SD) across all 925 springs. Where taxonomy assignment failed to classify to genus level, the closest assigned taxonomy is shown (f=family, o=order, p=phylum).
Phylum* Genus Abundance SD* Max Prevalence
Aquificae Venenivibrio 0.112 0.231 0.968 0.742
Proteobacteria Acidithiobacillus 0.111 0.242 0.994 0.629
Aquificae Hydrogenobaculum 0.100 0.235 0.999 0.608
Aquificae Aquifex 0.086 0.212 0.971 0.497
Deinococcus-Thermus Thermus 0.025 0.071 0.732 0.552
Proteobacteria Thiomonas 0.024 0.101 0.941 0.396
Proteobacteria Desulfurella 0.022 0.067 0.758 0.497
Crenarchaeota Sulfolobaceae (f) 0.020 0.091 0.951 0.416
Euryarchaeota Thermoplasmatales (o) 0.019 0.059 0.495 0.539
Proteobacteria Thiovirga 0.015 0.077 0.816 0.374
Proteobacteria Hydrogenophilaceae (f) 0.015 0.072 0.704 0.406
Thermodesulfobacteria Caldimicrobium 0.015 0.052 0.651 0.519
Proteobacteria Hydrogenophilus 0.013 0.045 0.432 0.484
Thermotogae Mesoaciditoga 0.011 0.033 0.286 0.410
Parcubacteria Parcubacteria (p) 0.010 0.024 0.193 0.608
*This table was generated before recent nomenclature changes to bacterial and archaeal phyla. For current valid names, please refer to Oren & Garrity, 2021198.
2.5.4 Aquificae and Proteobacteria are abundant and widespread
To determine whether individual microbial taxa favoured particular environmental conditions and locations, we first assessed the distribution of genera across all individual springs. Within 17 geothermal fields and 925 geothermal features, 21 phyla were detected with an average relative abundance >0.1 % (Figure 2.5). We found that two phyla and associated genera, Proteobacteria (Acidithiobacillus spp.) and Aquificae (Venenivibrio, Hydrogenobaculum, Aquifex spp.), dominated these ecosystems (65.2 % total average relative abundance across all springs), composing nine of the 15 most abundant genera >1 % average relative abundance (Table 2.1). Considering the broad spectrum of geothermal environmental conditions sampled in this study (we assessed microbial communities in springs across a pH
gradient of nine orders of magnitude and a temperature range of ~87 C) and the prevalence of these taxa across the region, this result was surprising. Proteobacteria was the most abundant phylum across all samples (34.2 % of total average relative abundance; Table 2.1), found predominantly at temperatures less than 50 °C (Figure B.9). Of the 19 most abundant proteobacterial genera (average relative abundance >0.1 %), the majority are characterised as aerobic chemolithoautotrophs, utilising either sulfur species and/or hydrogen for metabolism.
Accordingly, the most abundant (11.1 %) and prevalent (62.9 %) proteobacterial genus identified was Acidithiobacillus235, a mesophilic-moderately thermophilic, acidophilic autotroph that utilises reduced sulfur compounds, iron and/or hydrogen as energy for growth.
Aquificae (order Aquificales) was the second most abundant phylum overall (31 % average relative abundance across 925 springs) and included three of the four most abundant genera;
Venenivibrio, Hydrogenobaculum and Aquifex (11.2 %, 10.0 %, and 8.6 % respectively;
Table 2.1). As Aquificae are thermophilic (Topt 65–85 °C)236, they were much more abundant in warmer springs (>50 °C ; Figure B.9). The minimal growth temperature reported for characterised Aquificales species (Sulfurihydrogenibium subterraneum and S.
kristjanssonii)236 is 40 °C and may explain the low Aquificae abundance found in springs less than 50 °C. Terrestrial Aquificae are predominately microaerophilic chemolithoautotrophs that oxidise hydrogen or reduced sulfur compounds; heterotrophy is also observed in a few representatives236. Of the 14 currently described genera within the Aquificae, six genera were relatively abundant in our dataset (average relative abundance >0.1 %; Figure 2.5): Aquifex, Hydrogenobacter, Hydrogenobaculum, and Thermocrinis (family Aquificaceae); and Sulfurihydrogenibium and Venenivibrio (family Hydrogenothermaceae). No signatures of the Desulfurobacteriaceae were detected and is consistent with reports that all current representatives from this family are associated with deep-sea or coastal thermal vents236. Venenivibrio (OTUs; n=111) was also the most prevalent and abundant genus across all communities (Table 2.1). This taxon, found in 74.2 % (n=686) of individual springs sampled, has only one cultured representative, V. stagnispumantis (CP.B2T), which was isolated from the Waiotapu geothermal field in the TVZ179. The broad distribution of this genus across such a large number of habitats was surprising, as growth of the type strain is only supported by a narrow set of conditions (pH 4.8–5.8, 45–75 °C). Considering this, and the number of Venenivibrio OTUs detected, we interpret this result as evidence there is substantial undiscovered phylogenetic and physiological diversity within the genus. The ubiquity of Venenivibrio suggests that either the metabolic capabilities of this genus extend substantially
beyond those described for the type strain, and/or that many of the divergent taxa could be persisting and not growing under conditions detected in this study209,237.
2.5.5 Geochemical and geographical associations exist at the genus level
The two most abundant phyla, Proteobacteria and Aquificae, were found to occupy a characteristic ecological niche (<50 °C and >50 °C, respectively; Figure B.9). To investigate specific taxa-geochemical associations beyond just temperature and pH, we applied a multivariate linear model to determine enrichment of taxa in association with geothermal fields and other environmental data (Figure 2.5). The strongest associations between taxa and chemistry (Z-score >4) were between Nitrospira – nitrate (NO3-) and Nitratiruptor – phosphate (PO43-). Nitrospira oxidises nitrite to nitrate and therefore differential high abundance of this taxon in nitrate-rich environments is expected. Further, the positive Nitratiruptor-PO43- relationship suggests phosphate is a preferred nutritional requirement for this chemolithoautotroph238 and informs future efforts to isolate members of this genus would benefit from additional phosphate or the presence of reduced phosphorous compounds in the culture medium239,240. Thermus and Hydrogenobaculum were the only bacterial taxa to differentially associate (compared to other taxa) positively and negatively with pH respectively. This is consistent with the lack of acidophily phenotype (pH<4) reported in Thermus spp.241 and conversely, the preferred acidic ecological niche of Hydrogenobaculum242. Aquifex was the only genus to display above average association with temperature, confirming abundance of this taxon is significantly enhanced by hyperthermophily243.
Similar to the chemical-taxa associations discussed above, differential abundance relationships were calculated with respect to individual geothermal fields (Figure 2.5). A geothermal field, which contains springs across the pH scale (i.e., Rotorua), was closely associated with the highly abundant and prevalent Acidithiobacillus and Venenivibrio. On the other hand, a predominantly acidic geothermal system (i.e., Te Kopia), produced the only positive associations with ‘Methylacidiphilum’ (Verrucomicrobia), Acidimicrobium (Actinobacteria), Terrimonas (Bacteroidetes) and Halothiobacillus (Proteobacteria). These relationships are likely describing sub-community requirements that are otherwise not captured by conventional spatial-statistical analysis, therefore providing insight into previously unrecognised microbe-niche interactions.
Figure 2.4 - Alpha and beta diversity as a function of geographic distance. (a) Alpha diversity scales (via Shannon index) across individual springs (n=925), split by geothermal fields. Fields are ordered from north to south (H: Kruskal-Wallis test). (b) A distance-decay pattern of beta diversity (via Bray- Curtis dissimilarities of 925 springs) against pairwise geographic distance in metres, with linear regression applied (m=slope). Geographic distance is split into bins to aid visualisation of the spread.
Figure 2.5 - Taxonomic association with location and physicochemistry. The heat map displays positive (red) and negative (blue) association of genus-level taxa (>0.1 % average relative abundance) with each geothermal field and significant environmental variables, based on Z-scores of abundance log ratios. Each taxon is colour-coded to corresponding phylum* on the approximately maximum- likelihood phylogenetic tree.
*This figure was generated before recent nomenclature changes to bacterial and archaeal phyla. For current valid names, please refer to Oren & Garrity, 2021198.
2.5.6 Distance-decay patterns differ at local and regional scales
Environmental selection, ecological drift, diversification, and dispersal limitation all contribute to distance-decay patterns62. While several studies have shown microbial dispersal limitations and distance-decay patterns exist in diverse geothermal and non-geothermal environments (e.g., hot springs6, freshwater streams199, salt marshes244), the point of inflection between dispersal limitation and selection, at regional and local geographic scales, remains under-studied. We identified a positive distance-decay trend with increasing geographic distance between 925 geothermal spring communities across the TVZ region (slope=0.031, P<0.001; Figure 2.4). This finding strongly suggests dispersal limitation exists between individual geothermal fields. Increasing the resolution to within individual fields, distance-decay patterns are negligible compared to the regional scale (Table B.7).
Interestingly, the greatest pairwise difference (y=1) between Bray-Curtis dissimilarities was also observed in springs classified as geographically-adjacent (<1.4 m). In the 293 geothermal spring pairs separated by <1.4 m, temperature had a greater correlation with beta diversity than pH (Spearman’s coefficient: ρ=0.44 and 0.30 respectively, P<0.001). This result illustrates the stark spatial heterogeneity and selective processes that can exist within individual geothermal fields. Congruently, each OTU was detected in an average of only 13 springs (Figure B.4). We propose that physical dispersal within geothermal fields is therefore not limiting, but the physicochemical diversity of geothermal springs acts as a barrier to the colonisation of immigrating taxa. However, even between some neighbouring springs with similar (95% CI) geochemical signatures, we did note some dissimilar communities were observed (for example, Waimangu geothermal field; Figure 2.3 position E). These differing observations can be explained either one of three ways: Firstly, the defining parameter driving community structure was not one of the 46 physicochemical variables measured in this study (e.g., dissolved organic carbon or hydrogen); secondly, through the process of dispersal, the differential viability of some extremophilic taxa restricts gene flow and contributes to population genetic drift within geothermal fields86. We often consider
‘extremophilic’ microorganisms living in these geothermal environments as the epitome of hardy and robust. In doing so, we overlook that their proximal surroundings (i.e. immediately outside the host spring) may not be conducive to growth and survival245 and therefore the divergence of populations in neighbouring, chemically-homogeneous spring ecosystems is plausible. Thirdly, aeolian immigration107 from the non-geothermal surrounding environment could alter the perceived composition of a community, even when immigrants are not competing to survive. Future work could be undertaken to understand individual population
responses to community-wide selective pressures and the temporal nature of ecosystem functioning.