Introduction

Rapid advancements in sequencing technologies have revolutionized our view on microbial communities. While amplicon sequencing provides information on community composition and diversity, shotgun and whole genome sequencing allow us to reliably anticipate evolutionary and ecological relationships between microbes and to obtain functional information on communities. Computational models assessing the metabolic capacity of individual members, or an entire consortium, have become very popular and powerful^1–3. The major focus of such modelling approaches is typically on the primary metabolism of bacteria, as genes involved in core metabolic pathways are highly conserved and can be identified with relative ease^2,4. Conversely, analysis of the secondary metabolites has attracted less attention, even though they include compounds such as antibiotics, toxins, siderophores, biosurfactants, all known to have important implications for microbial community assembly^5,6 and to be important sources for pharmaceutical discoveries^7–9.

There are multiple challenges that currently prevent a detailed unravelling of secondary metabolism of bacteria based on genome data^10–12. First, most secondary metabolites are produced by pathways comprised of modular enzymes such as non-ribosomal peptide synthetases (NRPSs) or polyketide synthases (PKS)^13,14. Locating complete synthesis clusters and identifying all enzyme-encoding genes is challenging from highly fragmented metagenomic sequences or draft genomes with a high number of contigs. Second, functional predictions for coding regions within a cluster rely on homologous comparisons with experimentally characterized genes. Such information is often restricted to a limited number of model organisms, meaning that only a small portion of the existing secondary metabolism pathways is covered by current data bases. Finally, given the complex multi-modular synthesis machineries, it is challenging to precisely predict the secondary metabolites produced even with accurately annotated NRPS or PKS clusters. The main challenge is that a large pool of non-proteinogenic amino acids is used as substrates and the specificity of an enzyme’s A domain, connecting these unusual amino acids, is often poorly understood¹⁵. As a result, new computational methods are needed to accurately reconstruct bacterial secondary metabolism from sequence data.

Here, we present a new bioinformatic pipeline that overcomes these challenges. We specifically focus on a particular class of secondary metabolites (iron-scavenging siderophores) as a case study to develop a bioinformatic workflow that predicts the chemical structure of the produced metabolites with nearly 100% accuracy. Our pipeline is based on improved gene annotation combined with a phylogeny- and feature sequence-based substrate prediction techniques (Figure 1). In comparison with the currently available databases and bioinformatic tools¹⁶, the main advancement of our method is the more accurate prediction of synthesized products based on NRPS clusters identified in genome data. In addition, we show that our workflow is capable of correcting previous mis-annotations and is extendable to other secondary metabolites (e.g., toxins, biosurfactants, virulence factors) and microbial genera.

Figure 1

Among siderophores, we focus on NRPS machineries that are responsible for the synthesis of pyoverdines, a class of chemically diverse siderophores with high iron affinity, produced by Pseudomonas bacteria^17,18. While each Pseudomonas strain produces a single type of pyoverdine, an enormous structural diversity has been described across strains and species^19–22. Pyoverdine types differ in their peptide backbone, meaning that the diversity should be mirrored in NRPS enzyme diversity and their selectivity for the different amino acid substrates¹⁸. Based on this knowledge, our pipeline entails the following steps (Figure 1): (i) identification of the complete sequences of pyoverdine synthetase genes from fragmented draft genomes, (ii) building the pyoverdine synthesis machinery in silico by extracting the feature sequences for substrate specificity from motif-standardized NRPSs, and (iii) predicting the precise chemical structure of pyoverdines followed by empirical verification.

An additional element of iron metabolism is that when siderophores are secreted and bound to iron, bacteria rely on a specific receptor for their uptake into the cell. Pyoverdine receptors are annotated as FpvA and it is known that receptor diversity matches pyoverdine diversity^22,24. Moreover, FpvA belongs to the family of TonB-dependent receptors and a single Pseudomonas species often has many gene copies encoding these receptors. This poses an additional bioinformatic challenge: how to find the gene encoding the specific pyoverdine receptor among several potential receptor genes? To overcome this, we develop an algorithm that focuses on sequence regions involved in pyoverdine recognition and translocation across the outer membrane with supervised learning methods that locate the fpvA genes in the fragmented genomes based on these regions (Figure 1). Not only did our method unveil a substantial number of previously unrecognized pyoverdine receptors, but it also revealed that receptors with the potential to uptake pyoverdine are distributed across a broad spectrum of bacteria. This suggests the possibility of pyoverdine serving as another ‘shared language’ among microbes. Altogether, our bioinformatic pipeline uses knowledge-guided insights empowered by supervised learning to construct a first systematic sequence-to-function mapping of a family of secondary metabolites (pyoverdine) and their corresponding receptors. Our analysis unveils a yet unrecognized extraordinary diversity of iron-scavenging machineries in pseudomonads.

Results

Section 1 Improved annotation pipeline reveals a vast reservoir of pyoverdine synthetase genes

The first step of our bioinformatic pipeline was to improve the annotation of pyoverdine synthetase genes. The pyoverdine molecules is composed of a conserved fluorescent chromophore (Flu) and a peptide chain (Pep), which are both synthesized by NRPS enzymes²⁵. There are already existing tools, such as antiSMASH that can find and annotate NRPS clusters in microbial genomes²³. However, antiSMASH (and other popular annotation platforms^26,27) rely on accurate gene predictions, which are typically problematic for fragmented genomes. Consequently, while antiSMASH can recognize and annotate certain genes of an NRPS cluster, the precise reconstruction of a complete NRPS assembly line often fails. This is particularly problematic because most available genomes are drafts and any analysis suffers from the unavoidable issue of incomplete or misannotation of gene fragments.

To overcome these issues, we developed an improved four-step annotation pipeline starting with the raw annotation of the pyoverdine cluster obtained from antiSMASH²³ (Figure 2a). First, we implemented a NRPS Hidden Markov Model (HMM) to re-annotate and extract the entire nucleotide sequence of the pyoverdine synthetase cluster²⁸, including the genes missed by antiSMASH. For this step, the nucleotide sequences were converted into amino acid sequences to avoid erroneous gene predictions typically associated with antiSMASH. Second, we assembled the entire re-annotated pyoverdine coding region into a single sequence with a defined start (CAL) and/or end (TE) markers. Third, we used NRPSMotifFinder to identify the C, A, T, E and TE motifs that are characteristic for the NRPS structure of pyoverdine¹⁵. Finally, we applied a safety measure to ensure that the recovered NRPS assembly line is complete (contains both Flu and Pep) and is not truncated, which can occur when a synthetase coding region is at the edge of a contig. Consequently, we dismissed all pyoverdine synthesis clusters located within 100 bp proximity to contigs’ edges and either lacked Flu or Pep synthetic genes.

Next, we applied our improved pyoverdine synthesis annotation pipeline to 9599 Pseudomonas genomes (including 613 complete and 8986 draft genomes) retrieved from the Pseudomonas Genome Database²⁹. We found the pyoverdine synthesis machinery in 97% of the genomes (Figure 2b), indicating that the machinery is ubiquitous in Pseudomonas. However, since 94% of the analyzed genomes were in draft form, the pyoverdine synthesis machinery was likely truncated (i.e., on the edge of the contig) in 63.4% (6087) of the genomes. These genomes were excluded from further analysis. Around 3.1% of retained genomes (293) with high assembly completeness were missing pyoverdine synthetic genes, indicating that these Pseudomonas strains were not able to produce pyoverdine (‘non-producers’). The rest of the genomes (33.5%; 3219 genomes) were classified as ‘producers’ with complete pyoverdine NRPS assembly lines that meet all our quality controls. For these 3219 genomes, we used NRPSMotifFinder to find boundaries between the various synthesis domains and to determine amino acid length and the number of A domains. The lengths of pyoverdine synthetic genes ranged between 7690 and 21333 amino acids, and the number of A domains per synthetase ranged between 6 and 17, with a total of 35,281 A domains being present across all strains (Figure 2c and Figure S1). Overall, our analysis pipeline unveiled a vast diversity of pyoverdine synthetase that goes far beyond of what has previously been described in the literature.

Figure 2

Finally, we conducted a phylogenetic analysis based on 400 conserved genes with the 293 non-producers and the 3219 producers. We first removed redundant non-producers by retaining the most integrative genome among strains with high phylogenic similarity. Then, we removed redundant producers by retaining the most integrative genome among strains with high phylogenic and pyoverdine synthetase similarity (Figure 2d). This data cleaning yielded a total of 1928 Pseudomonas strains (403 complete and 1525 incomplete genomes), segregating into 1664 pyoverdine producers and 264 non-producers. The phylogenetic tree revealed that all major Pseudomonas species clades were present in our data set (Figure 2e). Moreover, the number of A domains varied widely among species and even between strains within species. For example, the number of Pseudomonas aeruginosa A domains ranges between 7 and 14. In summary, by improving the synthetase annotation method, we successfully obtained 1664 highly reliable pyoverdine synthetases (with a total of 18,292 A domains) and 264 non-producers.

Section 2 Phylogeny-focused substrate prediction for pyoverdine A domains

Our next goal was to precisely predict the molecular structure of the pyoverdines produced by the 1664 strains with complete synthetase gene clusters. The first essential step towards this goal was to reliably predict the substrate selectivity of all A domains in the NRPS assembly line. The A domain of each module selects for a single substrate among 22 proteinogenic and hundreds of non-proteinogenic amino acids^30,31. Moreover, whenever an E domain exists downstream of an A domain, the chirality of the amino acid incorporated into the peptide chain gets modified from L to D. Thus, the modularity combined with the selectivity of A domains can promote an enormous diversity of pyoverdine molecule structures. To date, 73 pyoverdine structures have been reported (Supplementary_table1) out of which 13 have their synthetase genes sequenced (Supplementary_table2). In order to make reliable predictions, two challenges must be addressed: (i) the extraction of relevant information from A domain sequences for which the substrate is known, and (ii) the effective application of this information to predict specificity of A domain sequences for which the substrate is unknown.

To address the first challenge, we built our analysis on the NRPS assembly lines of the known 13 pyoverdines to extract relevant information from the A domain sequences. From this dataset, we could identify 101 A domains that could be experimentally linked to 13 amino acid substrates (Supplementary_table3). We next performed multisequence alignment of the 101 A domains to determine the “feature sequence distance”, which is the most informative for the substrate selectivity. To this end, we tested three different A domain regions, three different sequence similarity measurements and seven different clustering methods for their predictive power (Figure 3a). We found that the full A domain sequence is not informative for substrate prediction (Figure 3b, left panel). Instead, our analysis indicated that information-rich positions start with motif A4 and end before motif A5, consistent with the known role of the A domain pocket in substrate selectivity¹⁵. Overall, the sequence region from motifs A4 to A5 (termed “Amotif4-5”), in conjunction with Jukes-Cantor distance and Ward linkage clustering, performed best in accurately distinguishing between different substrates and maintaining homogeneity for identical substrates (Figure 3b). Utilizing these parameters, we conducted a comprehensive analysis of all A domains and their respective substrate information within the Pseudomonas NRPS biosynthetic gene cluster available in the MIBIG database³² (Supplementary_table4). To our delight, our method can find and correct the incorrect substrate information corresponding to A domain in the MIBIG database (Figure S2 and Supplementary_table5). This demonstrates the effectiveness of performing substrate predictions for the A domain based on the selected parameter combination.

Figure 3

To address the second challenge, we developed a “phylogeny-focused method” to apply the feature sequence distance derived in the preceding paragraph to the 18,292 discovered A domains. We realized that a direct construction of a phylogenetic tree including all 18,292 query A domains and the 101 reference A domains would be computationally too demanding and impossible to scale up. Furthermore, such an approach would result in phylogeny-interference issues, where domains would cluster not only based on their substrate similarities but also based on overall species relatedness¹⁵. To minimize the effect of phylogeny and speed up calculation, we took each of the 18,292 query A domains and identified the two most similar A domain clusters within the 101 reference A domain set. We then compared the feature distance between each query A domain and the two most similar reference A domains in different clusters and assigned the query A domain to a substrate specificity using the following rules (Figure 3c). (1) If the feature distance is below the 0.7 threshold (corresponding to 50% identity) for only one of the two reference A domains, then the substrate of the query A domain is matched to the substrate of the more similar (lower distance) reference A domain. (2a) If the feature distance is below the 0.7 threshold for both reference A domains, then we considered the relative difference of the query A domain towards the two reference A domains. If the relative difference is larger than 0.2, the query A domain is matched to the substrate of the more similar reference A domain. (2b) If the relative difference is smaller than 0.2, the substrate of the query A domain cannot unambiguously be determined and is thus matched with both reference substrates. (3) If the feature distance is above the 0.7 threshold (below 50% identity) for both reference A domains, then the substrate of the query A domain is marked as “unknown”. For most query A domains, rule (1) could be applied (17880 cases), whereas rules (2) and (3) had to be used rarely (133 and 279 cases, respectively). We applied our methodology termed “phylogeny-focused method” to all following substrate and pyoverdine structure predictions.

Section 3 Experimental validation of the annotation and prediction pipeline

We tested whether our bioinformatic pipeline can reliably predict the structure of a set of yet uncharacterized pyoverdines. To achieve this objective, we selected 20 Pseudomonas strains, all known to produce pyoverdines, from a natural strain collection that was previously isolated from soil and water³³. We sequenced their genomes and subsequently applied our annotation and prediction pipeline to generate predicted pyoverdine structures for all 20 strains harboring a total of 237 A domains. Notably, none of the predicted structures matched any of the 13 reference pyoverdine structures, and 9 out of the 20 structures were predicted to be novel (not yet characterized) pyoverdines. To verify our predictions, we elucidated the chemical structure of the 20 pyoverdines using culture-based methods combined with UHPLC-HR-MS/MS³⁴.

We found a near-perfect match (96.2%) between the predicted and the observed pyoverdine chemical structures and were able to accurately assign amino acids in 228 out of 237 cases (Figure 3d). Our method demonstrated a substantial improvement comparing to the prediction accuracy of antiSMASH¹⁶ in pyoverdines (46.0%), which could accurately assign correct amino acids only in 109 out of 237 cases (Supplementary_table6). The nine non-matching cases in our analysis segregated into three groups. In three cases (1.3%), our algorithm could not distinguish between the substrates Lysine and Ornithine, as these two amino acids are highly similar both in terms of their chemical structures and corresponding A domain sequences. This is the only sensitivity issue that is associated with our approach. In four cases (1.7%), our technique assigned an “unknown” substrate to amino acids that turned out to be valine, citrulline and histidine. Indeed, these three amino acids have not been reported in pyoverdines before and are therefore not yet present in the reference dataset. These cases show that our analysis pipeline can be used to identify new substrates. Once experimentally verified, the new A domains and their substrates can expand the reference dataset, allowing targeted improvement of our phylogeny-focused prediction technique. Finally, there were only two cases (0.8%) that represented true mismatches between observed and predicted amino acids. Altogether, our phylogeny-focused method is highly accurate in predicting pyoverdine peptide structures and in identifying unknown substrates in Pseudomonas.

To assess the extendibility of our pipeline, we chose Burkholderiales as a test case, leveraging a manually curated dataset that we have accumulated from literatures (see Method for detail). This dataset contains 203 A domain sequences and corresponding substrate information for 34 NRPS metabolites (e.g., toxins, biosurfactants, virulence factors) across 7 genera of Burkholderiales (Supplementary_table7). We also integrated experimentally validated biosynthetic gene cluster data from Burkholderiales in MIBIG 3.0³² (Supplementary_table8) with our manual dataset. Upon analyzing these datasets, we observed that our pipeline consistently maintains high prediction accuracy within Burkholderiales (83% vs. 67% in antiSMASH¹⁶, Figure S3 and Supplementary_table9). This outcome underscores the robustness of our developed prediction algorithm, highlighting its potential extension to other NRPS secondary metabolites and microbial genera.

Section 4 Application of the annotation and prediction pipelines to a full dataset

After successful validation, we applied our bioinformatic pipeline to the 1664 complete NRPS assembly lines annotated in our genome analysis (Figure 2). Across all assembly lines, we were able to predict the substrates of 17,880 A domains (97.75%) without ambiguity, whereas 133 A domains (0.73%) were associated with two different substrates, and 279 A domains (1.52%) predicted an unknown substrate (similar to the case of valine above). After considering the presence/absence of an E domain in each module, we derived the structure of 1664 pyoverdines according to method at section 2 (Figure 4).

Figure 4

Our prediction yielded 188 different pyoverdine molecules, out of which only 37 structures had been previously reported. While these 37 reported structures were highly abundant across strains (1103 out of 1664), our pipeline was powerful in identifying many of the rarer pyoverdine variants. Agreeing with previous studies, we observed that the fluorophore is highly conserved among the 188 predicted structures. Moreover, our analysis confirmed that 13 amino acid substrates form the core of all the 188 pyoverdine structures, with most of the variation being attributable to different substrate combinations, peptide lengths, and substrate chirality (Figure 4). Notably, pyoverdine structural diversity was not strongly linked to phylogeny because the same pyoverdine structure could be found in completely unrelated species, while closely related species often had different pyoverdine structures (Figure 4). These observations suggest that there may be both frequent recombination and horizontal gene transfer of pyoverdine synthetase clusters between species. Taken together, the bioinformatics methods developed in our study can predict a suit of secondary metabolites (pyoverdines) from sequence data with high accuracy, revealing an unprecedented richness and evolutionary history of siderophores within pseudomonads and the discovery of 151 novel pyoverdine candidate variants.

Section 5 Development of a region-based identification method for annotation of the FpvA receptors

In pseudomonads, iron-loaded pyoverdines are recognized by FpvA, a TonB-dependent receptor, that transports the ferri-siderophore into the periplasm^19,35–37. The protein structure of characterized FpvA variants consists of three domains: The Secretin and TonB N-terminus short domain (STN), the Plug domain (Plug), and the TonB dependent receptor domain (TonB)¹⁹. While these domains are conserved across FpvA variants and other siderophore receptors, there is substantial variation at the sequence level. This makes it challenging to reliably identify FpvA receptors from sequence data by homologous search. As an example, we were unable to find FpvA genes (with a 60% identity threshold) by homologous search in several genomes although they had complete pyoverdine synthesis machineries. Moreover, there are many other TonB-dependent receptors with fairly high sequence identity to FpvA but that transport other siderophores than pyoverdine (e.g. FpvB, 55% identity, transporting pyoverdine, ferrichrome and ferrioxamine B³⁸. Therefore, it is imperative to develop a new comprehensive method for identifying FpvA receptors in Pseudomonas genomes.

We started our approach by comparing the sequences of 35 reported siderophore receptors, including 21 FpvA, 6 FpvB, and 8 TonB-dependent siderophore receptor sequences often found in Pseudomonas genomes, encoding receptors for the uptake of heterologous siderophores (Supplementary_table10). We found that all receptor sequences share a similar length of around 800 amino acids (FpvA and FpvB sequences: 809 ± 10 amino acids). We then used the complete sequences to calculate the pair-wise distances by global alignment before applying hierarchical clustering (Figure 5a). We found substantial divergence between FpvA variants to an extent that was comparable to the distance between FpvA and other siderophore receptors. Moreover, FpvB variants clustered with FpvA variants, showing that FpvA identification based on full sequence distances is unachievable. We hence focused on the three typical receptor domains (TonB, Plug, and STN, retrieved from the Pfam database) and applied Profile Hidden Markov Models (pHMM) to calculate the pHMM probability scores for each domain and reference sequence. The probability scores (calculated as the log-odd ratios for emission probabilities and log probabilities for state transitions) had reasonably high scores but no distinction was apparent between the three receptor classes (Figure 5b).

Figure 5

We next asked whether there are specific regions within the receptor sequences that are characteristic of FpvA. To address this, we conducted a multiple sequence alignment (MSA) with all 35 reference receptor sequences and mapped them onto the sequence of the well-characterized FpvA of P. aeruginosa PAO1 (Figure 5c). MSA allows to identify conserved sites (Figure 5c, black dots representing the top 10% most conserved sites) that are shared by the majority of the reference sequences. We then used these conserved sites to partition the MSA into variable regions which were flanked by two conserved sites. For each variable region, we assessed its predictive power to differentiate FpvA from non-FpvA sequences. For this we defined the “FpvA identification score” analogous to the intercluster-vs-intracluster Calinski-Harabasz variance ratio, as

where d_A:A is the sequence distance among all 21 FpvA sequences, and d_A:non is the sequence distance between all 21 FpvA and the 14 non-FpvA sequences.

Our analysis yielded two locations with noticeably high FpvA identification scores (Figure 5c). The region with the highest FpvA identification score (referred to as R1) locates at the intersection of the Plug domain and the barrel structure of the TonB domain (Figure 5d, between 258 Gly and 309 Gly in the PAO1 FpvA). According to the sequence distance matrix, the R1 region allows to distinguish heterologous siderophore receptors from FpvA and FpvB receptors (Figure 5e). The region with the second highest FpvA identification score (referred to as R2) was located in the C-terminal signaling domain (Figure 5d, between 59 Leu and 86 Lys in the PAO1 FpvA). The sequence distance matrix revealed that R2 allows to distinguish FpvB from FpvA receptors (Figure 5f).

We then constructed two pHMM by (i) the alignment of the 21 FpvA sequences in the R1 region, termed R1(FpvA), and (ii) the alignment of the 6 FpvB sequences in the R2 region, termed R2(FpvB). Running R1(FpvA) and R2(FpvB) against all 35 reference sequences revealed a clear separation between the three receptor categories (Figure 5g). Along the R1(FpvA) axis, FpvA and FpvB reference sequences have high R1 scores (minimal score 77.0) that separate them from other siderophore receptors (maximal score 38.1), whereas FpvAs references have substantially lower R2 scores (maximal 20.2) than FpvBs (minimal 49.0) along the R2(FpvB) axis.

Based on these insights, we developed a decision flow chart for annotating FpvAs in Pseudomonas genomes (Figure 5h): First, we considered sequences as Fpv-like receptors that share similar properties to the ones identified in our reference database. Particularly, protein coding sequence (CDS) length has to be between 750 and 850 amino acids and the phMM scores for the three typical receptor domains STN, Plug, and TonB have to be greater than 25, 50, and 80, respectively (Figure 5b, red dashed lines). Second, we used the pHMM threshold scores obtained for R1(FpvA) and R2(FpvB) (Figure 5g) to differentiate other siderophore receptors (R1(FpvA) score < 50) from FpvB receptors (R1(FpvA) score > 50 and R2(FpvB) score > 30) and FpvA receptors (R1(FpvA) score > 50 and R2(FpvB) score < 30). Our method effectively identifies FpvA receptors from sequence data and can be readily applied to the entire Pseudomonas dataset.

Section 6 Application of the receptor annotation pipeline to the full dataset

The region-based receptor identification pipeline was applied to all 1928 Pseudomonas genomes. The analysis identified 4547 FpvAs, 615 FpvBs, and 9139 other TonB-dependent Fpv-like receptors across the dataset (Figure 6a). The 4547 FpvA sequences clustered hierarchically into 114 groups, defined by an identity threshold of 60%. When comparing to the 21 reference FpvAs (Figure 6b), we found that 2293 FpvA sequences have close homologues in the reference data base, while 2254 FpvA sequences lack such close homologues (sequence identity < 50%). These latter sequences, termed as “dissimilar to reference”, may represent novel subtypes of FpvA receptors that could not be found by simple homology search. Our analysis further shows that many strains have more than one FpvA receptor.

Figure 6

We then asked whether the 4547 FpvAs are found in proximity of pyoverdine Pep synthetase genes as it is commonly the case for cognate FpvA receptors ³⁹. We thus calculated the proximity between pyoverdine Pep genes and the Fpv-like receptor genes by counting the number of base pairs between the two coding regions. All TonB-dependent receptors that have not been classified as FpvAs were more than 20 kb away from the Pep genes (Figure 6c). In contrast, 92% of the nearest FpvA genes were indeed located within 20 kb of their pyoverdine Pep genes (Figure 6d, called proximate receptors). These proximate receptors encompassed both those with close (66%) and more distant (34%) resemblance to the reference receptor types. Overall, this proximity analysis confirmed that our region-based gene identification method can reliably identify FpvA receptors.

We next explored the diversity among FpvA receptors in more detail by focusing on the 1534 strains that had proximate-receptors within 20 kb of the pyoverdine Pep genes (Figure 6d) and using high-confidence FpvAs for sequence feature extraction. When considering the whole gene sequences, these receptors segregated into 44 groups according to single-linkage clustering with an identity threshold of 60% (Figure S4a). To investigate which sequence regions were the most informative for reliable clustering, we used a similar approach as with FpvAs detection by quantifying the “group identification score” for variable regions flanked by highly conserved sites. The higher the score, the stronger a region’s capacity to discriminate between FpvA groups. We found that the four regions with the top discrimination capacities all located near the Plug domain surrounding the pyoverdine transmission channel (Figure 6e). The plug domain is known to undergo conformational changes and is involved in pyoverdine selectivity and import^40,41, suggesting that the four high-score regions are responsible for pyoverdine specificity.

Based on the above insights, we concatenated the four high-score regions (from 168 Pro to 295 Ala in PAO1) into a single “feature sequence”. The feature sequence could characterize 98% of the distance matrix compared to the whole sequence (1534 FpvAs, r=0.98, Figure S4a-b) and substantially reduced within-group distance. We applied the concatenated feature sequence approach to all the 4547 annotated FpvAs to calculate the sequence distance matrix. Single-linkage clustering with an identity threshold of 70% revealed a total of 94 groups, out of which 43 groups contained more than 10 members (Figure 6f). The diversity of receptors is hence much larger than currently anticipated as only 3 groups of FpvAs have previously been reported. Finally, we calculated the diversity of receptor FpvAs for each of the 13 phylogenetic clades with more than one strain by the Shannon entropy, which is similar to the alpha-diversity in microbial community (Figure 6g). We noticed large differences in FpvA diversity across the clades and species: clades with P. aeruginosa and P. syringae species had lower FpvAs diversity (1.55 and 1.60) than clades containing P. putida and P. fluorescens species (4.82 and 3.77). Taken together, the region-based identification method developed in our study can reliably mine the FpvAs (pyoverdines receptors) from genome data, revealing undiscovered diversity of FpvA pyoverdine receptors that are unequally distributed across the different phylogenetic clades of pseudomonads.

We further extended the search into the entire bacterial world, by applying the same pipeline in Figure 5h. Among the >30,000 complete bacterial genomes in NCBI, 31,936 receptors pass the R1 and R2 thresholds (Figure S5). To our surprise, these FpvAs are dispersed across 12,944 strains spanning 13 phyla, 22 classes, 54 orders, 139 families, and 468 genera, encompassing 2,598 species (Figure S6 and Supplementary_table11). While Pseudomonas remains the genus with the highest distribution of FpvA, other genera with notable distributions comprise Escherichia, Salmonella, Enterobacter, Acinetobacter, Xanthomonas, Bordetella, and Achromobacter. This study marks the first systematic exploration of FpvA’s distribution throughout the bacterial kingdom. The extensive distribution may implie that pyoverdine serves as another ‘shared language’ in the bacterial community and plays a more important role than previously thought.

Discussion

The rapid expansion of sequencing data offers exciting opportunities for microbiology^42–44. One key challenge of current research in the field is to infer biological functions of microbial communities from genome sequence data^45–47. While this endeavor is increasingly successful for biological functions involving the primary metabolism and the associated complex metabolic flux, reconstructing aspects of the secondary metabolism is much more challenging. The main issue is that neither the function of a secondary metabolite enzyme nor the resulting metabolite can be precisely predicted from gene sequence data. In our study, we tackled this challenge and developed a bioinformatic pipeline to reconstruct the complete secondary metabolism pathway of pyoverdines, a class of iron-scavenging siderophores produced by Pseudomonas spp. These secondary metabolites are synthetized by a series of non-ribosomal peptide synthetases and require a specific receptor (FpvA) for uptake. We combined knowledge-guided learning with phylogeny-based methods to predict with high accuracy: (i) the full pyoverdine assembly line, (ii) the substrate specificity for each enzyme within the assembly lines, (iii) the complete chemical structure of pyoverdines, and (iv) the FpvA receptors from genome sequences. After validation, we tested our pipeline with sequence data from 1664 phylogenetically distinct Pseudomonas strains and were able to determine 18,292 enzymatic A domains involved in pyoverdine synthesis, reliably predicted 97.8% of their substrates, identified 188 different pyoverdine molecule structures and 4547 FpvA receptor variants belonging to 94 distinct groups. The uncovered diversity is stunning and goes far beyond currently known levels of variation (73 pyoverdines and 3 FpvA groups). The molecular diversity of iron scavenging capacity highlights its importance among pseudomonads.

We show that knowledge-guided learning is an extremely powerful tool to predict enzyme, metabolite, and receptor properties. The establishment of our entire pipeline is based on only 101 previously known enzymatic A domains (from 13 known pyoverdine assembly lines) and 21 FpvA receptor sequences. Even with this limited amount of information, we were able to predict the substrates of almost all the 18,292 enzymatic A domains and to identify 4547 FpvA receptors from the sequence data. A key insight from our knowledge-guided learning is that comparisons based on the full gene sequences (e.g., for pyoverdine synthetase or receptor) are likely non-informative and unsuitable for obtaining functional information. This is because overall diversity does not stand for functional diversity, meaning that A domains recognizing the same substrate can diverge substantially in their full sequences. The same holds true for receptor sequences: whole-sequence alignments can neither accurately identify FpvA receptors nor reliably separate them into functional groups. Instead, it is imperative to extract informative feature sequences that are defined as sequence stretches within a gene whose diversity is tightly linked to variation in its functioning. We successfully extracted and applied feature sequence comparisons for both A domain substrate prediction and FpvA identification. It is important to note that a knowledge-guided pipeline does not have to be perfect right from the start. For example, our pipeline for pyoverdine structure prediction returned unknowns for several amino acid positions within the PEP. Our experimental verifications then revealed indeed new substrates such as valine and citrulline. This information can then be used to refine our prediction algorithm in a feedback loop.

Another main advantage of our bioinformatic pipeline is that it can be applied to draft genomes. This reflects a major improvement compared to existing annotation tools such as antiSMASH²³, which typically has difficulties in recognizing NRPS structures in fragmented genome assemblies. However, draft genomes are the most common data source in microbiology. While our pipeline shows high performance, we need to acknowledge that we still lose many genomes (6087 out of 9599, 63.4%). The reason for the loss is that the pyoverdine synthesis machinery is large, which increases the probability that it is positioned at the end of a contig. We decided to exclude those cases because the annotated synthesis machinery might be truncated and thus incomplete. Thus, the high loss rate of draft genomes is rather due to limitations in sequence quality (too many short contigs) and not due to a limitation of our bioinformatic pipeline. We believe that this limitation may disappear in the future as long-read sequencing technologies are quickly becoming cheaper and more reliable.

We further show that knowledge-guided learning combined with a phylogeny-focused approach is a powerful tool for predicting the substrate specificity of A domains of synthetases. It outperforms currently known bioinformatics prediction tools of NRPS substrates such as antiSMASH²³. Most current algorithms^48–52 perform poorly when applied to pyoverdines, particularly when encountering non-proteogenic amino acids. The high accuracy of our algorithm can largely be attributed to our reference set, composing only 13 pyoverdines from Pseudomonas spp., yet capturing most of the substrate diversity. Similarly accurate predictions based on a handful of known substrates among closely related species were observed in several fungal NRPS systems⁵³. It is worth noting that when the algorithm output is “unknown,” it actually signifies uncharacterized A domains not yet incorporated into the reference data set. This should prompt researchers to pay attention to these A domains, and like in our case, subject them to further experimental investigation. This approach helped us discover new substrates (valine, histidine, citrulline), which had not been previously documented in pyoverdines and were therefore absent from the reference A domains. The novel substrates identified through our structural assessment and mass spectrometry experiment can subsequently be used to enhance the precision of our phylogeny-centered substrate prediction technique in the future, creating a progressive feedback loop of expanding knowledge. Taken together, supervised learning based on a few known compounds produced by species from the same genus probably outperforms generalized prediction algorithms trained on many products from a diverse set of microbes for NRPS substrate predictions.

Our results show that both pyoverdine and receptor diversity has been vastly underestimated. While considerable pyoverdine diversity (n=73) has already been captured in previous studies, here we discovered 151 new variants. On the receptor side, the uncovered novel diversity is more dramatic. One reason for this is that research on receptors has mainly focused on the pathogen P. aeruginosa^19–22,54. For this species, three different pyoverdine types were described²¹ together with three structurally different FpvA receptor types that each recognize one of the pyoverdine types²². While our study confirmed that P. aeruginosa strains (n = 554) indeed have only 3 pyoverdine-receptor systems, we also discovered 91 new FpvA groups among environmental Pseudomonas spp. Our findings raise the question why there are so many different pyoverdine and receptor variants. One potential explanation is that the benefit of specific siderophores could be context-dependent and locally adapted to multitude of different environmental conditions pseudomonads are exposed to. For example, experimental work has revealed that pyoverdines can be cooperatively shared among strains with matching receptors^36,55, or conversely, pyoverdines can serve as competitive agents by locking away iron from species that have non-matching receptors⁵⁶. Given that bioavailable iron is limited in most natural and host-associated habitats^57–59, the unraveled functional diversity is likely a direct evolutionary consequence of the struggle and competition of microbes for iron. While experimental work is often restricted to a low number of strains, we propose that our bioinformatic pipeline can be used to predict pyoverdine-mediated interaction networks across thousands of strains and across different habitats. We will address this point in a future study.

We believe our pipeline could be easily expanded to study iron competition in multi-species communities in the future and perhaps in plant-microbe ecosystems, as siderophores exist ubiquitously and are shared among microbes⁶⁰. To move further, a key question is whether our knowledge-guided approach can be applied to other microbes and important secondary metabolites, such as antibiotics, toxins, biosurfactants and pigments? Our preliminary investigation, focused on the curated dataset in Burkholderiales as a test case, provides a highly promising response to this query. The outcomes from this evaluation predominantly showcased the method’s potential for extension to diverse microbial genera and metabolites, given sufficiently accurate dataset. Meanwhile, we observed a notable decline in accuracy when compared to pyoverdine. This indicates that our algorithm possesses potential for refinement, particularly concerning specific secondary metabolites or microbial genera. As soon as sufficient case-by-case knowledge on a specific system is available, the annotation strategies together with the feature sequence extraction and the phylogeny-focused approach developed in our paper can be applied. For most of the secondary metabolites listed above, there are no receptors as the compounds have purely extra-cellular functions, which substantially simplifies the development of bioinformatic pipelines. In the long run, it will certainly be possible to automate the steps implemented in our workflow so that the algorithms can be applied to a large set of secondary metabolites when fed with an appropriate training set.

Method

Construction of phylogeny tree

The phylogenetic tree depicted in Figure 1e was constructed utilizing the PhyloPhlAn3 pipeline⁶¹. PhyloPhlAn is a comprehensive pipeline that encompasses the identification of phylogenetic markers, multiple sequence alignment, and the inference of phylogenetic trees. In this analysis, we employed over 400 universal genes defined by PhyloPhlAn as our selected phylogeny markers. Subsequently, the taxonomic cladogram was generated using the iTOL web tool (http://huttenhower.sph.harvard.edu/galaxy/).

Apply our pipeline to identify errors in the substrate information associated with the A domain of Pseudomonas in the MIBIG database

We downloaded all A domain and substrate data from all NRPS secondary metabolism instances within Pseudonomas from MIBIG 3.0³²(Supplementary_table4). Following our pipeline, we initiate the process by utilizing NRPSMotifFinder to extract all Amotif4-5 sequences of the A domain. Subsequently, we calculate the Jukes-Cantor sequence distance between these Amotif4.5 pairs. Finally, employing Ward linkage clustering, we conduct a cluster analysis on all Amotif4.5 sequences, resulting in the generation of Figure S2. We observed that within the same cluster, substrates often exhibit some very subtle impurities. We speculate that such clusters might contain some errors in substrate information. Consequently, we focused on these potential errors, manually searched for their original literature, and verified that the A domain and its corresponding substrate information entered in MIBIG do indeed contain errors (Supplementary_table5).

Apply our pipeline to predict the structural composition of diverse NRPS metabolites across various genera of Burkholderiales

To evaluate the extendibility of our bioinformatic pipeline to various metabolites and microbial genera, we utilized a previously curated dataset of A domain sequences with experimentally confirmed substrates from Burkholderiales ⁶² (Supplementary_table7). This manually curated dataset, encompassing 7 genera, 34 secondary products, 203 A domains, and 25 substrates within Burkholderiales, serves as a valuable resource and mainly as a test set. Simultaneously, we compiled the training set by aggregating A domain and substrate data from all NRPS secondary metabolism instances within Burkholderiales, utilizing experimentally validated biosynthetic gene cluster data from MIBIG 3.0³² (Supplementary_table8). This training set covers 6 genera, 19 secondary products, 99 A domains, and 21 substrates. To ensure the inclusivity of all substrates from the test set, we randomly selected an A domain from each substrate in the manually curated dataset, incorporating it into our training set. Consequently, the finalized training set comprises 124 A domains (Supplementary_table12), while the test set encompasses 178 A domains (Supplementary_table9). Then, we concurrently employed our pipeline and antiSMASH¹⁶ to predict the corresponding substrate based on the sequence data of the A domain in the test set. Subsequently, we compared these predictions with the actual substrate information, calculating the accuracy rate for each algorithm.

Apply our pipeline to annotate the receptor FpvA across the bacterial domain

To explore the distribution of FpvA across the bacterial domain, we acquired all complete bacterial genomes (33,207) from NCBI (https://www.ncbi.nlm.nih.gov). Employing the pipeline depicted in Figure 5h to examine the CDS of these genomes, we identified a total of 357,790 TonB-dependent receptors, comprising 31,936 for FpvA and 623 for FpvB. Subsequently, leveraging the taxonomy information from NCBI, we conducted a comparative analysis of the distribution differences between FpvA and FpvB among various species. The findings revealed that FpvA exhibited a widespread distribution across different species, whereas FpvB was exclusive to the genus Pseudomonas, with a maximum occurrence of one in a given genome.

Data Availability

The source code and parameters used are available in the supplementary material.

Acknowledgements

We thank Richard Allen for genome sequencing of the 20 strains. We also appreciated Vera Vollenweider for sample preparation for the experiment of pyoverdine structure elucidation.

Funding

This work was supported by the National Key Research and Development Program of China (No. 2021YFF1200500, 2021YFA0910700), National Natural Science Foundation of China (No. 42107140, No. 32071255, No.41922053, No. T2321001), National Postdoctoral Program for Innovative Talents (No. BX2021012). R.K. is supported jointly by a grant from the Swiss National Science Foundation no. 310030_212266. V-P.F. is supported jointly by a grant from UKRI, Defra, and the Scottish Government, under the Strategic Priorities Fund Plant Bacterial Diseases program (BB/T010606/1), Research Council of Finland, and The Finnish Research Impact Foundation.

Author contributions

Shaohua Gu performed the majority of computational analysis in this research and drafted the manuscript. Yuanzhe Shao built the pipeline for retrieving synthetase sequence information and developed the phylogeny-centered method for predicting substrates. Karoline Rehm and Laurent Bigler performed the experiment of pyoverdine structure elucidation. Di Zhang depicted the circular plot. Ruolin He standardized of all NRPS structures into motif-intermotif format. Jiqi Shao and Ruichen Xu assisted in revising the manuscript. Xiaoying Bian provided curated dataset of A domain sequences with experimentally confirmed substrates from Burkholderiales. Xiaoying Bian, Alexandre Jousset and Ville-Petri Friman offered insightful comments and assisted in revising and writing of the manuscript. Rolf Kümmerli and Zhong Wei oversaw the project, designed experiments and revised the manuscript. Zhiyuan Li conceptualized and oversaw the project, conducted the analysis of the receptor annotated method, and revised the manuscript.

Competing interests

The authors declare no competing interests.