A Modern Framework for Identifying Novel Environmental Legionella Species

1. Introduction

The genus Legionella, the sole genus within the family Legionellaceae, comprises more than 60 validly described species (67), with new taxa continuing to be identified as environmental surveillance and genome-based taxonomy advance [1]. Since its first recognition during the investigation of a pneumonia outbreak in 1976, Legionella has been established as an environmentally ubiquitous bacterium of major public health relevance. Members of the genus are Gram-negative, strictly aerobic rods that require L-cysteine and iron salts for growth and are primarily associated with aquatic environments.

In natural and engineered water systems, Legionella species persist and replicate predominantly within free-living amoebae, an intracellular niche that enhances their ecological fitness and tolerance to environmental stressors such as disinfectants, temperature fluctuations, and nutrient limitation. Approximately one third of described Legionella species are considered opportunistic human pathogens capable of causing legionellosis, a spectrum of respiratory diseases ranging from Pontiac fever to severe Legionnaires’ disease. Infection typically occurs through inhalation of contaminated aerosols, linking disease risk directly to the presence and proliferation of Legionella in water systems. Consequently, proactive monitoring and effective control of Legionella in engineered environments remain essential components of public health prevention strategies.

Engineered water systems, including cooling circuits, potable water installations, and other complex aquatic infrastructures, provide favourable conditions for Legionella growth through biofilm formation, stable temperatures, and stagnant water. Routine environmental surveillance programs rely largely on standardized culture-based methods to detect and enumerate Legionella; however, these approaches offer limited resolution beyond genus-level identification. Phenotypic, serological, and proteomic methods, such as latex agglutination tests and MALDI-TOF mass spectrometry, enable rapid screening but frequently fail to reliably distinguish closely related species, particularly among non-Legionella pneumophila environmental isolates.

Recent genome-based studies, including our previous work, have demonstrated that environmental Legionella isolates may represent previously unrecognized species when assessed using whole genome sequencing and genome-wide relatedness metrics [2,3,4]. Despite the increasing accessibility of genomic tools, their application in routine environmental surveillance remains largely unstructured, and clear decision-making frameworks guiding their stepwise implementation are lacking.

The aim of this study was to develop a structured framework for the identification and taxonomic assessment of novel environmental Legionella species, similar like NOVA study [5]. The framework integrates standardized culture-based detection with proteomic, molecular, and genome-resolved analyses, providing a transparent and scalable pathway for species-level identification and recognition of candidate novel taxa. By aligning advanced genomic approaches with routine surveillance practice, this approach aims to enhance both the consistency of Legionella identification and the systematic exploration of environmental Legionella diversity.

2. Materials and Methods

A framework was developed and used for the identification and characterization of potentially novel environmental Legionella species (Figure 1).

Environmental water samples were collected from engineered water systems as part of routine Legionella surveillance programme.

The study focused exclusively on environmental Legionella isolates recovered from water systems and did not include clinical samples or bacterial genera other than Legionella.

Laboratory culture of samples followed the methods according to ISO/IEC 11731:2017 (Water quality – Enumeration of Legionella), as well as the procedures described in the ‘Procedures for the Recovery of Legionella from the Environment’ (January 2005, USDHHS, Public Health Service, CDC, Atlanta, GA) [4,5]. The samples were cultivated on buffered charcoal yeast extract (BCYE) media and glycinevancomycin-polymyxin-cycloheximide (GVPC) agar and incubated at 36 °C with 3% CO₂ for 10 days. Presumptive Legionella colonies were further confirmed on culture media; BCYE agar Cys+ and BCYE Cys- agar in accordance with the ISO 11731 method in an accredited microbiological laboratory (HAA 1550).

Species identification of bacterial isolates was performed using matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF MS) (Bruker Daltonics, Bremen, Germany) with full extraction protocol as previously described [6,7]. Briefly, a loopful of a bacterial colony was suspended in 300 µl of deionized water and vortexed, after which 900 µl of absolute ethanol (Kemika, Croatia) was added. The suspension was centrifuged at 13,000 rpm for 2 min, the supernatant was decanted, and the resulting pellet was resuspended in an equal volume of 70% formic acid (Sigma Aldrich, Germany) and 100% acetonitrile (Fisher Chemical, Spain), followed by a second centrifugation at 13,000 rpm for 2 min. The supernatant was spotted onto a 96-spot polished steel target plate (Bruker Daltonik, Germany), air-dried at room temperature, and overlaid with 1 µl of a saturated solution of α-cyano-4-hydroxycinnamic acid (10 mg/ml; Bruker Daltonik, Germany) prepared in 50% acetonitrile and 2.5% trifluoroacetic acid. Mass spectra were acquired in positive linear ion mode within a mass range of 2–20 kDa. Spectral analysis and identification were performed using MBT Compass HT software version 5.0 (Bruker Daltonik) by comparison with the Bruker Biotyper reference database version 11. Species identification was interpreted based on log score values, where scores of 2.00–3.00 indicated high-confidence identification, scores of 1.70–1.99 indicated low-confidence identification, and scores below 1.70 were considered unreliable. Isolates yielding low-confidence identification scores were selected for further molecular and genomic analysis within the framework which is partial 16S rRNA gene PCR and sequence analysis.

PCR was carried out to amplify the 16S rRNA according to [10]. PCR products were visualized by gel electrophoresis on 1% agarose gel stained with Xpert Green DNA stain (GRISP Lda., Porto, Portugal), purified and sent to Macrogen Inc. for Sanger sequencing.

The obtained sequences were compared against 16S rRNA gene nucleotide databases available through the National Center for Biotechnology Information (NCBI) using the BLAST network service [11]. Sequence similarity was assessed by comparison with the closest validly published bacterial species. Bacterial species were considered validly described only if they were listed as validly published in the List of Prokaryotic Names with Standing in Nomenclature (LPSN), maintained by the German Collection of Microorganisms and Cell Cultures [1].

DNA extraction, sequencing, de novo hybrid assembly and genome annotation of Legionella strains were performed following the methods outlined in our previous studies [2,10]. Genomic DNA was extracted from a single bacterial colony grown on BYCE agar. The colony was resuspended in phosphate-buffered saline (PBS), washed three times, and DNA was isolated using the DNeasy Blood and Tissue Kit (Qiagen, Germany) with RNase treatment, following the manufacturer’s instructions. Genomic libraries for short-read sequencing were prepared using the Illumina Nextera DNA library preparation kit and sequenced on the Illumina MiSeq platform to generate 150 bp paired-end reads. Long-read sequencing was performed using the MinION Nanopore system (Oxford Nanopore Technologies, United Kingdom), with libraries prepared using the rapid 96 barcoding kit and sequenced on an R9.4.1 flow cell until sufficient data yield was obtained. Basecalling was carried out in real time using Guppy v5.0.11 with the high-accuracy model, and reads were exported in FASTQ format. Quality control and trimming of Illumina reads were performed using fastp (v0.20.1), while adapter trimming and quality assessment of Nanopore reads were conducted using Porechop (v0.2.4) and NanoPlot (v1.28.2), respectively. Hybrid de novo genome assembly was performed using Unicycler (v0.4.1) with default settings. Assembly quality was evaluated using QUAST (v5.2.0), and quality control reports were summarized with MultiQC. Functional genome annotation was performed using Prokka (Galaxy v1.14.6). Genome completeness and contamination were assessed using CheckM (v1.0.18, KBase), yielding values within acceptable thresholds for high-quality genome assemblies.

In our previous study, several genome-based typing analyses were performed to determine average nucleotide identity (ANI) [2]. Also, the calculation of genome-to-genome distance by in silico DNA– DNA hybridization was conducted against Legionella type strains. These analyses revealed four novel species candidates: PATHC032, PATHC035, PATHC038, and PATHC039. One of these, PATHC038, has already been validly described, Legionella sheltonii [3]. This study expands on those findings by applying phylogenomic analyses, including protein-based overall genomic relatedness indices (OGRIs) - average amino acid identity (AAI) and percentage of conserved proteins (POCP) to further resolve species boundaries for the remaining three isolates: PATHC032, PATHC039, and PATHC035 [12]. These metrics were applied in combination to evaluate taxonomic distinctiveness, particularly in cases where nucleotide-based thresholds alone may be inconclusive.

Whole-genome sequence data were uploaded to the Type (Strain) Genome Server (TYGS) for a whole-genome-based taxonomic analysis [13]. The analysis was performed using the current TYGS pipeline, including recently implemented methodological updates and features [14,15]. Taxonomic nomenclature, synonymy, and associated literature were obtained through the TYGS's sister database, List of Prokaryotic names with Standing in Nomenclature [14,15].

The identification of closely related type strains was carried out using two complementary approaches. First, all submitted genomes were compared against type strain genomes available in the TYGS database using the MASH algorithm, which provides a rapid estimate of intergenomic relatedness [16]. For each query genome, the ten type strains with the smallest MASH distances were selected. Second, an additional set of closely related type strains was identified based on 16S rRNA gene sequences. These sequences were extracted from the query genomes using RNAmmer [17] and compared against 16S rRNA gene sequences of all type strains in the TYGS database using BLAST [18].

The 50 best-matching type strains, ranked according to bitscore, were retained for each genome, and precise intergenomic distances were subsequently calculated using the Genome BLAST Distance Phylogeny (GBDP) approach with the coverage algorithm and distance formula d5 [19]. From these results, the ten closest type strain genomes for each of the submitted genomes were selected.

For phylogenomic inference, pairwise genome comparisons among all selected genomes were conducted using GBDP with the trimming algorithm and distance formula d5 [19]. One hundred distance replicates were calculated for each comparison. Digital DNA–DNA hybridization (dDDH) values and confidence intervals were estimated using the recommended settings of the GGDC 4.0 [14,19].

Phylogenetic relationships were inferred from the resulting intergenomic distances using a balanced minimum evolution tree generated with FASTME v2.1.6.1, including SPR post-processing [20]. Branch support was estimated from 100 pseudo-bootstrap replicates. Tree was rooted at the midpoint and visualized using PhyD3 [21,22].

Species-level clustering was performed using a type-based approach with a 70% dDDH threshold around each type strain [13], while subspecies clustering was determined using a 79% dDDH threshold [23].

The genomes and proteomes of isolates were evaluated alongside 36 Legionella type strains selected based on the availability and quality of publicly accessible genomic and proteomic data, representing the major phylogenetic diversity within the genus. For phylogenomic analyses, we ensured that the closest matches, as identified by MALDI-TOF MS, were included, enabling us to compare the proposed novel species with their closest relatives. The genome sequence data of reference type strains were obtained from the NCBI whole-genome sequencing (Table S1). AAI was calculated using the EzAAI pipeline [24]. POCP was calculated using the POCP-nf pipeline [25,26]. The de novo phylogenomic analyses were conducted with proteomic data using PhyloPhlAn 3.0 [27]. To reconstruct the phylogenomic tree, we used the supermatrix approach on 400 universal marker genes and 100 bootstrap replicates [28], with the diversity parameter set to "low" and other parameters set to their default values. The resulting tree was graphically adapted using iToL [29].

Isolates showing ANI values below 95% and dDDH values below 70%, while maintaining AAI values above 60–65% and POCP values above 50%, and forming distinct clusters, were considered candidates for novel taxa within the genus Legionella. Additionally, as a supplementary novel approach, we considered POCP values of about 90% or lower as indicative of significant species-level divergence.

3. Results

All isolates were successfully cultured under conventional aerobic conditions in accordance with ISO/IEC 11731:2017 and confirmed as members of the genus Legionella during the initial identification steps of the workflow. Preliminary identification using MALDI-TOF mass spectrometry enabled rapid genus-level assignment for all isolates.

Partial 16S rRNA gene sequencing confirmed affiliation with the genus Legionella for all isolates but did not provide sufficient resolution to reliably distinguish closely related species (Table 1). Consequently, these isolates progressed to genome-based analysis within the proposed guidelines.

Whole genome sequencing followed by comparative genomic analyses revealed clear differentiation among the four isolates. All isolates, except PATHC039, which approached the ANI threshold, exhibited levels of genomic divergence that fell outside accepted species delineation thresholds when compared with their closest validly described relatives. (Table 2).

To put the observed genomic divergence into the phylogenetic context, we conducted a whole-genome-based TYGS analysis. The results, particularly TYGS species clustering, suggest that isolates PATHC032, PATHC035, and PATHC039 likely represent novel species (Figure 2). The TYGS analysis showed that none of the isolates clustered with any currently described Legionella species. Instead, the isolates formed distinct and well-separated clusters, displaying levels of genomic divergence from their closest relatives comparable to those observed between already established Legionella species, supporting their classification as novel taxa within the genus (Figure 2, Supplementary File 1).

Protein-based OGRIs further supported these findings. AAI values ranged from 94.3% to 96.6%, closely matching ANI values, while POCP values ranged from 86.5% to 91.0%, indicating a substantial fraction of proteins lacking detectable homologs between genomes, consistent with distinct species-level divergence. Protein-based OGRIs analyses indicated that the three unresolved isolates were distinct from all currently validly described Legionella species, while maintaining relatedness consistent with placement within the genus (Table 2, Figure 3).

To contextualize the observed proteome-level divergence, we reconstructed a whole-proteome phylogenomic tree of the available Legionella type strains (Figure 4). Even at this more conservative level, the total amount of inferred evolutionary change separating the proposed novel species from their closest relatives is relatively large, comparable to some other pairs of sister species, e.g. L. rubrilucens and L. taurinensis, and L. gormanii and L. qingyii.

The overall topology of the PhyloPhlAn (proteomic) and TYGS (genomic) trees was largely congruent at the level of shallow and intermediate nodes, with both trees showing highly similar groups of species. However, some differences in deeper branching patterns were observed between the two approaches, suggesting that phylogenetic relationships among more distantly related Legionella lineages may still require further resolution through expanded taxon sampling and additional genome-based analyses. Regardless, similar compositions of shallow nodes in both analyses support the robustness of the phylogenomics approach for separating species, even if they do not recover identical phylogenies at deeper timescales.

These results, together with genome-wide comparisons, suggest that the three isolates represent separate Legionella species that have not yet been validly described.

Overall, the application of our proposed framework enabled clear separation of a known Legionella species from three environmentally derived isolates representing putative novel species, demonstrating the effectiveness of the stepwise approach for resolving taxonomic uncertainty in environmental Legionella surveillance.

4. Discussion

Routine environmental surveillance is essential for monitoring the presence of Legionella in engineered water systems; however, species-level identification remains challenging when relying exclusively on conventional laboratory methods. In this study, we present a structured framework that integrates standardized culture-based detection with genome-resolved analyses to improve taxonomic resolution of environmental Legionella isolates. Application of this workflow highlights both the limitations of routinely used identification tools and the added value of a stepwise genomic approach within environmental microbiology.

Culture-based isolation in accordance with ISO/IEC 11731 remains indispensable for the detection of Legionella in water samples and provides the foundation for downstream analyses. Consistent with previous reports, phenotypic confirmation and MALDI-TOF MS enabled rapid species-level identification of all isolates recovered in this study.

Despite yielding high-confidence MALDI-TOF MS identification scores (≥2.0), isolates were shown by genome-based analyses to be genomically distinct from their assigned reference species. This apparent discrepancy reflects intrinsic limitations of proteomic identification when applied to environmental Legionella isolates and does not represent a contradiction between methods.

MALDI-TOF MS identification relies on spectral similarity of a limited set of the dominant ribosomal protein profiles and provides the best possible match within the constraints of the reference database [29,30]. Consequently, a high-confidence score indicates strong similarity to the closest available reference spectrum rather than definitive species-level relatedness [30,31]. The limitations of MALDI-TOF MS are further exacerbated by the incomplete representation of environmental Legionella diversity in current reference databases, which are predominantly populated by clinically relevant species and type strains. Environmental isolates representing previously uncharacterized lineages are therefore forced to cluster with the most similar available reference, resulting in apparently high-confidence identifications that mask underlying genomic divergence. In contrast, whole genome–based metrics, including ANI, dDDH, AAI and POCP, evaluate relatedness across the entire genome or proteome and are internationally accepted standards for prokaryotic species delineation [12]. The consistently subthreshold values obtained for these indices provide robust evidence that the investigated isolates represent distinct genomic lineages within the genus Legionella, despite their proteomic similarity to known species. These findings highlight the necessity of genome-resolved analyses for accurate species-level assessment and support the use of MALDI-TOF MS as a rapid screening tool rather than a definitive taxonomic method in environmental Legionella surveillance.

Partial 16S rRNA gene sequencing provided additional confirmation of genus affiliation, but as expected, offered limited discriminatory power among closely related Legionella species. The high degree of sequence conservation within the genus restricts the utility of this marker for species delineation, particularly in environmental isolates that may represent previously uncharacterized taxa. Within our proposed framework, 16S rRNA gene analysis therefore functions primarily as a screening and decision-making step rather than a definitive identification tool.

Whole genome sequencing constituted the critical resolution step of the NOVA algorithm, enabling comprehensive assessment of genomic relatedness across multiple levels [5]. Nucleotide-based comparisons using ANI and dDDH demonstrated that all investigated isolates fell below established species delineation thresholds relative to their closest validly described relatives. Notably, the ANI value for isolate PATHC039 approached the ANI threshold, a scenario that complicates species assignment as genomic databases expand and phylogenetically proximate taxa are described. These findings illustrate the limitations of relying on single genomic metrics and support the use of complementary approaches for robust taxonomic interpretation.

To address this complexity, protein-based overall genome relatedness indices (OGRIs), including AAI and POCP, were applied in combination with phylogenomic analysis based on conserved marker genes. AAI values between PATHC032 and PATHC039 are close to the species threshold, suggesting they may represent closely related taxa or variants of L. pneumophila. The discrepancy between identity-based metrics (ANI, AAI) and gene-content similarity suggests that Legionella species share a conserved core gene set while maintaining species-specific adaptations, reflected in relatively low pairwise POCP values. These methods provided consistent and independent evidence for the genomic distinctiveness of isolates PATHC032, PATHC035, and PATHC039. While all isolates clearly belonged to the genus Legionella, their AAI and POCP values, together with their phylogenomic placement, indicate that isolates PATHC032, PATHC035, and PATHC039 represent distinct taxa. The concordance of nucleotide-based and protein-based OGRIs and phylogenomic results strengthens confidence in species-level inference and demonstrates the utility of integrated genomic analyses in environmental microbiology.

Importantly, our proposed framework is designed as an extension of NOVA algorithm. By defining transparent criteria for progression to genome-based analyses—such as low-confidence MALDI-TOF MS identification the workflow enables targeted use of sequencing resources while maintaining compatibility with existing monitoring programs.

The recovery of genomically distinct Legionella isolates from diverse engineered water systems and geographical settings highlights the substantial and likely underestimated diversity of environmental Legionella. Although the pathogenic potential of these putative novel species remains unknown, their identification contributes to a more comprehensive understanding of Legionella ecology and evolution in engineered environments. Systematic recognition of such taxa is an important prerequisite for future studies addressing virulence-associated traits, host interactions and public health relevance.

In summary, our proposed guidelines, similar to the conclusions of the NOVA study, provide a reproducible and scalable framework for improving the taxonomic resolution of environmental Legionella isolates. By integrating standardized culture-based methods with genome-resolved analyses, this approach supports accurate species-level placement and facilitates the systematic detection of novel Legionella taxa. As genomic technologies become increasingly accessible, structured workflows such as the one described here may play an important role in advancing environmental microbiology and water safety research.

5. Conclusions

The proposed framework bridges routine environmental monitoring and high-resolution microbial systematics. By combining standardized isolation methods with genome-based taxonomy, it enables reliable detection, accurate taxonomic placement, and identification of novel Legionella species. This framework may provide a valuable basis for future research into the virulence potential and pathogenic characteristics of environmental Legionella isolates.