Association between Legionella species and humic substances during early summer in the northern Baltic Sea

Introduction

Climate change is projected to cause alterations in the Baltic Sea, such as humification that leads to decreasing salinity and increased water temperatures. This can lead to ecosystem imbalances and the emergence or establishment of new bacterial species. Several bacterial pathogens, such as Legionella pneumophila, has its natural lifecycle in the environment and can live inside protozoan host cells. This ecological niche of predation resistant bacteria are relevant for public health, since some species can cause disease in humans (Thomas et al., 2010).

Bacterial defense strategies has evolved due to the co-existence between bacteria and their predators (protozoa). Some of these strategies contribute to the ability of bacteria to cause human disease. Francisella tularensis, Legionella pneumophila, Vibrio cholerae and Yersinia pestis all are examples of pathogenic bacteria that have anti-predator strategies. The mechanisms bacteria use for both of these traits are often similar (Sun et al., 2018; Amaro and Martín-González, 2021). To understand the emergence of diseases from opportunistic or pathogenic bacteria, it is important to complement studies of bacteria in association to the susceptible host (human cells) with studies of bacteria in the environment of their natural lifecycle. It is likely that bacterial traits that result in pathogenicity to humans originally arose as a response to environmental stress such as predation (Park et al., 2020). For example, Legionella is a predation resistant bacteria that reside in the environment, and sporadically cause problems in water supplies (Hamilton et al., 2018; Leoni et al., 2018; Holsinger et al., 2022). Still, there is a knowledge gap regarding the occurrence and persistence of pathogenic and predation resistant bacteria like Legionella in natural environments including water.

What allows public health relevant species within the genus Legionella to survive in environmental water is unknown, even though more than 20 out of 71 reported species within Legionella have been recognized as human pathogens and the causative agent of Legionnaires’ disease and Pontiac fever (Muder and Yu, 2002; Legionella Species - Special Pathogens Laboratory, 2020; Parte et al., 2020). L. pneumophila is the most studied species in human disease as well as in environmental sampling. Other Legionella non-pneumophila species that cause disease are less studied and undiagnosed (Pascale et al., 2021). Thus, there is a need to investigate non-pneumophila species (Duron et al., 2018).

Bacterial species within Legionella has been isolated from various types of environments such as soils and drinking water (Chambers et al., 2021). Even though freshwater is considered the major reservoir of L. pneumophila (Boamah et al., 2017), the bacterium is mostly studied in man-made water systems such as cooling towers (Brigmon et al., 2020), where the thermophilic bacteria have been shown to dwell (Lesnik et al., 2015). Only a few studies have reported occurrence of DNA sequences specific to the genus Legionella in natural aquatic environments such as in lakes, rivers, brackish and marine waters (Carvalho et al., 2007; Parthuisot et al., 2010; Shimada et al., 2021). Thus, there is a need to increase our understanding of the natural diversity of Legionella spp. in aquatic systems. One approach to understand driving mechanisms of aquatic bacteria is mapping them across dynamic habitats. Freshwater and marine systems are especially dynamic in coastal zones that function as filters for chemical compounds and organisms entering the ocean through rivers (Lamb et al., 2017). Due to this inflow, the coast is generally more nutrient rich compared to the offshore areas. Thus, bacteria that inhabit these interfaces between freshwater and marine systems are exposed to stresses/disturbances such as nutrient rich river inflow, which may lead to fluctuations in selection dynamics.

The Baltic Sea is at risk for large fluctuations, maybe even to the scale of the structure and function of the ecosystem (Blenckner et al., 2015). Projections indicate that the Baltic Sea will suffer alterations due to climate change, both due to increased temperature and due to increased inflow of organic matter caused by precipitation (Andersson et al., 2015; Harvey et al., 2015; Meier et al., 2022). Bacteria can use humic rich organic substances as food resource, giving them an ecological advantage over other plankton organisms. This, in turn, will cause lowered planktonic food quality as for example high C to N ratios and poorer fatty acid composition (Bandara et al., 2022; Guo et al., 2022). Previous research show that in the north of the Gulf of Bothnia phosphorus is the most limiting nutrient, while in the south nitrogen limits organismal growth (Andersson et al., 1996; Tamminen and Andersen, 2007). These habitat variations can select for different bacterial species. The brackish Baltic Sea provides a unique opportunity to study drivers along natural environmental gradients, where the northern parts are more similar to freshwater and the southern parts are more saline.

In the middle of the Gulf of Bothnia, along a 10 km coastal area, a previous study reported widespread occurrence of Legionella spp. in early summer (Eriksson et al., 2022). To map occurrence of Legionella spp. on a larger spatial scale in the Gulf of Bothnia and with the aim to find environmental drivers, we here analyzed samples in the coast-offshore and north-south gradients in early summer. We used two complementary methods for a general bacterial identification to capture the wide ecological niche of predation resistant bacteria. Then, we particularly focused on Legionella, with the aim to also resolve the species. The methods used were sequencing the V4 region of the 16S rRNA gene (using Illumina sequencing) and the V1-9 region (using Nanopore sequencing). To support the species classification, we amplified a L. pneumophila specific macrophage infectivity potentiator (mip) gene that is associated with virulence (Ballard et al., 2000; Nazarian et al., 2008; Collins et al., 2015; Yin et al., 2022). We expected that relatively high concentrations of humic substances in the coast and the open northern basin of the Gulf of Bothnia (Bothnian Bay) would promote bacterial production and their predators (protozoa), which in turn constitute a niche for Legionella spp. In the southern open basin, Bothnian Sea, the humic substance concentration is lower and thus we expected lower bacterial production, as well as lower abundance of protozoa and Legionella spp.

Materials and methods

Study sites and environmental data

In June 2019, water samples were collected at 9 different stations (5m depth) in the Gulf of Bothnia: Råneå 1 (RA1), Råneå 2 (RA2), A5, A13, Örefjärden (B3), Gavik (GA1), C3, C14 and C24 (Figure 1). Stations RA1, RA2, A5 and A13 are located in the Bothnian Bay basin (northern Gulf of Bothnia). Stations B3, C14, GA1, C3 and C24 are located in Bothnian Sea basin (southern Gulf of Bothnia). Stations RA1, RA2, GA1, B3 and C14 are coastal sampling sites. Coordinates can be found in Supplementary Material 2.

Figure 1 Map of the stations in the Gulf of Bothnia. Northern stations RA1, RA2, A5 and A13 are located in the Bothnian Bay basin (northern Gulf of Bothnia). Southern stations B3, C14, GA1, C3 and C24 are located in the Bothnian Sea basin (southern Gulf of Bothnia). Costal stations: RA1, RA2, B3, C14 and GA1. Offshore stations: A5, A13, C3 and C24.

Umeå Marine Science Centre (UMSC) collected the environmental data via the Swedish National Environmental monitoring program. Data were analyzed according to the HELCOM Manual for Marine Monitoring in the COMBINE report (Cooperative Monitoring in the Baltic Marine Environment).

The included variables for this study were conductivity, temperature and depth (CTD), bacterial production, humic substances, nutrient salts, dissolved organic carbon (DOC), total nitrogen and total phosphorous, collected at 5m depth. For Chlorophyll-a and primary production, a pooled fraction of the surface water (1-10m) was used. Procedures for the measurements are presented at HELCOM (HELCOM, 2017). The data was accessed at SHARKdata (SHARKdata - Datasets, 2022). Differences in environmental concentration levels were estimated using linear models fitted within R, comparing the Bothnian Bay (north) and the Bothnian Sea (south).

Bacterial DNA sampling and extraction

At each sampling site, 75-500 ml water was gently filtered (≤20kPa) onto 0.2 µm filters (Pall Corporation sterilized filters, Supor ^® 0.2 µm, 47 mm, S-pack white gridded). For samples from RA1 and RA2, <500ml was filtered due to limitations with the filtering of turbid water (100 ml and 250 ml, respectively). A total of 9 water samples were collected. The filters were folded using cleaned tweezers and placed in 2 ml Eppendorf tubes. The samples were stored in -80°C until the DNA extraction. Filters were thawed and placed into Powerwater bead beating tubes. The DNA was extracted using a DNeasy^® PowerWater^® kit (Qiagen, Hilden, Germany) according to a modified DNeasy PowerWater protocol, where the samples were treated with an additional heating step (horizontal water bath for 30 min, 65°C) and with a bead-beating step in each direction of 20 Hz, 3x3min with a Tissuelyser II (Qiagen, Hilden, Germany). The concentration of the DNA extracts from the stations ranged from 6-46 ng/µl. Two filter blanks were used as negative control during DNA extraction. The final DNA was frozen before subjected to PCR analysis. Sample preparation, thermal cycling and PCR product preparation were performed in separate rooms.

Illumina amplicon preparation and sequencing

For the 16S rRNA amplicon preparation, DNA was amplified and sequenced as previously described (Hägglund et al., 2018), part from the PCR purification kit used. In short, the DNA was amplified using the No. 5 Hot Mastermix 2.5x kit (5 PRIME, Qiagen, Hilden, Germany) with bacteria/archaeal primers 515F/806R specific for the hypervariable V4 region of the 16S rRNA gene (Caporaso et al., 2012). BSA has been shown to reduce PCR inhibition from humic substances (Opel et al., 2010; Sidstedt et al., 2020). To avoid inhibition, 8 µg of BSA was therefore used for each reaction. The forward and reverse primers were modified to incorporate a 12 bp Golay error-correcting barcode that enables sample multiplexing (Caporaso et al., 2012). All samples were amplified in triplets and pooled after PCR amplification (94°C for 3 min; 35 cycles of 94°C for 45 s, 50°C for 1 min, 72°C for 1.5 min; 10 min rest to finish). For PCR, 3 template blank reactions were included by using nuclease free water instead of template.

For the 18S rRNA amplicon preparation, DNA was amplified using the No. 5 Hot Mastermix 2.5x kit (5 PRIME, Qiagen, Hilden, Germany) with eukaryotic primers V6F/V8R (V6F: 5’- AATTYGAHTCAACRCGGG, V8R: 5’-GACRGGCGGTGTGNACAA) specific for the hypervariable V6-V8 region of the 18S rRNA gene. These primers were designed by Eriksson et al. (2022), using higher degeneracies than earlier publications (Edgcomb et al., 2011; Hadziavdic et al., 2014). The benefits of this region is described in detail by Latz et al. (2022). All samples were amplified in triplets and pooled after PCR amplification (94°C for 3 min; 35 cycles of 94°C for 45 s, 57°C for 1 min, 72°C for 1.5 min; 10 min rest to finish).

The PCR product was run on a 1% agarose gel and the DNA concentration was estimated with a Qubit fluorometer (Invitrogen, Carlsbad, CA, United States). The amplicons were pooled at equimolar concentrations and purified with the QIAquick PCR purification kit (Qiagen, Hilden, Germany) following the supplier’s instructions. The DNA concentration of the pooled amplicon product was measured with a Qubit fluorometer and adjusted to 2 nM. The library was denatured and diluted as described by Illumina (MiSeq System User Guide, Part # 15027617 Rev. C), before it was loaded onto a MiSeq cartridge (Illumina, San Diego, CA, United States) and sequenced using a 2x300 bp paired-end sequencing protocol (Hägglund et al., 2018).

Illumina quality control and raw data processing

Demultiplexing was done using deML (Renaud et al., 2015). The Quantitative Insights Into Microbial Ecology (QIIME2) pipeline, version 2019.1, was used for processing raw sequence data (Bolyen et al., 2019). Greengenes version 13.8 (McDonald et al., 2012) and PR2 version 4.12.0 (Guillou et al., 2012) databases were used for bacteria and eukaryote taxonomic assignment. Quality filtering was done using dada2 default parameter values. For the 16S rRNA sequence raw data processing, the reads were trimmed to 275 bp for the forward read and to 250 bp for the reverse read. For the 18S rRNA sequence raw data processing, the forward read were trimmed to 290 bp allowing a maximum of 5 expected errors per read. For 18S rRNA, the reverse read were not used due to low quality.

The Illumina 16S rRNA amplicon sequencing yielded 10549128 reads, 5211659 (49.4%) of which remained after the data was processed into amplicon sequence variants (ASVs). The data was filtered, denoised, merged and chimeras were removed, resulting in 4497 ASVs; the number of reads per sample ranged from 242562 to 562960. The Illumina 18S rRNA amplicon sequencing yielded 323690 reads, 63790 of which remained after the data was processed (19.7%, forward read only), resulting in 97 ASVs; the number of reads per sample ranged from 11815 to 2593.

Nanopore amplicon preparation and sequencing

To increase the taxonomic resolution and to compare with Illumina generated data, longer amplicons were also prepared and sequenced using the Oxford Nanopore Technology (ONT). In short for the 16S rRNA amplicons, the DNA was amplified using the Terra PCR Direct Polymerase Mix (Clontech, Mountain View, CA) with bacteria/archaeal primers with ONT overhang (S-D-Bact-0008-c-S-20 (27F): 5’-AGRGTTYGATYMTGGCTCAG-3’, S-D-Bact-1492-a-A-22 (1492-R): 5’-TACGGYTACCTTGTTACGACTT) (Lane, 1991) amplifying the V1-V9 region of the 16S rRNA gene. To avoid inhibition, 8 µg of BSA was used for each reaction. All samples were amplified in duplicates using 30 and 33 cycles (98°C for 2 min; 30 or 33 cycles of 98°C for 10 s, 57°C for 15 s, 68°C for 1.5 min; 68°C for 5 min rest to finish).

Single stranded DNA were removed from the PCR products using 2 µl ExoSAP-IT (Applied Biosystems, USA) to ~150 fmol/sample and each purified sample were then quantified using the Qubit 1×HS Assay Kit (Thermo Fisher Scientific).

The PCR Barcoding Expansion kit (EXP-PBC096) together with Ligation Sequencing Kit (SQK-LSK109) were used to sequence the 16S rRNA amplicons PCR. After the PCR-barcoding step in the protocol a purification of PCR products were performed using 0.5× of AMPure magnetic beads (Beckmann) following the manufacturer’s protocol, thereafter, the purified PCR products were quantified using the Qubit1×HS Assay Kit (Thermo Fisher Scientific), and were pooled in equimolar amounts (200 ng/sample) before proceeding with the protocol. Some modifications were made within the protocol to adjust the volumes for the following sequencing on a MinION Mk1B instrument (ONT) with a Flongle Flow Cell (R9.4.1 chemistry, FLO-FLG001). The MinION instrument was ongoing for approximately 48 hr, or until no further reads could be collected.

Nanopore quality control and raw data processing

Basecalling and demultiplexing were performed using Guppy v6.0.1 (Oxford Nanopore Technologies). The data processing and read classification steps were implemented using the Snakemake workflow management system (Koster and Rahmann, 2012) together with Bioconda (Grüning et al., 2018). For visualization, nanoplot was used (version 1.29.1, https://github.com/wdecoster/NanoPlot/releases). Quality control was done using nanoq (version 0.8.6, https://github.com/esteinig/nanoq). The nanopore data was quality trimmed using Porechop (version 0.2.4, https://github.com/rrwick/porechop). Reads were then filtered by length (0-1 Mb) with Nanofilt (version 2.8.0, https://github.com/wdecoster/nanofilt). We assessed read statistics using NanoStat (version 1.4.0, https://github.com/wdecoster/nanostat) (De Coster et al., 2018). Overall, we obtained 834343 reads, with a mean read quality of 11.7. Otherwise, default settings were used. Taxonomic classification of the consensus sequences was done by kraken2 version 2.1.2 (Wood et al., 2019) and bracken version 2.6.2. (Lu et al., 2017) using a custom database and GTDB v207 (a list of included GCF files and taxonomy can be found in Supplementary material 2). EMU (version 2.0) was used using the original database (Curry et al., 2022).

Target organisms

To include as many of the possible species of Legionella, the ASVs were selected at a higher taxonomic level for the Illumina dataset (at family level Legionellaceae). Phylogenies were estimated based on FastTree 2 (Price et al., 2010). For the Nanopore dataset, the genus Legionella was selected.

In order to compare the co-occurrence of Legionella and include heterotrophic and mixotrophic protists, other phylogenetic groups were excluded from the eukaryotic data set based on knowledge of feeding style (Olenina et al., 2006). The supergroup Archaeplastida were excluded. For the supergroup Opisthokonta, Choanoflagellates were kept as they are filter feeders (Richter and Nitsche, 2017), while the Divison of Metazoa and Fungi were excluded. For the supergroup Stramenopiles, heterotrophic and mixotrophic Chrysophyceae was included since they have phagotrophic feeding style. From the Division Ochrophyta, class Chrysophyceae, Dictyochophyceae and unknown assignment (NA) were kept while non-relevant classes were excluded (Bolidophyceae, Bacillariophyta, MOCH-2, Picozoa, Phaeothamniophyceae, MOCH-5 Raphidophyceae, Synurophyceae, Eustigmatophyceae, Phaeophyceae, and Xanthophyceae). The sequence data was visualized using R package phyloseq (McMurdie and Holmes, 2013) and ggplot2 (Wickham, 2016).

Microscopic analysis

To obtain the protist biomass, microscopic counting and analysis of phytoplankton from the environmental monitoring program was included in this study. Phytoplankton were analyzed microscopically using the Utermöhl method according to the Helcom Combine manual (Utermöhl, 1958; Olenina et al., 2006; HELCOM, 2017). Autotrophic and mixotrophic phytoplankton were identified to the highest possible taxonomic level. Most heterotrophic protists were however not counted. The size of the cells was measured and their biovolume calculated (Olenina et al., 2006). Carbon biomass was calculated using biovolume to carbon biomass conversion factors (Menden-Deuer and Lessard, 2000). This microscopic method is constrained to nano- (2-20 µm) and micro scale (20-200 µm), and thus picoeukaryotes (<2 µm) are not included in the analyses.

Phylogenetic analysis of sequences assigned to Legionella pneumophila

To investigate the validity of the classification of the ASVs assigned to Legionella pneumophila by kraken, a multiple sequence alignment and phylogentic tree was created for the environmental ASVs and Legionella spp. The similarity of Legionella spp. 16S rRNA sequences was determined using the Basic Local Alignment Search Tool (BLASTN) at the National Centre for Biotechnological Information (NCBI) and EzTaxon-e server (Kim et al., 2012). The 16S rRNA of the closest hit was retrieved from the IMG-JGI database (Chen et al., 2021). The SBDI Sativa curated 16S GTDB database (version 5) (Lundin and Andersson, 2021) was downloaded and sequences annotated as Legionella were extracted to construct the phylogeny. The extracted sequences were concatenated with sequences selected in this study and the sequence of Francisella tularensis subsp. tularensis SCHU S4 was selected as the outgroup. All concatenated sequences were aligned using MAFFT (Katoh and Standley, 2013) and alignment was trimmed using ClipKIT with the smart-gap option (Steenwyk et al., 2020). The evolutionary relationship was inferred by the Maximum Likelihood method using RAxML-NG (Kozlov et al., 2019) with model GTR+G with Transfer Bootstrap Expectation (TBE) calculation (Lemoine et al., 2018) and visualized using Figtree (version 1.4.4, https://github.com/rambaut/figtree/).

Generalized linear latent variable model

A generalized linear latent variable model (GLLVM) was used to estimate the importance of the environmental factors on the occurrence on Legionella spp. To include as many of the possible species, the ASVs were selected at a higher taxonomic level (at family level Legionellaceae). For accurate frequency distribution of ASVs, the Illumina dataset was used for the model (a more error prone nanopore dataset lead to many unique ASVs). The GLLVM was fitted to the data using the R-package gllvm (Niku et al., 2019). Gaussian distribution function (i.e. identity link) for the responses was selected: environmental variables was included in the model as continuous predictors: humic substances, TDP and temperature. To avoid including highly correlated variables as predictors in the model, a pre-selection step was performed based on the PCA scoreplot (Figure 2), of the environmental variables (i.e. only one of a group of tightly clustered variables was selected, Figure 2). This resulted in a model with fewer parameters that, according to the theory of Occam’s razor, is preferred since they are easier to fit and less prone to convergence issues compared to complex models. For example, temperature and bacterial production were highly correlated (correlation: 0.90), thus bacterial production was excluded from the statistical analyses. In the analysis, it is therefore difficult to disentangle the effect of these variables on the abundance of Legionellaceae. To identify how much of the ASV variation each variable could explain, a separate model was first made for each variable (Table 1). For the selection of the final model, a sensitivity analysis was performed, where the model with the best model of fit and lowest AIC and BIC scores was chosen, i.e. the most optimal trade-off between model fit and complexity (Table 1). Finally, based on the model selection results, humic substances, TDP and temperature were included in the model along with a factor for coastal habitat and a factor for filter volumes (to compensate for when less than 500 ml of water could be filtered).

Figure 2 PCA of the spatial distribution of environmental variables at the stations. The prcomp function was used in the R package stats with default settings, using normalized and scaled data.

Table 1 Variation explained in the Amplicon Sequence Variant (ASV) communities for each environmental variable (according to the ratio of traces).

The continuous variables were normalized so that the mean value and standard deviation were set to zero and one, respectively, for each included predictor. To improve convergence, jitter variance for starting values of latent variables were set to 0.5, and the number of initial runs were set to five. Since additional latent variables did not improve the model, one latent variable was selected to capture the residual variation. Otherwise, default values of the GLLVMs function were used. The identified Legionella occurred at relatively low abundances in the sequence data, constituting of <1% of the entire bacterial community. To manage rare sequence variants and avoiding bias, and to reduce noise in the models, rare sequence variants that only occurred in one sample were excluded in the model.

Real time PCR - Legionella specific amplification

Primer and probe sequences for the L. pneumophila-specific mip gene target were used (Nazarian et al., 2008). The L. pneumophila specific TaqMan^® probe had a 6-carboxy-fluorescein (FAM) as the fluorescent reporter on the 5′ end and a 6-carboxytetramethylrhodamine (TAMRA) as the quencher dye on the 3′ end of the probe (Nazarian et al., 2008). SsoAdvanced Universal Probes Supermix was used (Bio-Rad). To avoid inhibition, 8 µg of BSA was used for each reaction. For each reaction, 1 μL of template DNA was used for a total volume of 20 μL. Thermal cycling, fluorescent data collection, and data analysis were carried out on an iCycler (Bio-Rad Laboratories, Hercules, CA). The thermal cycling profile consisted of an initial incubation at 95°C for 3 min, followed by 45 cycles of 95°C for 5 s and 60°C for 30 s.

Primer targeting a 16S rRNA gene region specific for Legionella were used. For each reaction 1 µl of template DNA was used in a mix of SsoAdvanced Universal Probes Supermix (Bio-Rad), BSA and primer iQLeg1F 5’-CACTGTATGTCAAGGGTAGGTAAG-3’ and iQLeg1R5’-TTAGTGGCGCAGCAAACGCGAT-3’ (this study) for a total volume of 20 µL. The thermal cycling profile consisted of an initial incubation at 95°C for 3 min, followed by 40 cycles of 95°C for 5 s and 60°C for 5 s. Water samples spiked with L. pneumophila concentrations ranging from 10 to 10⁵ colony forming units per mL were analyzed in duplicates to test the limit of detection and generate a standard curve for assessing target DNA concentrations. The standard curve was evaluated with the iQLeg1 primer pair. Positive control mixtures, using DNA from L. pneumophila and negative control mixtures without a template, were included in each PCR run.

Results

Environmental conditions

Principal component analysis (PCA) of the environmental variables show that the northern parts of the Gulf of Bothnia was characterized by high DOC and humic substances while the southern parts was clearer, had relatively higher salinity and higher pH (Figure 2). In addition, the spatial distribution of environmental variables was larger in the northern parts of the Gulf of Bothnia compared to the south. There were several significant environmental differences between the north and the south: for humic substances, pH, TDN, chlorophyll-a, and salinity (Tables S1, S2 and Supplementary Figure S1).

The concentration of humic substances were highest at station RA1, RA2, and B3 (Figure 3), where the bacterial production were also high compared to the rest of the stations (>8 compared to <6 µgC/(dm3*day)). Nitrogen concentrations were higher in the north (mean north 19.21 +/- 2.65, mean south 15.22 +/- 1.00, p<0.05, Supplementary Figure S1; Tables S1, S2), especially at the coastal stations (RA1: 23.42 µM and RA2: 19.42 µM), followed by the offshore stations in the north (A13: 17.21 µM and A5: 19.28 µM). Phosphorus concentrations were especially low at station A13 (0.18 µM), while station RA1 and B3 had the highest concentrations of phosphorous (>0.4 µM). The temperature were the highest in the coastal northern RA1 and RA2 stations, followed by the southern C24, B3, C14 and C3 stations.

Figure 3 Environmental variables at the stations.

Occurrence of Legionella in the Gulf of Bothnia

Short amplicon (Illumina) sequence variants assigned to genus Legionella were only detected in the northern samples (RA1, RA2, A13, A5 and B3), (Supplementary Figure S2, top). A total of 2 amplicon sequence variants (ASVs) were identified that were assigned to genus Legionella (Supplementary Figure S3, top). Legionellaceae was however detected in all samples from north to south. In total, 24 ASVs were detected (Supplementary Figure S2, bottom), 14 occurred in more than one sample (Figure 4). Station A5 had the highest relative abundance of ~0.3%, followed by station A13 and B3 with relative abundances of about 0.25% (Supplementary Figure S2, top). The relative abundance in the remaining samples was <0.15%. The relative abundance at A5 was about four times higher compared to the stations with the lowest relative abundance (RA1, C24 and C3).

Figure 4 Filtered Illumina Legionellaceae ASVs, detected in more than one sample (ASVs used in GLLVM).

For the long amplicon sequencing (the nanopore dataset), species were assigned by classification algorithm emu, kraken and bracken. Some differences between the classifiers were observed. For the genus Legionella, 5 species were suggested by EMU (Figure 5, top), where L. pneumophila were the most frequently suggested species (sample A13, A5, B3, GA1 and RA2). According to kraken, 16 species were suggested (Figure 5, middle), where L. waltersii were the most frequently suggested species (sample A13, A5, B3 and RA2). Bracken assigned 20 species (Figure 5, bottom), where L. pneumophila were one of the most frequently suggested species (sample A5, B3, C14 and GA1). Correlations between the two 16S rRNA datasets indicated that an assigned species often correlated strongly with several Illumina ASVs (Figure 6). For the turbid samples RA1, Legionellaceae was not detected by the classifier EMU. Also, the classifier did not detect Legionella in sample C3 and C24. The relative abundance of ASVs assigned to Legionella species was highest at B3 (EMU: ~0.35% of total reads, kraken: ~100 reads out of 48468 and bracken: 0.45% of total reads), about 4-9 times higher relative abundance, according to the different classifiers, compared to the stations with the lowest relative abundance (RA1, C24 and C3). This is supported by qPCR analysis performed targeting a Legionella specific region of the 16S rRNA. Occurrence of Legionella specific DNA was confirmed for all samples and ranged from 2 to 40 bacteria/mL sea water (low versus high abundance). The A5, B3, RA1 and RA2 samples all had Legionella corresponding to >20 bacteria/mL sea water, while samples from C14, C24, C3 and GA1 all had Legionella corresponding to 2 bacteria/mL sea water (Supplementary Figure S4).

Figure 5 Relative abundance of nanopore Amplicon Sequence Variant reads (ASV) of genera Legionella at each station, assigned by classification algorithm EMU (top), kraken (middle) and bracken (bottom). Each sample was sequenced twice, after 30 and 33 PCR cycles.

Figure 6 Correlations (≥0.5) between ASVs in the family Legionellaceae (Illumina) and species classification in the Legionella genus (Nanopore, bracken classification).

The ASVs that were assigned to L. pneumophila by kraken was aligned along with other 16S rRNA sequences for other species of Legionella (Figure 7). Phylogenetic analysis reveals that most sequences are most similar to L. sp011764505 (deposited in NCBI as L. antarctica). Legionella sp. GA1-1 is closest related to Legionella A oakridgensis while Legionella sp. A5-1 and Legionella sp. A5-2 forms a novel branch in the tree without previous characterized species representative. None of the sequences generated in this study were most similar to L. pneumophila. Also, the Legionella specific amplification of the mip gene were negative for all samples (data not shown), indicating that the Legionella species detected in the Gulf of Bothnia samples does not have this exact mip variant that is associated with virulence, and suggest that the detected species aren’t the virulent L. pneumophila.

Figure 7 Phylogenetic relationship of 16S rRNA gene sequences of species of Legionella and ASVs assigned to L. pneumophila by kraken.

Potential hosts: Phytoplankton and phagotrophic protozoa

The estimated biomass of phagotrophic protozoa using microscopy are presented in Supplementary Figure S5. Station B3 had the largest phytoplankton biomass, over 70 µgC/L (70000 µgC/m3), while station A13, GA1, C3 and C24 had a biomass around 40 µgC/L. Station RA1 and RA2 had the lowest biomass, around 5 µgC/L (station A5 and C14 were not counted). Diatoms were abundant at the A13, B3 and RA1 stations. For the phagotrophic protists, the biomass were higher in the south, where the ciliate Mesodinium rubrum was abundant.

For the genetic method, amplification of the 18S rRNA gene (region V6-V8) indicate that the ciliate Mesodinium rubrum was most relatively abundant compared to other phagotrophic protozoans, except in the RA-stations. At the A-stations and at station B3, Woloszynskia halophila was detected, an organism that had positive co-occurrence with several Legionella ASVs (Supplementary Figures S6, S7). Unlike the 18S V9 region, this longer 18S V6-V8 region were able to assign the majority of the taxa to species level (>70%), which enabled investigation of specific species co-occurrence (Supplementary Figure S7). In general, the co-occurrence varied between Legionella ASVs and phagotrophic protozoa. Legionella ASVs number 5, 6 and 14 had positive co-occurrence with Mesodinium, Protaspa lineage, Botuliforma, and Telonema.

GLLVMs analysis identify potential drivers for the occurrence of Legionella spp.

To identify drivers for the occurrence of Legionella spp. we corrected for sampling location while taking dependencies among taxa into account (Figures 2, 6). Included in the model were ASVs that occurred in more than one sample, 14 in total (Figure 4). Significant parameter estimates for each organismal group show that humic substances had a positive effect on the ASV abundance for 6 of the Legionellaceae ASVs (Figure 8). Humic substances were the predictor with the most significant positive associations for Legionella ASVs.

Figure 8 Coefficient estimates by GLLVM analysis showing average effect and uncertainty (as 95% CI) for each organismal group.

The contribution of each independent variable in the model is found in Table 1, where temperature explain 9% of the total variation of the Legionellaceae ASVs. TDP explains 15%, coastal habitat explain 19%, bacterial production 10%, chlorophyll-a explains 11% of the total variation, salinity explain 6% and humic substances explain 3%. The full model explains 48%, including the predictors Humic substances, TDP, Coast, Temperature and Filtered volume. Salinity were not included in the model since it had a strong negative correlation to humic substances (correlation: -0.94). One ASVs abundance were negatively influenced by TDP and one were positively influenced. Four ASVs abundances were negatively influenced by temperature. Five ASV abundances are significant for a coastal habitat, some ASVs were negatively influenced by reduced filtered volumes (when < 500 ml water were filtered, Supplementary Table S3), only one ASV were significantly negatively impacted. Overall, three potential niches were identified. Legionella ASV number 4 and 13 co-occurred. Legionella ASV number 5, 6, 8 and 10 were not strongly co-occurring with other taxa. The remaining Legionella ASVs were strongly positively correlated (Supplementary Figure S8). This remaining subset of Legionella ASVs were positively correlated with humic substances (Figure 8, Supplementary Table S4).

Discussion

Data on occurrence and persistence of potentially pathogenic bacteria like Legionella spp. in natural waters are missing, which hinders assessment of the risk of exposure. In this study we identified humic substances to be a potential environmental driver for occurrence of Legionella species ASVs (6 out of 14) in the northern Baltic Sea, where the ASVs were more relatively abundant in the humic northern waters. Here, we for the first time present data of the species biodiversity of Legionella from long 16S rRNA sequencing using the Nanopore technologies. This rising technology unravelled opportunities and challenges for describing the occurrence of this public health relevant genus in natural waters.

Geographic variation of Legionella species sequences in the Gulf of Bothnia

A previous study identified Legionella in a 10 km coastal area in the Bothnian Sea during June 2018 using short 16S rRNA Illumina sequencing (Eriksson et al., 2022). Here, we identify Legionella both using standard short 16S rRNA Illumina sequencing and by using longer Nanopore sequencing. With Illumina, Legionella was mainly detected in the Bothnian Bay where Legionella was confirmed at all stations. In addition to detection in the Bothnian Bay, Legionella was detected at station B3, confirming the occurrence in this region. Interestingly, at B3 the relative abundance of Legionella ASVs was the highest of all stations, constituting ~0.15% of the reads (Supplementary Figure S2, top). Regarding the ASVs of the family Legionellaceae, these were detected in all samples (Supplementary Figure S2, bottom). Some specific ASVs were mainly detected in the offshore in the north (e.g. ASV no. 14), while others were mainly detected in the southern samples (e.g. ASV no. 6), indicating a spatial variation of Legionella spp. in the Gulf of Bothnia.

Previous studies have suggested temperature as a driver for Legionellales in freshwater (Graells et al., 2018). Generally, for the Gulf of Bothnia, the south-north differences in temperature and the salinity gradient makes the Bothnian Bay more similar to freshwater. A study of bacteria (16S rRNA amplicon sequences) along the salinity gradient of the Baltic Sea argued that the Bothnian Bay is an easier habitat to exploit by species adapted to freshwater compared to more saline water (Herlemann et al., 2011), such as Legionella. This is consistent with our finding of the highest relative abundance of Legionella spp. at the northern stations where the salinity was lower. However, in our model, the relative abundance of only one ASV (no. 5) were significantly correlated to higher temperature. Instead, the occurrence of several ASVs were significantly correlated to lower temperatures (Figure 8). Temperature was thus not found to be a driver for Legionella spp. in this study.

The microbial food web link: humic substances, iron and Legionella

In contrast to our expectation, predators (protozoa) were more abundant in the south of the Gulf of Bothnia (Supplementary Figure S5), where there were less humic substances. It is, however, notable that the total phytoplankton biomass was highest at station B3, where the relative abundance of Legionella spp. were also the highest among the samples (Figure 5, Supplementary Figure S4). Most mixotrophic algae can consume bacteria through phagocytosis (Stoecker et al., 2017). Non-phagotrophic algae such as diatoms have other endocytic mechanisms for ingesting nutrients like iron (Sutak et al., 2020). Interestingly, by taking advantage of the process of endocytosis, it has been shown that bacteria can enter eukaryotic cells (Bonazzi and Cossart, 2006). Since Legionella mainly is an intracellular organism, this observation opens for the hypothesis that Legionella also can use phytoplankton as host. However, the standard monitoring of the Baltic Sea does not include counting of strictly heterotrophic protists, which is limiting the assessment of potential hosts. DOC and predation pressure from unicellular eukaryotes, such as amoebas, together drive the evolution of bacteria. If the amounts of DOC and predators increase, this bacterial resilience to stress may also rise (Fang et al., 2022). Therefore, it can be argued that bacteria that already are resistant to predation might have an ecological advantage.

Another factor that may impact the distribution of Legionella is availability of iron, which has shown to play an important role for the growth, survival and virulence of L. pneumophila. Several studies showed that the growth of this bacterium depends on the presence of iron (Cianciotto, 2015; Portier et al., 2016). Several intracellular bacteria has mechanisms to carefully take up iron without killing their host (Ratledge and Dover, 2000; Wandersman and Delepelaire, 2004; Leon-Sicairos et al., 2015; Yan et al., 2021). Iron has been described as the most valuable metal from a metabolic perspective, since a multitude of enzymes and co-factor require this vital compound (Klebba et al., 2021). Interestingly, there is a correlation between iron sequestration and pathogenicity, and many gram-negative bacterial pathogens have elaborate pathways to acquire iron (Klebba et al., 2021).

Importantly, iron is bound to humic substances and is transported into the ocean together along with the humic substances, providing iron beyond estuaries, especially during spring when the water flow is higher (Herzog et al., 2020). When the pH is lower, iron can be released. In soils, iron has shown to bind humic acids at pH 7 twice as much compared to at pH 5 (Boguta et al., 2019). Humic substances can thus control the biological availability of iron, also in seawater (Laglera and Van Den Berg, 2009). Several studies have investigated the role of humic substances complexed with iron in seawater (Yang et al., 2017; Whitby et al., 2020), but humic substances potential role in the uptake of iron by Legionella in marginal seas is unknown. However, a new iron-regulated gene (IroT/MavN) involved in ferrous iron transport has been identified in Legionella (Cianciotto, 2015; Portier et al., 2015). In a climate change perspective, a further understanding of the microbial degradation of humic substances may be of great importance, especially in humic-rich waters like boreal marginal seas.

Along the gradients of the Gulf of Bothnia, the pH ranged from 7.25 (RA1) to 8.31 (C3). Station RA1 is located near a river mouth (Råne river), and is characterized as an estuarine area (Soerensen et al., 2017). Station B3 is also situated 5 km from a river mouth (Öre river). Of the Legionellaceae, several ASV abundances (6 out of 14) were positively linked to humic substances, perhaps as an indirect effect of iron availability. Since climate change most likely will cause increased precipitation in this region of northern Europe (Meier et al., 2022), freshwater inflows of iron-containing terrestrial humic substances can be expected to rise in the northern Baltic Sea. Potentially, this may lead to promotion of iron-demanding Legionella spp.

Standard short 16S rRNA Illumina sequencing: Potential misleading conclusions regarding the diverse genus Legionella

The study of particular bacterial species in their natural habitat is challenging due to the limitations with the standard methods used to study these organisms. Therefore, bacteria can most often only be resolved to genus level by using the standard short Illumina amplicon sequencing. For Legionella, only two ASVs were assigned to the genus, all detected in the northern stations (RA1, RA2, A13, A5 and B3, Supplementary Figure S2). These two ASVs are not representative of the whole family tree of Legionellaceae (Supplementary Figure S3). Thus, potential misleading conclusions could easily be drawn when studying this genus based solely on short amplicon Illumina data, since the majority of the diversity is not considered.

To complement the standard Illumina sequencing, we used Nanopore sequencing. The advancement of the Nanopore sequencing technologies are providing new possibilities to sequence longer amplicons (Benítez-Páez et al., 2016; Heikema et al., 2020; Rodríguez-Pérez et al., 2021), which has the potential to classify the bacteria to species level (Johnson et al., 2019; Urban et al., 2021). We performed full-length nanopore sequencing with the aim to resolve the genus of Legionella to species level. To verify the robustness of the technique, we sequenced each sample twice, with 30 and 33 PCR cycles. We used the software EMU that is specifically developed for erogenous ONT sequencing data, since it has been suggested to outcompete other software considering specificity (Curry et al., 2022). In addition, we used classifier kraken and bracken for comparison. Nanopore sequencing resulted in comparable relative abundances to the Illumina dataset for the family Legionellaceae. It therefore seems that the Nanopore sequencing results suggested species for all of the Legionellaceae reads detected using Illumina sequencing (Figure 5, Supplementary Figure S2).

Error prone Nanopore sequencing remain a challenge for accurate species classification of environmental strains

Nanopore sequencing is a promising technique for bacterial species classification. However, current limitations on base accuracy together with limitations in database completeness, limits classification specificity, especially when working on environmental samples with many uncharacterized environmental species (Latorre-Pérez et al., 2021). Therefore, caution is needed when interpreting read classification data based on ONT sequencing. Comparisons between short and long read 16S rRNA classification has shown that there are different biases in the two methods which complicates direct comparison of the results. A question that arise is how to determine if a particular classified ASV from a Nanopore sequencing run belong to an already sequenced species, or if it comes from a non-sequenced relative? For example, Legionellaceae is a diverse family and the majority of environmental variants remain uncultured (Graells et al., 2018). EMU classified all ASVs to species; the level of uncertainty allowed for classification of taxa in environmental samples could therefore be a target for future optimization.

Classification-bias towards virulent L. pneumophila

In addition to the more error prone sequences, the pathogenic strain L. pneumophila is well studied and hence overrepresented in databases. Even though L. pneumophila is the most virulent among the Legionellaceae, more than 20 other species of Legionella have shown to cause disease. The remaining strains are considered to be non-pathogens until the opposite has been proven (Legionella Species - Special Pathogens Laboratory, 2020). Thus, future efforts in the sequencing of isolates from the field and patients could aid in reducing this bias.

The Gulf of Bothnia may harbor understudied species of pathogenic Legionella spp.

For non-L. pneumophila strains, cases of disease are estimated to be underrepresented due to the standard testing methods to confirm Legionnaires’ disease and Pontiac fever (Chambers et al., 2021), and alternative methods for detection is debated, often including specific gene amplification to facilitate time consuming culture based methods (Vittal et al., 2022).

By combining short and long 16S rRNA amplicon sequencing to map the distribution of Legionella spp. in the Gulf of Bothnia, this study highlights that there is an environmental niche for several species of Legionella in the Gulf of Bothnia. This study also suggest that these public health relevant bacteria that previously has shown the capability to cause disease such as L. oakridgensis, L. parisiensis, L. longbeachae, L. cherrii and L. anisa (Bibb et al., 1981; Tang et al., 1985; Bornstein et al., 1989; Fang et al., 1989; Lo Presti et al., 1997) can be part of the natural diversity. It is known that pathogenic strains live in natural aquatic environments, and that they are being selected in manmade water systems such as drinking water pipes and cooling towers (Lesnik et al., 2015; Tsao et al., 2019). Therefore, these bacteria could be of importance considering climate change and future water management. Future field studies should thus focus on how species of Legionella are influenced by natural gradients such as humic substances, iron, temperature and salinity over time. Genus specific targets such as the non-Legionella pneumophila target rpoB gene (Pascale et al., 2021) may also be useful to further resolve the pathogenic and non-pathogenic species of Legionella spp. In addition, experimental studies could provide direct evidence of environmental drivers for the occurrence of Legionella, preferably separating the effect of humic substances and salinity.

Conclusion and outlook

In conclusion, this study highlights that there is an environmental niche for the occurrence of several species of Legionella in the Gulf of Bothnia, and that pathogenic bacteria can be part of the natural diversity. Partial 16S rRNA sequencing did not capture the diversity within Legionella while the third generation long Nanopore 16S rRNA sequencing was shown promising for bacterial species classification, even though harmonization between classification methods is still needed. Importantly, the Nanopore classification methods had a bias toward the well-studied L. pneumophila, limiting classification of environmental strains. Thus, sequencing environmental isolates could improve the classification using this powerful and rising sequencing method. These bacteria are facultative intracellular, can be pathogenic and have anti-predator strategies. Therefore, knowledge regarding potential hosts for these bacteria such as phagotrophic protozoa could be crucial, since they might provide a niche for emerging pathogens. Currently, taxonomic information of bacteria and heterotrophic protozoa is not included in the environmental monitoring. Here, we demonstrate the significance of sequencing techniques to map public health relevant bacteria such as Legionella spp. in the environment, using a method that can be used to study all bacteria. By applying generalized linear latent variable models (GLLVMs) we show that the occurrence of opportunistic Legionella spp. were linked to high concentrations of humic substances and low salinity. Considering climate change projections in this area, with increased humification and freshwater inflow, this bacterial niche might become more widespread in the future Baltic Sea. Environmental monitoring including sequencing of bacteria and protozoa could thus help identify climate change induced ecosystem imbalances and to appropriately respond to emerging diseases.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA901258.

Author contributions

KE and AA conceived the study. AA organized the data collection and KE processed the samples for sequencing. KE, JA and DS processed the data after sequencing. KE and JA conducted the statistical data analyses. KR and AS performed the phylogenetic analysis. AS and LK developed the Nanopore 16S rRNA amplicon workflow. JT organized the qPCR. KE wrote the manuscript together with all co-authors. All authors contributed to the article and approved the submitted version.

Funding

This project was financed by the Industrial doctoral school at Umeå University, the marine research environment EcoChange (FORMAS) and the Swedish Ministry of Defence (grant A4001).

Acknowledgments

We thank Umeå Marine Sciences Centre for chemical analyses, Ingrid Dacklin and Emelie Näslund Salomonsson for help with the sequencing. We thank Stina Bäckman and Mattias Jonsson for the help with the qPCR. Umeå Marine Science Centre is acknowledged for providing data on physicochemical variables and phytoplankton microscopy, via the Swedish National Marine Monitoring Program.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.1070341/full#supplementary-material

Supplementary data related to this article can be accessed at the Short Read Archive (https://www.ncbi.nlm.nih.gov/sra), project accession: PRJNA, ID: 901258. Supplementary Material 1 contains all associations between ASVs and predictors included in the GLLVM. ASV tables are provided in Supplementary Material 2 and 3.

References