Algae-Bacteria Community Analysis for Drinking Water Taste and Odour Risk Management

1. Introduction

The predominant source of customer complaints to drinking water companies worldwide is the earthy-musty flavours imparted by trans-1,10-dimethyl-trans-9-decalol (geosmin) and 2-methylisoborneol (2-MIB) (Devi et al., 2021; Pochiraju et al., 2021; Wu et al., 2022). These taste and odour (T&O) compounds exhibit low odour thresholds, with geosmin and 2-MIB ‘events’ defined as water concentrations of ≥10 ng L^-1 and ≥5 ng L^-1, respectively (Hooper et al., 2023). Due to these low concentration thresholds and the high cost of treatment processes such as activated carbon and ozonation (Mustapha et al., 2021), water companies aim to limit the production of these compounds in drinking water reservoirs.

Cyanobacterial blooms (i.e., large increases in biomass) have been considered the leading cause of T&O production (Adams et al., 2021; Bruder et al., 2014; Hayes & Burch, 1989; Newton et al., 2015). However, biomass alone does not explain all potential variations in T&O compound concentrations: significant relationships exist between periods of high cyanobacterial productivity and T&O metabolite production, irrespective of total biomass (Zimmerman et al., 1995). Elevated cyanobacterial productivity could explain how some studies have been unsuccessful in correlating T&O compounds with cyanobacterial biomass (Graham et al., 2010; Watson et al., 2008). For example, Kim et al. (2018) found that concentrations of geosmin per unit biomass were maximal when biomass was lowest, but productivity was high. While we know certain environmental conditions influence T&O metabolite production (i.e, high ammonia and low nitrate; (Hooper et al., 2023; Perkins et al., 2019), it remains unclear how biological factors at the community level affect cyanobacterial productivity-driven T&O events.

Cyanobacterial productivity is influenced by a variety of factors. A substantial amount of research has focussed on the abiotic drivers (e.g., nutrients, temperature, light; (Harris et al., 2016; Montgomery, 2014; Robarts & Zohary, 1987), however, recent work has shown that the wider bacterial communities associated with Cyanobacteria are very influential (Louati et al., 2015; Salomon et al., 2003). Heterotrophic bacteria can change the physiochemical properties of the water environments (Havens, 2008) creating favourable/ unfavourable environmental conditions for cyanobacterial productivity-driven T&O production. Previous studies revealed that cyanobacterial bloom-associated microbial communities are highly varied at the phylum level (Louati et al., 2015; Steffen et al., 2012). However, gene function could be highly conserved despite differences in taxonomic assemblages. It may therefore be possible to define cyanobacterial T&O risk categories based on indicator taxa assigned ecological traits such as nitrogen fixation and assimilation (Louati et al., 2015), organic matter degradation (Berg et al., 2009), and iron uptake (Ploug et al., 2011). Former research has focussed on the indirect and direct effects of bacterial communities on cyanobacterial T&O metabolite production (Berg et al., 2009; Graham et al., 2010; Youn et al., 2020). In doing so, a range of interactions have been observed, all revealing links between certain bacterial communities and increases in the cyanobacterial production of T&O metabolites, including geosmin and 2-MIB. Bacterial communities, however, can also degrade the high molecular weight organic compounds formed by Cyanobacteria (Cai et al., 2014), constituting key roles in nutrient recycling. Bacterial communities associated with Cyanobacteria can also lower oxygen tension near cyanobacterial cells where oxygen-sensitive biochemical processes (photosynthesis, N-fixation) occur (Paerl & Kellar, 1978). These interactions may enhance cyanobacterial productivity and, in turn, the production of T&O compounds. Additionally, these bacterial communities participate in other processes that have the potential to promote or inhibit T&O metabolite production by supporting cyanobacterial growth or altering environmental conditions favourable to T&O metabolite release (Eiler & Bertilsson, 2004; Salomon et al., 2003; Simon et al., 2002).

Relationships also exist between Cyanobacteria and algae, some of which have the potential to affect productivity. Cyanobacteria provide fixed nitrogen (N) to diatoms in return for amino acids and organic carbon (Foster et al., 2011; Hilton et al., 2013; Thompson et al., 2012). They also help release phosphorus (P) by triggering green algae (Chlorophyta) to secrete alkaline phosphatase (APase), which converts organic P into usable phosphate (Bar-Yosef et al., 2010). Additionally, algal communities help detoxify extracellular iron (Fe) compounds produced by Cyanobacteria, which could otherwise be toxic to the Cyanobacteria in low-oxygen waters (Paerl & Pinckney, 1996). Shared functional traits of algal communities associated with Cyanobacteria could thus accelerate nutrient cycling, facilitate cyanobacterial productivity, and potentially cause a T&O event risk. Although there has been considerable research on communities associated with Cyanobacteria, no studies so far have explored interactions within the broader planktonic communities, including bacteria and algae, or their role in T&O events. This is a crucial gap in our knowledge, as these interactions may affect cyanobacterial T&O metabolite production.

Identifying relationships within communities associated with Cyanobacteria is critical for understanding the mechanisms underlying T&O metabolite production and, ultimately, moving towards predicting T&O events. Advances in molecular methods, particularly correlation-based network analyses, have been successful in exploring microbial community co-occurrence and co-exclusion patterns (Barberán et al., 2012; Eiler et al., 2012; Ju & Zhang, 2015). These molecular techniques, which use the relative abundances of marker genes (e.g., 16S rRNA and rbcL), allow the identification of organisms (and thus potential inter-specific interactions) that were previously impossible to detect, particularly within cyanobacterial and algal communities. Such molecular methods provide novel insights into the complex dynamics of these communities and their roles in T&O events (Faust & Raes, 2012). The co-occurrence of organisms reflects potential positive interactions (e.g., facilitation and mutualisms), and negative interactions (e.g., co-exclusion due to competition of resources or differing environmental preferences) or indicates taxa that share similar environmental niches (Faust, 2021). The knowledge of both community dynamics and potential interactions afforded by these methods can offer a potential to explore and understand the relationships between communities associated with Cyanobacteria. Investigating such relationships can highlight taxa potentially implicated in T&O production and possibly pinpoint indicative taxa related to increases in geosmin and 2-MIB concentrations.

In this study, next-generation sequencing of 16S rRNA and rbcL genes was employed to assess bacterial and algal community dynamics across three drinking water reservoirs in Wales, UK to explore ecological relationships and identify biotic drivers of T&O events. One reservoir exhibited elevated geosmin concentrations, another showed elevated 2-MIB levels, and the third had neither geosmin nor 2-MIB concentrations exceeding the event threshold, as defined by Hooper et al. (2023) (geosmin > 10 ng L⁻¹ and 2-MIB > 5 ng L⁻¹). The research aimed to 1) identify bacterial and algal communities associated with elevated T&O metabolite concentrations; 2) determine whether there are taxa indicative of low and high T&O concentrations; and 3) investigate the co-occurrence networks surrounding T&O-producing Cyanobacteria and understand relationships with environmental factors. The aim was to identify taxa associated with T&O event risk, offering valuable insights for water companies and environmental managers.

2. Materials and Methods

2.1. Reservoir Locations and Sample Collection

Between July 2019 and August 2020, 500 mL of reservoir water (maximum depth 0.5 m) was collected monthly using bankside sampling for molecular analysis from multiple sites in three different reservoirs across Wales, UK (map available in Online Resource 1). Additional samples were taken simultaneously from the same locations using the same sampling method and used for chemical analysis. Utilising multiple sites in each reservoir reflects a holistic understanding of ecosystem dynamics, capturing the spatial heterogeneity of the reservoirs. This approach allowed for a comprehensive assessment of T&O risk by collating community compositions across the reservoirs. Moreover, it aligns with water quality management practices, which typically examine the reservoir as a single system. The exception to this sampling regime was a sampling break between March and April 2020 due to COVID-19 restrictions. Samples were transported by a refrigerated vehicle and kept at 4 °C before laboratory analysis. During the sampling period, Reservoir 1 (area covered = 360 ha) experienced geosmin events (≤ 520 ng L^-1, n = 32, sites = 5), Reservoir 2 (area covered = 102 ha, n = 30, sites = 2) experienced 2-MIB events (≤ 60 ng L^-1) and Reservoir 3 (area covered = 176 ha, n = 25, sites = 4) did not experience elevations in either T&O compound (≤1.80 ng L^-1 for both geosmin and 2-MIB). A summary of geosmin and 2-MIB concentrations for each studied reservoir is provided in Online Resource 2.

2.2. Geosmin and 2-MIB Concentrations

Gas Chromatography-Tandem Mass Spectrometry (GC–MS/MS) (Agilent, Santa Clara, USA) was used to determine concentrations of geosmin and 2-MIB. Water samples were extracted using a solid phase extraction (SPE) technique at a Dŵr Cymru Welsh Water accredited laboratory (ISO/IEC 17025:2017). Full details of the methodology are listed in (Hooper et al., 2023).

Categorising Geosmin and 2-MIB concentrations

Low, medium, and high geosmin and 2-MIB concentration categories were proposed to aid data analysis and support the prediction of T&O risk. Drawing distinctions for T&O risk has previously been challenging without clearly defined categories, which are essential for effective modelling and forecasting. By establishing these categories, we can better interpret complex data sets and identify trends or patterns that may indicate forthcoming T&O events. Categories were assigned based on the event-level thresholds set out by Hooper et al. (2023) (geosmin > 10 ng L^-1 and 2-MIB >5 ng L^-1) and are presented in Table 1.

2.3. Water Chemistry

Water chemistry measurements were conducted at Dŵr Cymru Welsh Water accredited laboratory in line with ISO/IEC 17025:2017. Ammonium (NH₄⁺), nitrate (NO₃^-), nitrite (NO₂^-) and total phosphorus (TP) were analysed using an Aquakem 600 discrete analyser (Thermo Scientific, Waltham, USA) with methods described in (Hooper et al., 2023). pH was measured via a pH meter using the Skalar SP2000 robot (Skalar Analytical B. V., Breda, Netherlands). For manganese (Mn) and dissolved iron (DFe), water samples were digested in 1% hydrochloric acid at 85 ± 5°C overnight before being introduced to the Agilent 7700 Inductively Coupled Plasma Mass Spectrometer (ICP-MS) with a Cetac ASX-500 Autosampler (Agilent Technologies, Santa Clara, USA).

2.4. eDNA Extraction and Next-Generation Sequencing (NGS)

2.4.1. Extracting eDNA from Reservoir Water

Water samples (500 ml) were filtered through a Sterivex filter (0.2 µM) using a vacuum manifold. 1 mL of ATL buffer (Qiagen, Germany) was added to the Sterivex filter before storing each sample at −20 °C for later use. Before extraction, 100 µL of proteinase K was added to each filter and left on a turntable for 2 hours. DNA was extracted using ~100 µL of the sample removed from each Sterivex filter, following methods described by Fawley and Fawley (2004). Briefly, samples (~100 μL) were mixed with 200 μL of extraction buffer (5M NaCl, 30 mM NaEDTA and 70 mM tris pH 7.0), 25 μL of 10% DTAB and 200 μL of chloroform and added to a tube containing 0.1 mm glass beads. Samples were then agitated using an MP Biomedical, FastPrep-24™ at 5 ms^-1 for 30 seconds twice (with a 5-minute incubation at room temperature), before being centrifuged (Eppendorf, 5417C) at 14000 rpm for 2 minutes. The top phase of samples was transferred to a sterile collection tube (200 μL), and 200 μL of buffer AL (QIAGEN, Germany) added along with 200 μL of ethanol before being vortexed for 30 seconds. DNA was then purified using the procedure and reagents provided within the DNAeasy^® Blood & Tissue Kit (QIAGEN, Germany).

2.4.2. Next Generation Sequencing (NGS) – 16S rRNA and rbcL genes

Two gene markers were selected to target algae (rbcL) and bacteria (16s rRNA). Initial 16S rRNA and rbcL PCRs were carried out on the extracted eDNA samples using primers listed in Table 2 with incorporated Nextera tags on the 5’ end of the primers. Compositions of 20 µL PCR reactions were consistent for both target genes, containing 5 µL of 4x AllTaq Mastermix (QIAGEN, Hilden, Germany), 0.30 pmol µL^-1 of each primer mentioned above and 1 µL of template DNA. PCR reactions were carried out on a SimpliAmp™ Thermal Cycler (Thermo Fisher Scientific, Waltham, USA), using the following conditions: initial denaturation at 95 °C for 2 minutes followed by 40 cycles of 95°C for 5 seconds (denaturation), 55°C for 15 seconds (annealing), and 72°C for 10 seconds (extension) with a single final extension at 72°C for 5 minutes. All reactions were performed in triplicate to account for PCR bias. PCR success was determined using a QIAxcel (QIAGEN, Hilden, Germany). The successful triplicates were subsequently pooled together and underwent a cleaning stage using a Zymo Research 96 DNA clean-up kit (Zymo Research, Cambridge, UK) according to the manufacturer’s instructions.

Cleaned amplicon samples underwent a secondary amplification to anneal Illumina-Nextera indexes (Illumina, San Diego, USA). 12.5 µL KAPA HiFi HotStart ReadyMix, 5 µL of nuclease-free water and 2.5 µL of each unique Nextera index combination along with 2.5ul of cleaned primary product. This underwent PCR amplification on a MultiGene™ OptiMax Thermal Cycler (Labnet International, Edison, USA) using the following conditions: initial denaturation at 95 °C for 3 minutes followed by eight cycles of 95°C for 30 seconds (denaturation), 55°C for 30 seconds (annealing), and 72°C for 30 seconds (extension) with a single final extension at 72°C for 5 minutes.

SequelPrep normalisation plates (Invitrogen, Waltham, USA) were used to normalise the concentrations of the 16S rRNA and rbcL amplicons as per the manufacturer’s instructions. Before sequencing, normalised samples were pooled together. Samples were sequenced for both genes (16S rRNA and rbcL) on separate runs with either an Illumina MiSeq-nano or v2 cartridge using 2 x 250 bp paired-end reads.

2.4.3. Bioinformatic analysis and databases for taxonomic assignment

Sequence data were processed using QIIME2 v2021.2 (Bolyen et al., 2019). Fastq files from replicate samples were concatenated, and paired-end sequences were demultiplexed and imported into QIIME2. Sequences were denoised, merged, and chimera-checked using the DADA2 pipeline (Callahan et al., 2016). Run-specific truncation was applied based on read quality scores, and samples from separate runs were merged via amplicon sequence variant (ASV) tables.

Taxonomic assignments were made using two databases: (1) 16S rRNA V4 region (515-806 bp) extracted from the SILVA v132 database (Callahan, 2018); and (2) a custom rbcL database, compiled by merging sequences from Diat.barcode v9.2 (accessed 25 October 2021) (Rimet et al., 2019) and GenBank (with overlapping entries removed). GenBank sequences were cross-referenced with NCBI taxdump files for taxonomic classification. Both databases were formatted for QIIME2 and processed using the extract-reads function with a 0.8 identity threshold to filter for primer-containing sequences.

ASV assignment was performed using a Naive Bayes classifier in QIIME2. Post-assignment, features were filtered to retain those with a minimum relative abundance of 0.1% in at least 0.1% of samples.

2.5. Statistical Analyses

2.5.1. Datasets – Proportional and Binomial Manipulation

All statistical analyses were completed in R 4.1.0 (R_Core_Team, 2021). Sequence data for 16S rRNA and rbcL were exported at a genus level (or the closest available lineage) and converted into proportions separately by calculating the percentage of total reads for each genus based on the total number of reads per sample.

Joint analyses of algae and bacterial communities required the 16S rRNA and rbcL proportional datasets to be merged. To do so, we transformed datasets into binomial data to prevent biases towards either group resulting from differences in read depth or amplification efficiencies between bacterial and algal sequences. This transformation allows for a balanced representation of both groups in the analysis, mitigating artefacts that might skew community interactions. Merging different data types in ecological studies is a common approach to reveal broader community dynamics, and similar methods have been applied in metabarcoding studies to combine datasets while addressing such biases (Cuff et al., 2022) Presence (1) was recorded as a taxon proportion ≥1% to reduce noise from rare taxa and to focus on ecologically relevant taxa and absence (0) <1%.

Environmental variables (pH, NH₄⁺, NO₃^- –, TP – total phosphorus, DMn – dissolved manganese and DFe – dissolved iron), were converted into the categories based on the quartiles of the distribution range of each variable (Q1 ≤ 0.25, Q2 0.25-0.50, Q3 0.50-0.75 and Q4 > 0.75). Samples with values for each variable that fell into a specific quartile range were marked present (1) and absent (0) for the remaining quartile ranges. This binomial categorisation facilitated analyses alongside the binomial community dataset, enabling comparisons and interpretations of community dynamics in relation to environmental conditions.

2.5.2. Non-metric Multidimensional Scaling of 16S rRNA and rbcL communities

Non-metric Multi-Dimensional Scaling (NMDS) was used to assess changes in the community structure of 16S rRNA and rbcL communities (using proportional data) for reservoirs experiencing geosmin (Reservoir 1) and 2-MIB (Reservoir 2) events. NMDS was calculated using the ‘vegan’ package (Oksanen, 2013) using a Bray-Curtis distance, with Wisconsin double standardisation to account for the effects of common and rare taxa. Results were plotted using ‘ggplot2’ (Wickham, 2016).

2.5.3. Random Forest (RF) Models

Random Forest^TM (RF) models were used to identify taxa that could predict low and high T&O risk by classifying samples into low, medium, and high geosmin (for Reservoir 1) and 2-MIB (for Reservoir 2) concentration categories. The purpose of this analysis was to determine which bacterial and algal taxa were most indicative of elevated T&O levels. Using the R package ‘randomForest’ (Breiman, 2001), taxa were treated as independent variables, with the T&O concentration categories as the dependent variable.

The number of trees was first optimised for both reservoirs by running preliminary RFs with 10⁵ trees; the point at which the error rate plateaued was selected for subsequent models (10⁴ trees for Reservoir 1 and 2 x 10³ trees for Reservoir 2). The number of variables used in each tree was optimised by cross-validation using 100 folds. After refining the number of trees and variables (Reservoir 1 – n = 197 and Reservoir 2 – n =144), a final model was run for each reservoir. To assess the models' performance, multi-dimensional scaling plots were created from the proximity matrices for the geosmin and 2-MIB RFs.

To identify indicative organisms relating to T&O levels, taxa in the RF models were ranked in descending order of the Mean Decrease Accuracy (MDA) (the decrease in model accuracy from permutations of each taxon), and the top 20 taxa extracted. The top 20 taxa were used to reduce model complexity, improve model performance and highlight key indicator taxa. RFs were performed for this subset of taxa for Reservoir 1 (geosmin) and Reservoir 2 (2-MIB), optimising the number of trees (Reservoir 1 – 2 x 10³ trees and Reservoir 2 – 200 trees) and the number of variables used (Reservoir 1 – 20 taxa and Reservoir 2 – 5 taxa) as described above. Multi-dimensional scaling plots of the proximity matrices from the RFs were again constructed. MDA values associated with the classification of low and high T&O concentrations were further investigated to detect taxa that can differentiate between the two categories (Online Resource 9, geosmin and Online Resource 10, 2-MIB); total mean abundances for each taxon were also calculated to determine average abundances during low and high T&O concentrations.

2.5.4. Co-occurrence Analyses

To investigate the potential linkages between taxa in the algal and bacterial communities, as well as associations with specific environmental conditions, co-occurrence analyses were performed separately for data from each Reservoir using binomial data and the ‘cooccur’ package (Veech, 2013) (parameters listed in Online Resource 3). This followed the approach outlined in Faust (2021). Environmental categories (see Section 2.5.1) were included alongside community data to account for the influence of environmental conditions (i.e., a common response to environmental factors) as well as the direct effects of taxa interactions or connectivity (Faust, 2021). The probabilistic co-occurrence model calculated P-values for pairwise associations (Griffith et al., 2016). Associations are characterised by the probability that the co-occurrence between two taxa is either greater or less than expected by chance (i.e., in randomly assembled communities). Specifically, some taxa pairs co-occur more frequently than expected, indicating a higher-than-expected co-occurrence, while others co-occur less frequently, indicating a lower-than-expected co-occurrence. Significant pairwise associations between taxa and environmental categories (α = 0.01) were extracted and represented as association networks. The nodes in the networks represent taxa and environmental categories and edges reflect the probability of the pairwise associations. Datasets with all nodes and edges used to construct the networks are available in Online Resource 4.

Subnetworks of T&O-producing Cyanobacteria and their associated communities were created to visualise these higher and lower-than-expected co-occurrences. These networks were visualised using ‘igraph’ (Csardi & Nepusz, 2006), ‘ggnetwork’ (Briatte, 2023), and ‘ggplot2’ (Wickham, 2016). Interactive visualisations of the full networks, created using the ‘visNetwork’ package (Almende et al., 2019), are available in Online Resource 5.

2.6. Data Availability

16S rRNA and rbcL Next Generation Sequencing data from this study were submitted to the NCBI Sequence Read Archive (SRA) under the accession number PRJNA1080335. All additional data and code are available on Zenodo (10.5281/zenodo.10671684).

3. Results

3.1. T&O Categories in NMDS Community Analysis

16S rRNA genera composition in Reservoir 1 (geosmin) revealed distinct communities associated with geosmin levels (Figure 1a). Communities sampled during periods of low and medium geosmin levels did not cluster as distinctively as high geosmin levels and overlapped. Reservoir 2 (2-MIB) 16S rRNA communities also exhibited distinct clustering within the 2-MIB concentration categories (Figure 1b). The community composition for rbcL in Reservoir 1 did not appear to cluster according to geosmin levels (Figure 1c). However, In Reservoir 2 (2-MIB), rbcL genera composition can be seen to form a cluster for data points subjected to high 2-MIB levels (Figure 1d).

Random Forest (RF) analysis, using the top 20 taxa as independent variables, correctly classified 100% of samples for high geosmin (>20.00 ng L^-1) and 80% of samples for low geosmin concentrations (<5 ng L^-1) in Reservoir 1. For 2-MIB concentrations in Reservoir 2, RF analysis correctly classified 100% of samples for high (>10 ng L^-1) and low concentrations (<2.50 ng L^-1) (confusion matrices provided in Online Resource 6 and multi-dimensional scaling plots of RF proximity matrices provided in Online Resource 7). For both T&O compounds, the classification of samples for medium concentrations was less accurate (geosmin: 5.00 – 20.00 ng L^-1 and 2-MIB: 2.50 – 10.00 ng L^-1). RF classification of samples for medium geosmin concentrations had a 50% error rate for Reservoir 1, and medium concentrations of 2-MIB in Reservoir 2 had a 100% error rate.

T&O Risk Categorisation of Indicative Taxa Identified by RF Analysis

The importance of each taxon from both models was assessed by exploring overall mean decrease accuracy (MDA) values (Online Resource 8). MDA values associated with the classification of low and high T&O concentrations were further investigated to detect taxa that can differentiate between the two categories (Online Resource 9, geosmin and Online Resource 10, 2-MIB); total mean abundances for each taxon were also calculated to determine average abundances during low and high T&O concentrations.

Three Cyanobacteria (Microcystis PCC-7914, Aphanizomenon NIES81 and Cyanobium PCC-6307) were influential in discriminating low and high geosmin concentrations (Reservoir 1), although of these only Aphanizomenon NIES81 was itself capable of geosmin production (Driscoll et al., 2018). Interestingly, Cyanobium PCC-6307 is a known 2-MIB producer (Jakubowska & Szeląg-Wasielewska, 2015) and Microcystis PCC-7914 is a known toxin producer, producing microcystins (Graham et al., 2010). Aphanizomenon NIES81 was more important in classifying low than high geosmin concentrations; its mean relative abundances show that small proportions of Aphanizomenon NIES81 (0.12 ± 0.4%) are characteristic of low geosmin concentrations, and higher proportions (14.96 ± 19.82%) are typical of high geosmin concentrations. No relative abundance of Microcystis PCC-7914 was recorded during low geosmin concentrations, but this genus received the third-highest MDA value for high geosmin concentrations.

In Reservoir 2, only one 2-MIB-producing Cyanobacteria (Cyanobium PCC-6307) was important in classifying low and high 2-MIB concentrations. Three taxa in this model had overall MDA values of 0 (Legionella, Gammaproteobacteria Class and the Comamonadaceae Family) and were disregarded from interpretation.

3.3. Co-occurrence of Bacterial (16S rRNA) and Algal (rbcL) Communities

In Reservoir 1 (geosmin), 311 pairwise comparisons yielded statistically significant associations between taxa and environmental variables; 87 pairwise associations exhibited lower-than-expected co-occurrence, and 224 pairwise associations exhibited higher-than-expected co-occurrences (Figure 2a). In this network, there were four significantly co-occurring Cyanobacteria taxa, only two of which are known geosmin producers: Aphanizomenon NIES81 and Nostocaceae (Driscoll et al., 2018). Aphanizomenon NIES81 (Figure 2b) showed higher-than-expected associations with Fluviicola, the 2-MIB producer Cyanobium PCC-6307 (Jakubowska & Szeląg-Wasielewska, 2015), and low concentrations of NO₃⁻ and NO₂⁻ and high geosmin concentrations (Figure 2c). Aphanizomenon NIES81 also had 13 lower-than-expected associations, including Bacillus, a known geosmin-degrading bacterium (Zhi et al., 2016). Additionally, Aphanizomenon NIES81 was associated with lower-than-expected co-occurrence in low pH, high NO₃⁻/NO₂⁻ and low geosmin concentrations. The geosmin-producing Nostocaceae showed one higher-than-expected association with the Chlamydomonadaceae Family and no significant lower-than-expected co-occurrence associations.

In Reservoir 2 (elevated 2-MIB site), 179 pairwise comparisons yielded statistically significant co-occurrences: 73 pairwise associations showed lower-than-expected co-occurrence and 106 pairwise associations showed higher-than-expected co-occurrence (Figure 2d). The 2-MIB-producing Dolichospermum NIES41 (Zhang et al., 2017) had one higher-than-expected association with the diatom Melosira. The 2-MIB producer, Cyanobium PCC-6307 (Jakubowska & Szeląg-Wasielewska, 2015) exhibited lower-than-expected co-occurrence with five taxa: Rubinisphaeraceae Family, R7C24, Gomphonema, Ulnaria, and Gemmataceae Family, and higher-than-expected co-occurrence with Asterionella, Sporichyaceae Family, and Rhizobiales A0839 (Figure 2e). Cyanobium PCC-6307 also showed lower-than-expected associations with low pH and high dissolved manganese (DMn), and one higher-than-expected co-occurrence with low DMn concentrations. Interestingly, high concentrations of 2-MIB were not associated with Cyanobacteria or any taxon capable of T&O production at a higher or lower-than-expected co-occurrence (Figure 2f).

Additional details on taxa and environmental variables showing higher-than-expected or lower-than-expected co-occurrence with T&O-producing Cyanobacteria in Reservoirs 1 (geosmin) and 2 (2-MIB) are provided in Online Resource 11 and 12.

In Reservoir 3 (no concentrations of geosmin or 2-MIB above event level threshold), 224 pairwise comparisons revealed statistically significant co-occurrences; 66 pairwise associations were lower-than-expected, while 158 pairwise associations were higher-than-expected (Figure 2g). In Reservoir 3, only two Cyanobacteria displayed significant co-occurrences: the 2-MIB-producing Cyanobium PCC-6307 and Aphanizomenon MDT14a (not a known geosmin producer according to genome analysis by Driscoll (2016)) (Figure 2h). Interestingly, the Cyanobacteria present in Reservoir 3’s network exhibited more higher-than-expected pairwise co-occurrences compared to the T&O-producing Cyanobacteria in Reservoirs 1 and 2, which had a higher frequency of lower-than-expected associations. For Cyanobium PCC-6307, higher-than-expected co-occurrences existed with Aphanizomenon MDT14a, Sporichthyaceae Family, Limnohabitans, and Fragilaria, as well as high pH and low TP and DFe concentrations. This community network had one lower-than-expected co-occurrence with Schlesneria and low pH and NO₂^- concentrations. The potential importance of these taxa and environmental variables associated with the 2-MIB-producing Cyanobium PCC-6307 is summarised in Online Resource 13.

3.4. T&O Risk Indicator Categories and Associated Taxa

Literature searches of taxa identified by RF models for Reservoir 1 (geosmin) and Reservoir 2 (2-MIB) (Subsection 3.2) were conducted and ecological functions were examined. T&O risk categories were then established based on the shared ecological functions of key taxa identified in the RF models (Table 3). These categories were further applied to taxa significantly co-occurring with T&O-producing Cyanobacteria, organising taxa into the same categories (Table 3).

4. Discussion

This study is the first to explore T&O risk associated with Cyanobacteria by integrating bacterial and algal communities in combination with environmental datasets to identify indicator taxa organised into T&O-risk categories. Identifying bacterial and algal signature communities of low, medium, and high T&O concentration categories indicates that T&O risk can be inferred from these communities. In addition, RF analysis revealed T&O-producing Cyanobacteria are influential in distinguishing between low and high T&O concentrations – critical threshold extremes for water companies in managing water quality. Community network analysis was consistent with findings from the RF analysis, identifying T&O-producing Cyanobacteria from Reservoir 1 (geosmin) and Reservoir 2 (2-MIB) to be significantly associated with high T&O concentration categories, thereby identifying the potential source of T&O metabolite production. T&O-producing Cyanobacteria in Reservoir 1 (geosmin) and Reservoir 2 (2-MIB) experienced more lower-than-expected co-occurrences with other taxa, suggesting that they are less frequently associated with other taxa in a stable or cooperative manner; this could indicate a more selective or unstable ecological role of these T&O-producing Cyanobacteria in the community. The unstable community structure associated with T&O-producing Cyanobacteria may influence their physiology, with T&O metabolite production potentially linked to stress-induced conditions. However, T&O production as a stress response can only occur during the right environmental conditions (low NO₃^- and NO₂^-). In contrast, significantly co-occurring T&O-producing Cyanobacteria with other taxa present in Reservoir 3 (no T&O events) exhibited a greater frequency of higher-than-expected co-occurrences with other taxa. T&O-producing Cyanobacteria in Reservoir 3 may display a more stable or integrated role, indicating stronger or more predictable associations with other taxa within Reservoir 3’s ecosystem. In identifying taxa from RF and community network analysis, we propose T&O-risk categories based on the shared function-specific traits of the indicator taxa (T&O-risk categories are discussed below).

4.1. Categorisation of indicative taxa for T&O risk

Ordination analysis revealed that bacterial (16S rRNA) and algal (rbcL) communities clustered by T&O concentration categories (low, medium and high) across Reservoir 1 (geosmin) and Reservoir 2 (2-MIB), indicating that the biology consortium as a response to water quality can reflect T&O concentrations. Both bacteria and algae display T&O signature communities, where specific taxa are associated with T&O concentration categories. From combining data from these communities, indicative taxa associated with T&O-risk were identified in a reservoir experiencing extreme geosmin concentrations (Reservoir 1 – ≤520 ng L-1) and in a reservoir with moderate 2-MIB concentrations (Reservoir 2 – ≤60 ng L-1). Indicative taxa for the differentiation between low and high T&O metabolite concentrations were identified and T&O risk categories were constructed based on the shared ecological functions of each taxon. These categories were further applied to include the taxa associated with T&O-producing Cyanobacteria from community network analysis. Here, five T&O risk categories are proposed and are summarised below; the taxa associated with each category are detailed in Table 3 for reference.

4.1.1. Cyanobacterial T&O Producers

For a T&O event to occur, Cyanobacteria that possess the genes encoding geosmin synthase (geoA) and/or 2-MIB cyclase (mic) must be present within the community under appropriate environmental conditions, such as low NO₃^- and NO₂^- concentrations as observed in this study. Aphanizomenon NIES81 was an important variable in the RF model analysis for categorising geosmin concentrations. RF results from Reservoir 2 indicated that Cyanobium PCC-6307 was associated with the presence of 2-MIB, suggesting that this Cyanobacteria could be responsible for 2-MIB production. In this study, both Aphanizomenon NIES81 (Reservoir 1) and Cyanobium PCC-6307 (Reservoir 2) exhibited more lower-than-expected co-occurrences with other community members than higher-than-expected co-occurrences. These lower-than-expected associations imply biotic stress (induced by unstable ecological associations) may trigger the production of T&O metabolites, with geosmin and 2-MIB potentially serving as a defence mechanism, as sesquiterpenoids have previously been implicated in cellular defence (Pattanaik & Lindberg, 2015). In Reservoir 1, Aphanizomenon NIES81 was associated with low NO₃^-, NO₂^-, and high geosmin concentrations, consistent with previous findings indicating maximal geosmin production when oxidised nitrogen forms are scarce (Harris et al., 2016). Thus, the presence of T&O-producing Cyanobacteria is likely influenced by biotic stress and/or specific nutrient conditions, such as a low NO₃^-:NH₄⁺ ratio, which may stimulate both cyanobacterial growth and T&O metabolite production (Hooper et al., 2023).

4.1.2. T&O Degraders

Elevations in T&O metabolite compounds are hypothesised to select for different bacterial species capable of degrading these compounds, thereby increasing the likelihood of such degradation processes occurring. Clercin et al. (2021) observed elevations in geosmin concentrations in a eutrophic reservoir were associated with Novosphingobium hassiacum and Sphingomonas oligophenolica, whereas increased concentrations of 2-MIB were linked to Flavobacterium species. The roles of Pseudomonas and Bacillus in T&O compound degradation remained unclear (Clercin et al., 2021), despite evidence of their capacity to degrade both compounds (Eaton & Sandusky, 2010; Ho et al., 2007; Ma et al., 2015). In this study, the geosmin-degrading Pseudomonas had a greater influence on the RF classification of high rather than low geosmin concentrations in Reservoir 1 (geosmin). The mean relative abundance of Pseudomonas increased during times of high geosmin concentrations compared to low, implying that increases in this degrader could indicate a T&O risk. Three strains of the T&O degrader, Pseudomonas, have been reported to also contain the geosmin synthase gene (geoA), although it is not known if they are effective geosmin producers (Churro et al., 2020). It is unknown why most species of Pseudomonas lack the genes to produce geosmin. In contrast, Bacillus in Reservoir 1 (geosmin) occurred at a lower-than-expected frequency with Aphanizomenon NIES81 and was an important variable in the RF model analysis for categorising geosmin concentrations. T&O degraders like Bacillus have been shown to display an antagonistic relationship with a previously described geosmin-producing strain of Streptomyces (Ma et al., 2015), which agrees with previous findings from natural communities (Bertos-Fortis et al., 2016; Salmaso et al., 2018). The increased abundance of Bacillus would thus be characteristic after periods of high T&O. A potential 2-MIB degrader identified in this study that needs further investigation is Limnobacter. In Reservoir 2 (2-MIB), RF analysis revealed Limnobacter as the most informative genus in categorising 2-MIB concentrations into low and high categories, with heightened relative abundance in Limnobacter during high 2-MIB concentrations. It has been postulated that Limnobacter species can degrade a vast array of aromatic compounds (Pérez-Pantoja et al., 2012; Vedler et al., 2013), although there are no reports to date of this genus degrading T&O compounds.

T&O degraders have the potential to cause fluctuations in geosmin and 2-MIB concentrations in water columns, depending on the relative rates of cyanobacterial T&O production and the biodegradation of the T&O metabolites.

4.1.3. Nutrient Enrichment Indicators

Nutrient enrichment indicators are taxa that tend to be present during periods of nutrient enrichment. For instance, the Peptostreptococcaceae family has been identified as an indicator of nutrient enrichment in urban lakes (Numberger et al., 2022). In Reservoir 1 (geosmin), Peptostreptococcaceae exhibited a significantly lower-than-expected co-occurrence with Aphanizomenon NIES81, despite having the greatest influence on the RF model’s performance in classifying low versus high geosmin concentrations. The enrichment of nutrients such as NH₄⁺ could trigger the cyanobacterial production of T&O metabolite compounds like geosmin (Perkins et al., 2019). The presence of Peptostreptococcaceae before a T&O event suggests a potential early indicator of T&O risk. The lower-than-expected co-occurrence with Aphanizomenon NIES81 in Reservoir 1’s community network could be explained by niche overlap, as both taxa are capable of nitrogen fixation (Paster et al., 1993). This overlap may lead to competitive exclusion, where these organisms compete for similar ecological resources, potentially reducing their co-occurrence within the same habitat. Another indicator of nutrient enrichment is the diatom genus Melosira, which is commonly associated with eutrophic and polluted aquatic systems (Van de Vijver & Crawford, 2020). In the studied reservoirs, Melosira played a significant role in the RF analysis for classifying low and high geosmin concentrations. Additionally, Melosira exhibited a higher-than-expected co-occurrence with the 2-MIB-producing cyanobacterium Dolichospermum NIES41 in Reservoir 2 (2-MIB), suggesting a potential ecological interaction between these taxa under nutrient-enriched conditions. Armatimonas was also significant in the RF classification of low and high geosmin concentrations in Reservoir 1, consistent with previous studies linking this taxon to geosmin production after increases in TP concentrations (Xu et al., 2018). Elevations in TP concentrations favour cyanobacterial growth and productivity, and consequential shifts in the total nitrogen (TN) to TP ratio – particularly TN ratios below 30:1 – have been associated with increased geosmin production in reservoirs (Harris et al., 2016).

Additionally, the diatom Nitzschia, known for its association with nutrient-enriched waters and sensitivity to nutrient depletion (Reynolds et al., 2002), was influential in the RF classification of low and high 2-MIB concentrations in Reservoir 2. In Reservoir 2 declines in Nitzschia abundance were indicative of heightened 2-MIB risk. Additional diatoms characteristic of nutrient enrichment are Aulacoseira and Asterionella (Reynolds et al., 2002). Aulacoseira displayed lower-than-expected co-occurrence with Aphanizomenon NIES81 in Reservoir 1 (geosmin), however, according to Reynolds et al. (2002), this is to be expected due to their sensitivity to rising pH values associated with increased cyanobacterial productivity (Slavin, 2020). In contrast, Asterionella had a higher-than-expected co-occurrence with Cyanobium PCC-6307 in Reservoir 2 (2-MIB), rapid increases in Asterionella species can deplete major anions (NO₃^-, PO₄^3- and silicates), leading to altered nutrient ratios (Krivtsov et al., 2000).

4.1.4. Nutrient Exchange/Recycling

Taxa capable of nutrient exchange and recycling can benefit Cyanobacteria by creating favourable conditions that promote the production of secondary metabolites, such as geosmin and 2-MIB (Wu et al., 2024). For example, in Reservoir 2, Limnohabitans had a higher-than-expected co-occurrence with Cyanobium PCC-6307 and was highly influential in classifying low and high 2-MIB concentrations. The primary growth substrates for Limnohabitans are autochthonous algal/ cyanobacterial organic matters and products from the photolysis of dissolved materials (Ye et al., 2016); culture-independent analyses show this genus to constitute 11% of the cyanobacterial phycosphere (Kim et al., 2020). The NS9 marine group was also an important variable in categorising low and high 2-MIB concentrations in Reservoir 2. The growth of these prokaryotes is associated with phytoplanktonic organic matter (Takebe et al., 2024) and have been identified to degrade polysaccharides (Reintjes et al., 2023). Similarly, previous work has identified Gemmataceae to contain wide repertoires of carbohydrate-active enzymes allowing them to utilise some polysaccharides, such as xylan, laminarin, lichenin and starch (Kulichevskaya et al., 2020). Gemmataceae had a lower-than-expected co-occurrence with Cyanobium PCC-6307 in Reservoir 2 but was influential in categorising 2-MIB concentrations into low and high categories. Another taxon affiliated with the nutrient exchange/ recycling for 2-MIB-producing Cyanobium PCC-6307 was Rhizobiales A0839. Previous studies link this taxon with NH₄⁺ oxidation and removal in wastewater treatment (Rodriguez-Sanchez et al., 2019). This is an unusual pairing as both taxa compete for NH₄⁺, although previous evidence suggests that NH₄⁺-oxidising bacteria deny Cyanobacteria access to NH₄⁺ and force diazotrophic Cyanobacteria into nitrogen fixation mode (Boyett et al., 2013). However, forced cyanobacterial nitrogen fixation in Cyanobium PCC-6307 cannot be assumed to be fuelling 2-MIB production in Reservoir 2 without further investigation, as the majority of 2-MIB-producing Cyanobacteria are thought to be non-heterocyst producing (Jeong et al., 2021).

The only higher-than-expected co-occurring taxon with Aphanizomenon NIES81 in Reservoir 1 was Fluviicola. Fluviicola is affiliated with the family Cryomorphaceae which belongs to the class Flavobacteriia and contains many species that play an integral role in the recycling of carbon and energy in freshwater environments (Woyke et al., 2011).. Flavobacteriia are major decomposers of high-molecular-mass organic matter in water bodies (Cottrell & Kirchman, 2000), like the 2-MIB degrading genus Flavobacterium (Clercin et al., 2021). Guedes et al. (2018) have previously shown Fluviicola to have positive correlations with Synechococcus; following previous findings from natural communities (Bertos-Fortis et al., 2016; Salmaso et al., 2018). The Chlamydomonadaceae family co-occurred at a higher-than-expected probability with the Nostocaceae family in Reservoir 1, this family contains the genus Chlamydomonas. Although no reports suggest that the Nostocaceae family can produce cylindrospermopsin, Kust et al. (2018) stated that Nostocaceae are a rich source of unknown secondary metabolites that require further investigation. These metabolites could have a similar effect to cylindrospermopsin on members of the Chlamydomonadaceae family, inducing the algae to produce APases for the Cyanobacteria to use.

4.1.5. Cyanobacterial Benefactors

Sediminibacterium has been shown to offer Cyanobacteria protection against the colonisation of opportunistic bacteria by producing bacteriocin and toxoflavin (Sethuraman et al., 2022). Sediminibacterium can also form biofilm consortia with Cyanobacteria (Velichko et al., 2015), and its relationship with Cyanobacteria in wastewater has also been previously reported (Lee et al., 2013). In this study, Sediminibacterium was considered more influential on the RF model’s performance in classifying high than low 2-MIB concentrations. Another taxon influential in categorising 2-MIB concentrations was the genus OM190. This taxon produces a diverse range of secondary metabolites including antimicrobial compounds (Pushpakumara et al., 2022), which have the potential to protect Cyanobacteria from undesired microbes. Protection of T&O-producing Cyanobacteria would in turn enhance T&O risk. For the categorisation of low and high geosmin concentrations, the OM27 clade was highly influential in the model’s performance, especially for high geosmin concentrations. Members of the OM27 clade are known to prey on opportunistic bacteria (Orsi et al., 2016) and have been shown to dominate the bacterioplankton community during cyanobacterial blooms (Wang et al., 2020). We hypothesise that members of the OM27 clade predate on opportunistic bacteria alongside Cyanobacteria, enabling them to become more resilient in the community, which in turn can present a T&O risk.

5. Conclusion

This study marks the first to investigate T&O risk linked to Cyanobacteria by integrating bacterial and algal data to identify indicator taxa for T&O-risk categories. The analysis of bacterial and algal communities in reservoirs with elevated geosmin and 2-MIB concentrations revealed distinct T&O signature communities. The identification of bacterial and algal communities associated with low, medium, and high T&O concentration categories suggests that T&O risk can be inferred from these communities. By combining bacterial and algal datasets, key taxa distinguishing between low and high T&O concentration categories were identified.

T&O-producing Cyanobacteria in Reservoir 1 (geosmin) and Reservoir 2 (2-MIB) showed fewer dependable associations with other taxa, indicating a selective or unstable ecological role. This instability may trigger T&O metabolite production as a stress response under specific environmental conditions, such as low NO₃⁻ and NO₂⁻ levels. Conversely, Cyanobacteria in the reservoir without T&O events exhibited more stable interactions with other taxa, suggesting a more integrated ecological role.

These findings underscore the importance of cyanobacterial interactions within the broader microbial community and the potential of using indicator taxa to assess T&O risk. This approach provides valuable insights for water companies and environmental managers, enabling more accurate predictions and management of T&O events. Importantly enhanced predictive capability also better identifies triggers of T&O metabolite production prior to events occurring after a temporal lag (between metabolite production and T&O event at the Water treatment works) and are hence of key importance in drinking water resource management.