Integrated FTIR and Whole-Genome Sequencing Reveal Scale-Dependent Genotype–Phenotype Relationships in Multidrug-Resistant Pseudomonas aeruginosa

TL;DR

This study combines genome sequencing and infrared spectroscopy to explore how genetic differences in drug-resistant Pseudomonas aeruginosa relate to observable traits.

Contribution

The study introduces a novel integrated approach using WGS and FTIR to reveal scale-dependent genotype-phenotype relationships in multidrug-resistant P. aeruginosa.

Findings

Genomic analysis showed a conserved resistance backbone supplemented by acquired resistance genes and variable virulence factors.

FTIR spectroscopy captured phenotypic variation linked to regulatory and metabolic adaptations rather than deep phylogeny.

Scale-dependent agreement between genomic and phenotypic data supports FTIR as a rapid screening tool for outbreak surveillance.

Abstract

Multidrug-resistant Pseudomonas aeruginosa is a major cause of healthcare-associated infections, particularly in high-burden clinical settings where rapid tools to capture clinically relevant resistance and virulence phenotypes are needed. In this study, we applied an integrated whole-genome sequencing (WGS) and Fourier-transform infrared (FTIR) spectroscopy approach to evaluate genotype–phenotype relationships in multidrug-resistant P. aeruginosa isolates collected during the COVID-19 pandemic. High-quality WGS data were used to characterize antimicrobial resistance determinants, mobile genetic elements, and virulence gene repertoires, while FTIR spectroscopy provided culture-based phenotypic fingerprints reflecting cell envelope composition. Genomic analyses revealed a conserved efflux-centered intrinsic resistance backbone, variably supplemented by acquired β-lactamases and aminoglycoside-modifying enzymes, alongside a largely conserved core virulome with heterogeneity driven primarily by type III secretion system effector profiles. Comparison of FTIR- and WGS-derived distance matrices revealed a weak but statistically significant global association, indicating a non-linear relationship between genomic relatedness and phenotypic similarity. Cluster-level concordance was strongly scale-dependent, with high agreement emerging only at finer clustering resolutions, consistent with FTIR capturing phenotypic variation linked to regulatory, metabolic, and cell envelope adaptations rather than deep phylogenetic structure. Together, these findings show that multidrug resistance and virulence in P. aeruginosa are shaped by a modular genomic architecture that manifests as distinct, measurable phenotypic states. The observed scale-dependent concordance supports FTIR spectroscopy as a rapid, cost-effective phenotypic screening tool for outbreak-oriented surveillance, complementing WGS in integrated antimicrobial resistance monitoring workflows.

Figures

Full text

1. Introduction

Pseudomonas aeruginosa is a highly adaptable Gram-negative opportunistic pathogen and a prominent member of the ESKAPE group, representing one of the leading causes of healthcare-associated infections worldwide [1]. Its clinical relevance is largely driven by a combination of intrinsic multidrug resistance and an exceptional capacity for adaptive evolution, which together severely limit therapeutic options. During the COVID-19 pandemic, the incidence of multidrug-resistant P. aeruginosa (MDR P. aeruginosa) increased markedly, coinciding with prolonged hospitalization, extensive antibiotic exposure, and the widespread use of invasive medical devices that facilitated the selection and persistence of resistant lineages [2].

The success of P. aeruginosa as both a persistent colonizer and an invading pathogen is rooted in its genomic plasticity and complex resistance architecture [3]. Core intrinsic resistance mechanisms include multiple resistance–nodulation–division (RND) efflux systems (such as MexAB-OprM, MexXY-OprM, MexCD-OprJ, and MexEF-OprN), which together form a conserved efflux-centered backbone mediating baseline multidrug resistance [4,5]. These systems are tightly regulated by global and local transcriptional regulators, including ParRS, CpxR, and ArmR, allowing rapid phenotypic adaptation under antimicrobial pressure rather than simple binary resistance switching [6,7,8]. In parallel, P. aeruginosa deploys a broad array of virulence determinants, including the type III secretion system effectors exoU and exoS, the alginate biosynthesis gene algD implicated in biofilm formation, the flagellin gene fliC essential for motility, and iron-acquisition systems such as pvdA that support survival in iron-limited host environments [9,10,11].

Whole-genome sequencing (WGS) provides high-resolution insight into these determinants and is widely used for clonal tracking, outbreak investigation, and antimicrobial resistance surveillance [12]. However, genomic similarity does not necessarily translate into phenotypic similarity, as gene regulation, metabolic state, and environmental adaptation substantially shape bacterial phenotypes [13]. Fourier-transform infrared (FTIR) spectroscopy offers a rapid, phenotypic typing method that detects biochemical variation in cell envelope composition [14]. Previous studies have shown that while FTIR-based clustering may partially reflect genomic relatedness, discordance between FTIR and WGS often reveals phenotypic and regulatory diversity rather than analytical inconsistency [15,16]. Compared to whole-genome sequencing, FTIR spectroscopy offers a substantially lower-cost and faster alternative for bacterial typing, making it particularly attractive for routine diagnostic and surveillance settings where sequencing capacity or turnaround time may be limited.

Integrating WGS with FTIR-based phenotyping, therefore, provides a robust framework to link genetic determinants with observable cellular characteristics [14]. In this study, we applied an integrated WGS–FTIR approach to MDR P. aeruginosa isolates collected during the COVID-19 pandemic to (i) characterize the architecture of antimicrobial resistance and virulence and (ii) assess the scale-dependent concordance between genomic relatedness and FTIR-derived phenotypic similarity.

2. Materials and Methods

2.1. Bacterial Isolates and Study Setting

A total of 26 MDR P. aeruginosa clinical isolates were included in this study. All isolates were recovered from routine diagnostic specimens at the Department of Medical Microbiology, University of Szeged, Hungary, a tertiary-care academic diagnostic laboratory serving a large regional hospital network. The laboratory provides routine microbiological diagnostics for both inpatient and outpatient services, including intensive care units, surgical wards, hematology–oncology departments, and transplant units. Species identification was performed as part of routine clinical diagnostics using matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF MS). Isolate metadata are provided in Supplementary Table S1.

One Pseudomonas mosselii isolate was also included as a non-P. aeruginosa internal control to validate species-level discrimination in downstream genomic and phenotypic analyses.

2.2. Genomic DNA Extraction and WGS

WGS was performed following the workflow described earlier by Kaszab et al. [17]. Briefly, isolates were cultured aerobically on blood agar (Biomerieux, Lyon, France), and genomic DNA was extracted using the Quick-DNA fungi/bacterial kit (Zymo Research, Irvine, CA, USA), according to the manufacturer’s instructions. DNA quantity and quality were assessed spectrophotometrically.

Illumina-compatible sequencing libraries were prepared using the Nextera XT DNA library preparation kit, and paired-end sequencing was carried out on an Illumina platform [18].

2.3. Genome Assembly, Quality Assessment, and Annotation

Raw sequencing reads were subjected to quality control using FastQC [19], followed by adapter trimming and quality filtering [20,21,22]. Low-quality reads and adapter sequences were removed, and read correction was applied prior to assembly. Default parameters were used unless otherwise specified.

High-quality reads were assembled into contigs using MEGAHIT v1.2.9 (automatic k-mer size selection) and SPAdes v.3.15.3 (error correction turned off) [19]. The assemblies were merged using GAM-NGS [20]. Assemblies generated by the two algorithms were merged to improve contiguity and robustness. Assembly quality metrics, including total genome size, contig number, N50, and GC content, were calculated with QUAST v5.0.2 [21]. Genome completeness and contamination were assessed using BUSCO v.5.8.2 and CheckM v1.2.4, respectively [22,23]. Only high-quality draft assemblies were retained for downstream analyses. Genome annotation was performed using Prokka [24]. Antimicrobial resistance (AMR) genes and virulence factors were predicted using ABRicate (CARD and VFDB databases) [25]. Contigs harboring at least one ARG were screened for integrative mobile genetic elements (iMGE) using MobileElementFinder located on the same contig and defined them as linked if their distance was less than 1000, 10,000, and 50,000 base pairs (bps) [26]. Distance thresholds of 1 kb, 10 kb, and 50 kb were used as commonly applied windows to assess physical proximity between ARGs and mobile genetic elements. To evaluate the chance of spreading ARGs, the plasmid origin of contigs with ARGs linked to iMGEs was predicted using PlasFlow 1.1 [27].

All analyses described in this section focus exclusively on genome assembly, quality assessment, and annotation workflows; interpretation of resistome and mobilome composition is provided separately below.

2.4. Resistome and Mobilome Analyses

This section focuses on the composition and structural organization of the resistome and mobilome inferred from the annotated genome assemblies. Predicted antimicrobial resistance (AMR) genes were identified using CARD-based annotation and analyzed at both the antibiotic class level and individual gene occurrence level. Resistance determinants were grouped according to antibiotic class and resistance mechanism, and gene occurrence frequencies were normalized across isolates to enable comparative analyses [28,29].

To investigate the contribution of genetic mobility to resistance diversification, predicted AMR genes were assessed for their genomic localization and association with mobile genetic elements. Resistance genes were classified as chromosomally encoded or plasmid-associated based on their genomic context, and their linkage to integrative and other mobile genetic elements was evaluated to distinguish genome-embedded resistance from mobility-associated resistance components.

The chromosomal versus plasmid origin of resistance gene-containing contigs was inferred using PlasFlow based on sequence composition and genome signatures. To minimize misclassification arising from contig fragmentation or plasmid integration events, plasmid assignments were interpreted conservatively and limited to contigs with high-confidence plasmid probability scores. Ambiguous contigs were retained as chromosomal by default, and results were interpreted at the population level rather than at single-contig resolution.

Gene copy number was estimated from short-read whole-genome sequencing data based on normalized read depth across resistance gene loci. ‘Copy number’ refers to gene copy counts detected in draft genome assemblies and does not imply true genomic amplification. Read depth values were normalized to the median genomic coverage of each assembly to account for differences in sequencing depth between isolates. For visualization and comparative analyses, copy number was discretized into categorical values (0 = absence, 1 = single-copy presence, ≥2 = putative amplification), reflecting relative gene dosage rather than absolute copy number.

2.5. Virulence Factor Profiling

Virulence factor profiling was conducted using VFDB annotations [30]. Presence/absence matrices were generated for all the detected virulence-associated genes. To facilitate functional interpretation, virulence genes were aggregated into major functional categories (e.g., motility, secretion systems, iron acquisition, immune modulation), and category-level representation was normalized to the total number of genes assigned to each category.

2.6. FTIR Spectroscopy and Phenotypic Clustering

FTIR spectroscopy was performed following the protocol described earlier [31]. Strain typing was performed using the IR Biotyper (IRBT) system (Bruker Daltonics GmbH & Co. KG, Bremen, Germany) with spectral analyses performed in the mid-infrared fingerprint region (1800–900 cm^−1^), corresponding to the default spectral window for bacterial typing in the IR Biotyper platform. All isolates included in the FTIR analysis were cultured under standardized conditions to minimize phenotypic variability. For each isolate, FTIR analysis was performed using three independent biological replicates, each measured once, in accordance with the manufacturer’s recommended IR Biotyper workflow to ensure robust and reproducible spectral typing.

Briefly, strains were grown on Mueller–Hinton agar plates at 37 °C for 24 h. For spectral acquisition, approximately 1 µL of bacterial biomass was collected using a sterile inoculation loop and suspended in 50 µL of 70% ethanol in a 1.5 mL microcentrifuge tube using sterile metal rods supplied by the manufacturer. The suspension was homogenized by vortexing for 1 min, followed by the addition of 50 µL sterile distilled water and a further 1 min of vortexing to ensure uniform cell disruption and dispersion. Aliquots of the bacterial suspension (15 µL) were spotted onto an IRBT silicon target plate together with 12 µL of infrared test standards (IRTS1 and IRTS2; Bruker Daltonics GmbH & Co. KG, Bremen, Germany), which served as internal quality and calibration controls. The spotted samples were dried at 37 °C for 30 min until complete film formation was achieved. The prepared silicon plate was then inserted into the IRBT spectrometer, and spectra were acquired using default analytical parameters.

Spectral acquisition and preprocessing were performed using OPUS software version 7.5 (Bruker Daltonics GmbH & Co. KG, Bremen, Germany), while strain-level spectral clustering and analysis were conducted using the IR Biotyper software (client version 4.0.3; Bruker Daltonics). Only spectra fulfilling the manufacturer-defined quality criteria were included in downstream analyses. These criteria comprised acceptable absorbance intensity (0.4 < D value < 2.0 arbitrary units), signal-to-noise ratio (<150 × 10^−6^ AU), signal-to-water ratio (<300 × 10^−6^ AU), and fringe values (<100 × 10^−6^ AU). Spectra meeting all quality thresholds were classified as “quality pass” and retained for further evaluation.

Pairwise spectral distance matrices were computed from quality-pass FTIR spectra using the OPUS software, based on normalized infrared absorbance profiles within the manufacturer-recommended spectral window. FTIR-derived distance matrices were exported for downstream statistical analysis. WGS-based distance matrices were generated independently from genomic data as described above.

2.7. Distance- and Clustering-Based Concordance Analysis

Pairwise distance matrices were generated independently from WGS data and FTIR spectral profiles, resulting in two symmetric distance matrices with identical isolate ordering. The WGS-based distance matrix was calculated using wgMLST allelic distances, defined as pairwise allele differences across shared loci, based on allele profiles generated by the PGAdb-builder pipeline. To evaluate global concordance between genomic and phenotypic distance structures, a Mantel test was performed by calculating the correlation between corresponding elements of the two distance matrices. Distance-level correlation and cluster-level agreement calculations were implemented in Python (version 3.10; Python Software Foundation) using standard statistical definitions of the Mantel test and Adjusted Rand Index (ARI) [32,33]. Correlation was assessed using Spearman’s rank correlation coefficient to account for potential non-linear relationships between genomic and phenotypic distances. Statistical significance was evaluated using a permutation-based approach, in which isolate labels were randomly permuted to generate a null distribution of correlation values under the assumption of no association between the distance matrices [32].

To assess agreement at the level of discrete classifications, hierarchical clustering was applied independently to the WGS- and FTIR-derived distance matrices using identical linkage criteria to ensure comparability. Clustering was based directly on the precomputed pairwise distance matrices and performed using the average linkage (UPGMA) method. Cluster assignments were compared across a range of clustering resolutions (k values). Agreement between clustering solutions was quantified using the ARI, which measures the similarity between two partitions while correcting for agreement expected by chance. ARI values range from 0 (random agreement) to 1 (perfect concordance). No additional distance transformation or optimal leaf ordering was applied, as cluster membership rather than dendrogram visualization was the primary focus of the concordance analysis. WGS-based pairwise distances represent genome-wide allelic distances derived from wgMLST profiles, whereas FTIR-based distances correspond to normalized spectral dissimilarities computed from IR Biotyper fingerprint spectra.

Graphical visualization of distance matrices, ARI profiles, and Mantel scatter plots, shown in Figure 5, was performed in Python within a Google Colab environment [34]. All other figures and graphical elements in the manuscript were generated using GraphPad Prism (version 10.0; GraphPad Software, San Diego, CA, USA).

2.8. Ethical Considerations

All bacterial isolates were obtained as part of routine diagnostic procedures. No patient-identifiable data were used, and no additional sampling was performed specifically for this study.

3. Results

3.1. Genome Sequencing and Assembly Quality

All 26 MDR P. aeruginosa clinical isolates produced high-quality draft assemblies using the SPAdes-MEGAHIT hybrid workflow (Supplementary Table S2). Assembly sizes ranged from 6.32 to 7.16 Mb, corresponding to the expected genomic range of P. aeruginosa. The number of contigs ≥ 200 bp varied between 596 and 1031 (median: 821), while the largest contigs were between 45.7 kb and 99.4 kb. GC contents were tightly clustered around the species baseline (65.6–66.3%).

N50 values ranged from 13.5 kb to 21.0 kb, with most isolates falling between 15 and 19 kb, consistent with typical draft-quality assemblies derived from short-read Illumina sequencing. The total length of contigs ≥ 0 bp closely matched the contig ≥ 200 bp dataset, indicating a negligible presence of small contigs or assembly artifacts. No isolate contained ambiguous nucleotides (N’s per 100 kbp = 0 for all samples), further supporting assembly integrity.

BUSCO and CheckM analyses confirmed > 98% genome completeness and <1% contamination for all P. aeruginosa isolates (data are shown in Supplementary Table S3). Mean sequencing depth ranged from 68× to 121×, ensuring robust detection of resistance genes, virulence loci, and mobile genetic elements (Comprehensive Antibiotic Resistance Database [CARD], Virulence Factor Database [VFDB], Integrated Mobile Genetic Elements database [IMGE]). wgMLST analysis revealed that all 26 Pseudomonas aeruginosa isolates exhibited unique allele profiles, indicating the absence of clonal duplication and a high level of genomic diversity within the dataset.

The internal sequencing control, the Pseudomonas mosselii isolate, displayed a substantially smaller genome (~5.50–5.52 Mb) and lower GC content (63.98%), clearly separating it from P. aeruginosa and validating both the sequencing and quality-control workflow.

3.2. Antibiotic Class-Level Distribution of Antimicrobial Resistance Determinants

At the antibiotic class level, the resistome of the isolates was dominated by β-lactam-, polymyxin-associated determinants, and efflux-mediated multidrug resistance systems (Figure 1). Nearly all P. aeruginosa isolates displayed high proportional coverage in these three categories, reflecting the conserved intrinsic resistance backbone of the species.

Aminoglycoside resistance showed pronounced inter-isolate variability, with several strains harboring a substantial complement of aminoglycoside-modifying enzymes, while others exhibited minimal representation, consistent with heterogeneous acquisition of accessory resistance genes. In contrast, determinants linked to sulfonamides/trimethoprim resistance were restricted to a limited subset of isolates, indicating sporadic and non-uniform distribution across the cohort.

The P. mosselii control strain (isolate #1) demonstrated markedly reduced representation across all antibiotic classes, clearly separating it from the P. aeruginosa isolates and validating the specificity of class-level resistome profiling.

3.3. Global Architecture of the Antimicrobial Resistome Revealed by Copy Number Analysis

Copy number-based analysis of CARD-annotated resistance determinants revealed a highly structured antimicrobial resistome across the 26 MDR P. aeruginosa isolates (Figure 2). Resistance genes clustered into clearly defined functional modules corresponding to efflux systems, regulatory elements, β-lactamases, aminoglycoside-modifying enzymes, and additional resistance-associated pathways.

Efflux-related determinants constituted the dominant and most conserved component of the resistome. Core RND efflux pump systems—including mexAB-oprM, mexXY-oprM, mexEF-oprN, and associated membrane fusion proteins (oprM, oprJ, opmB, opmE, and opmD)—were detected at high copy numbers in nearly all isolates. Regulatory genes controlling efflux activity and stress adaptation (armR, parR, parS, soxR, rsmA, and cprRS) were likewise broadly conserved, underscoring the central role of efflux-mediated intrinsic resistance in this cohort.

In contrast, β-lactam resistance determinants exhibited marked heterogeneity. Chromosomal ampC variants (PDC alleles) were present in all isolates, often with multiple copies, consistent with their role in baseline β-lactam resistance. Class D oxacillinases showed isolate-specific copy number patterns, with several strains harboring expanded OXA-family repertoires. Metallo-β-lactamase genes (including VIM-2, VIM-4, VIM-43, and IMP-79) were restricted to a subset of isolates but, where present, appeared at elevated copy numbers, highlighting their potential clinical relevance for carbapenem resistance. Elevated copy numbers of metallo-β-lactamase genes therefore indicate an increased genetic potential for carbapenem resistance, but do not imply a linear relationship with phenotypic resistance levels.

Aminoglycoside-modifying enzymes displayed substantial variability across isolates. Genes encoding acetyltransferases, nucleotidyltransferases, and phosphotransferases (e.g., aac, aad, ant, and aph variants) were unevenly distributed, with some isolates exhibiting multiple copies while others lacked acquired aminoglycoside resistance genes altogether. This pattern suggests that aminoglycoside resistance in the cohort arises from a combination of acquired enzymatic mechanisms and intrinsic efflux activity.

Additional resistance-associated pathways formed smaller but distinct modules. Polymyxin-related determinants (arnA, basS) were detected in most isolates, consistent with conserved lipid A modification capacity, while sulfonamide and tetracycline resistance genes (sul1, tet(A)) were confined to a limited number of strains. The disinfectant resistance gene qacEΔ1 occurred sporadically and at low copy number, indicating accessory acquisition rather than a core feature of the resistome.

Notably, the Pseudomonas mosselii control isolate (isolate #1) showed a markedly reduced resistance gene copy number across all functional categories, clearly separating it from the P. aeruginosa isolates and validating both species assignment and analytical specificity.

Overall, the copy number-resolved resistome analysis demonstrates that multidrug resistance in the P. aeruginosa cohort is driven by a conserved efflux-centered core architecture, upon which isolate-specific expansions of β-lactamase and aminoglycoside-modifying enzyme genes are superimposed. This modular organization highlights the balance between intrinsic resistance mechanisms and variably acquired determinants in shaping clinically relevant MDR phenotypes in P. aeruginosa.

3.4. Mobile Genetic Element-Associated Antimicrobial Resistance Determinants

Analysis of mobile genetic element-associated antimicrobial resistance determinants revealed a pronounced asymmetry between chromosomal and plasmid localization (Table 1). IMGE-linked resistance genes were predominantly associated with chromosomal regions, spanning a broad range of antibiotic classes, including β-lactams (including carbapenems), aminoglycosides, fluoroquinolones, sulfonamides, tetracyclines, and multidrug-associated resistance categories. This widespread chromosomal distribution suggests that insertion sequences, transposons, and integrative mobile elements play a major role in restructuring intrinsic resistance loci and facilitating local genomic diversification.

In contrast, plasmid-associated IMGE signals were markedly reduced in both abundance and spectrum. Plasmid-linked mobile resistance determinants were largely confined to a limited subset of antibiotic classes, indicating selective acquisition rather than broad horizontal dissemination. Notably, resistance categories associated with multidrug efflux regulation and disinfectant tolerance were almost exclusively chromosomal, underscoring the importance of genome-embedded mobile elements in modulating regulatory and permeability-related resistance mechanisms.

Overall, the IMGE landscape demonstrates that while plasmids contribute discrete, epidemiologically relevant resistance traits, chromosomal mobile elements constitute the dominant driver of resistance gene mobility and diversification in this isolate collection. This pattern supports a model in which genomic plasticity in Pseudomonas aeruginosa is primarily fueled by chromosomal IMGE activity rather than plasmid-mediated gene flux.

Given the use of short-read assemblies, these localization patterns should be interpreted as dominant trends at the cohort level rather than as absolute assignments for individual contigs.

3.5. Virulence Gene Profiles Revealed by VFDB Analysis

VFDB-based screening of the 26 MDR P. aeruginosa isolates revealed a dense and largely conserved virulence gene repertoire characteristic of the cohort. Analysis of the full VFDB presence/absence heatmap comprising 317 virulence-associated loci demonstrated that most classical virulence systems were uniformly represented across the isolate set. These included genes involved in motility (e.g., fliC, motAB, pilA, and pilT), adherence via pili and type IV fimbrial clusters, and alginate biosynthesis (algD, algU, and algR), as well as iron acquisition systems (pvd, fpv, and pch siderophore operons), quorum-sensing regulators (las, rhl, and pqs systems), and multiple secretion pathways (T2SS, T3SS, and T4SS-associated structural components). This functional virulence profiling revealed a clear gradient from highly conserved to variably represented virulence systems across the isolate set. Core functional categories—including effector delivery systems, exoenzymes, exotoxins, motility, and biofilm-associated functions—showed near-complete representation in most P. aeruginosa isolates, reinforcing their role as foundational components of the virulence architecture (Figure 3A). In contrast, virulence functions related to immune modulation exhibited marked inter-isolate heterogeneity, indicating substantial variability in host–pathogen interaction potential within the cohort.

Despite this overall conservation, discrete sources of variability were evident. The most prominent heterogeneity was associated with the type III secretion system effector repertoire, where exoU and exoS displayed a mutually exclusive distribution pattern (Figure 3B). A subset of isolates carried exoU, a cytotoxic genotype commonly linked to acute and severe infections, whereas the remaining isolates encoded exoS, which is more frequently associated with invasive behavior and chronic colonization. Additionally, though less pronounced, variability was observed in selected secretion- and motility-related loci, consistent with strain-level phenotypic diversity.

The Pseudomonas mosselii control strain (isolate #1) displayed a sharply reduced and distinct virulence profile in both gene-level and functional analyses, lacking the majority of P. aeruginosa-specific virulence determinants. This clear separation validates species assignment and confirms the specificity of the VFDB-based virulome profiling approach.

Taken together, these results demonstrate that the MDR P. aeruginosa isolates harbor a highly conserved core virulome complemented by a variably represented functional component, with heterogeneity driven primarily by T3SS effector composition and immune modulation-associated functions. This architecture reflects the genomic plasticity and adaptive potential characteristic of MDR P. aeruginosa in clinical settings.

3.6. Comparison of FTIR- and WGS-Derived Distance Matrices

Pairwise distance matrices derived from WGS and FTIR spectroscopy were visualized as heatmaps to compare isolate relatedness across the two approaches (Figure 4). The WGS-based genomic distance matrix exhibited a broad dynamic range of pairwise distances, with marked contrasts between closely related and more divergent isolate pairs. This wide distance spectrum reflects the high discriminatory power of genome-wide comparisons, capturing both fine-scale relatedness and deeper genomic divergence within the isolate collection.

In contrast, the FTIR-based spectral distance matrix displayed a narrower overall distance range, with values distributed within a more constrained interval. Despite this compressed scale, the FTIR matrix revealed pronounced local variability, indicating subtle but structured differences in phenotypic similarity among isolates. These patterns suggest that FTIR spectroscopy is particularly sensitive to fine-grained phenotypic variation, even when large-scale genomic divergence is limited.

Importantly, neither matrix exhibited a random organization. Instead, both showed structured distance patterns, supporting the presence of non-random relationships among isolates at both the genomic and phenotypic levels. To quantitatively evaluate the correspondence between these distance structures, distance-level and cluster-level concordance metrics were applied, as described below.

3.7. Comparison of FTIR- and WGS-Based Distance Structures and Clustering

Pairwise distance matrices derived from WGS and FTIR spectroscopy were compared at both the global distance level and the cluster level to assess concordance between genomic and phenotypic relatedness. A weak but statistically significant correlation was observed between FTIR- and WGS-based distance matrices, as assessed by the Mantel test using Spearman’s rank correlation (r = 0.217, p = 0.0001).

To further explore agreement at the level of discrete classifications, hierarchical clustering was performed independently on the WGS- and FTIR-derived distance matrices using identical linkage criteria. Cluster-level concordance was found to be strongly dependent on clustering resolution (Figure 5A). At low cluster numbers (k = 2–5), agreement between the two methods was minimal, with ARI values close to zero, indicating near-random overlap of cluster assignments. In contrast, a marked increase in ARI was observed at higher resolutions, with substantial agreement emerging at k ≥ 6 and reaching high ARI values at k = 8–10.

The global relationship between genomic and spectral distances underlying this scale-dependent clustering behavior is illustrated by the Mantel scatter plot (Figure 5B). Pairwise comparisons revealed a weak but consistent positive monotonic trend accompanied by substantial dispersion across the full range of distances, demonstrating that genomic proximity does not translate linearly into spectral similarity. Together, these analyses show that FTIR and WGS share a weak but significant global association while exhibiting pronounced concordance only at finer clustering resolutions. Inspection of pairwise distance relationships indicated that the observed dispersion was not driven by a small number of extreme outlier isolate pairs but rather reflected broadly distributed deviations across the dataset, consistent with heterogeneous phenotypic adaptation superimposed on shared genomic backgrounds.

Although the correlation coefficient was low, the statistically significant Mantel result indicates a non-random global association between genomic and phenotypic distance structures, consistent with a non-linear genotype–phenotype relationship rather than direct proportionality.

Consequently, FTIR-based clustering could not be directly benchmarked against repeated sequence types, and the observed spectral similarities reflect phenotypic convergence rather than a shared clonal background.

4. Discussion

This study provides an integrated genomic, resistome, mobilome, virulome, and phenotypic comparison of multidrug-resistant Pseudomonas aeruginosa clinical isolates within the multidrug-resistant cohort. By combining high-quality WGS with copy number-resolved resistance profiling, mobile genetic element mapping, virulence factor analysis, and FTIR-based phenotypic clustering, we aimed to dissect the structural organization of antimicrobial resistance and pathogenic potential while also assessing the extent to which genomic relatedness is reflected at the phenotypic level. From a practical perspective, FTIR may serve as a rapid and cost-effective first-line phenotypic screening tool, while WGS remains indispensable for high-resolution genomic characterization, resistance mechanism elucidation, and epidemiological reconstruction. This multi-layered approach enables a nuanced interpretation of P. aeruginosa adaptation in clinical settings, beyond single-mechanism or single-technology perspectives.

The consistently high genome completeness (>98%) and low contamination (<1%) observed across all Pseudomonas isolates provide a solid foundation for comparative resistome, mobilome, and virulome analyses. Assembly metrics, including genome size, GC content, and N50 values, fall well within the ranges reported for short-read Illumina-based P. aeruginosa draft genomes in large surveillance and population genomic studies [35]. The absence of ambiguous bases and the high sequencing depth ensured reliable detection of resistance determinants, regulatory elements, and mobile genetic elements, minimizing technical noise in copy number-based analyses.

The clear genomic separation of the Pseudomonas mosselii control strain further validates the analytical pipeline and confirms the discriminatory power of both WGS-based and functional profiling approaches. Similar control-based validation strategies have been emphasized as critical in comparative resistome studies to avoid species-level confounding [36].

At both the antibiotic class level and gene copy number resolution, the resistome of the MDR P. aeruginosa cohort was dominated by efflux-associated determinants, particularly RND-type systems (Figure 1). The near-universal presence and high copy number of mexAB-oprM, mexXY-oprM, mexEF-oprN, and associated membrane fusion proteins are consistent with the well-established paradigm that multidrug efflux constitutes the core intrinsic resistance mechanism of P. aeruginosa [5].

The restriction of metallo-β-lactamases to a subset of isolates mirrors global epidemiological patterns, where carbapenemase-positive P. aeruginosa lineages remain relatively uncommon but clinically critical [37,38]. Similarly, the uneven distribution of aminoglycoside-modifying enzymes reflects sporadic horizontal acquisition, consistent with previous population-level analyses demonstrating that aminoglycoside resistance in P. aeruginosa often results from a combination of intrinsic efflux activity and accessory enzymatic mechanisms [39,40]. Together, these findings support a modular resistome architecture, in which a conserved intrinsic core is variably augmented by high-impact acquired determinants—a model increasingly recognized in MDR P. aeruginosa epidemiology [41].

A key observation in this study is the predominance of chromosomal over plasmid-associated mobile genetic elements in shaping the resistance landscape (Table 1). While plasmids are often emphasized in discussions of antimicrobial resistance dissemination, P. aeruginosa is increasingly recognized as a species in which integrons, transposons, and insertion sequences embedded within the chromosome play a dominant role [42,43]. The broad chromosomal distribution of iMGE-linked resistance determinants across multiple antibiotic classes in the present cohort is consistent with this emerging view. Previous studies have shown that class 1 integrons and transposon-associated resistance cassettes are frequently chromosomally integrated into P. aeruginosa, enabling stable inheritance while retaining adaptive flexibility [42,43]. The limited plasmid-associated resistance signal observed here further supports the notion that, unlike Enterobacterales, plasmid-mediated gene flux plays a secondary role in P. aeruginosa resistance evolution. This has important epidemiological implications, suggesting that resistance persistence may be driven more by clonal expansion and chromosomal remodeling than by rapid plasmid exchange.

VFDB-based analysis revealed a dense and highly conserved core virulome across the MDR P. aeruginosa isolates, encompassing motility, quorum sensing, iron acquisition, secretion systems, and biofilm-associated functions (Figure 3A). This observation aligns with previous genomic surveys demonstrating that virulence determinants in P. aeruginosa are largely species-defining and deeply embedded in the core genome [44,45]. Notably, the mutually exclusive distribution of exoU and exoS mirrors well-characterized virulence stratification patterns. exoU-positive strains are frequently associated with acute, cytotoxic infections and poorer clinical outcomes, whereas exoS-positive strains are more commonly linked to invasive and chronic infection phenotypes [46,47]. The presence of both genotypes within the cohort highlights clinically relevant heterogeneity that may not be captured by resistance profiling alone (Figure 3B). The functional virulence category analysis further suggests that variability is concentrated in immune modulation-related functions rather than in foundational virulence systems. This supports a model in which P. aeruginosa maintains a stable pathogenic core while tuning host interaction strategies in a strain-specific manner.

From a biochemical perspective, the phenotypic variability captured by FTIR spectroscopy most likely reflects composite alterations in cell envelope composition, including lipopolysaccharide remodeling, membrane lipid composition, and extracellular polysaccharide production. In P. aeruginosa, adaptive resistance mechanisms such as lipid A modification associated with polymyxin resistance, increased alginate production linked to biofilm formation, and altered outer membrane protein expression driven by efflux pump regulation are all expected to influence infrared absorbance profiles. Consequently, FTIR-based clustering is not expected to mirror phylogenetic relatedness linearly, but instead integrates regulatory, metabolic, and structural adaptations shaping the expressed phenotype under clinical selective pressure (Figure 4). Although colony morphology can influence FTIR-based similarity assessments, colony phenotype information was not systematically recorded for the present isolate collection and could therefore not be evaluated as a contributing factor.

The weak yet statistically significant Mantel correlation between FTIR- and WGS-derived distance matrices observed in this study is consistent with prior work evaluating phenotype–genotype concordance in bacterial typing (Figure 5). Several studies have reported similarly modest correlations between whole-genome relatedness and spectroscopic or proteomic fingerprints, reflecting the complex, non-linear relationship between genotype and phenotype [48,49]. Importantly, this correlation should be interpreted as weak but non-random rather than indicative of a strong genotype–phenotype correspondence, in line with the known limitations of Mantel tests in high-dimensional biological datasets. The strong scale dependence of cluster-level agreement provides critical methodological insight. The emergence of high ARI values only at finer clustering resolutions suggests that FTIR spectroscopy captures subtle phenotypic variation potentially associated with recent genomic divergence or regulatory differences, rather than broad phylogenetic structure. Clustering resolutions were explored in an agnostic, non-phylogenetic manner to evaluate scale-dependent concordance rather than to assign explicit biological or epidemiological meaning to individual k values. Such distributed dispersion is consistent with regulatory and metabolic plasticity shaping FTIR phenotypes beyond what is captured by core-genome distance alone. This observation aligns with previous IR Biotyper and FTIR-based studies, which have emphasized the utility of these methods for outbreak-level discrimination rather than deep evolutionary inference [50,51,52].

The functional impact of resistance gene copy number variation is context-dependent and influenced by regulatory control, promoter strength, and cellular physiology. In the absence of quantitative expression or MIC correlation data, copy number expansion should be interpreted as an indicator of resistance potential rather than a direct proxy for phenotypic resistance.

From a clinical microbiology perspective, the scale-dependent concordance observed between FTIR- and WGS-based analyses has important practical implications. The emergence of high cluster-level agreement only at finer resolutions suggests that FTIR spectroscopy is particularly well suited for rapid, outbreak-oriented phenotypic stratification rather than for broad phylogenetic inference. In high-burden clinical settings, such as intensive care units during periods of increased antimicrobial use, FTIR-based screening could enable timely identification of phenotypically related multidrug-resistant P. aeruginosa isolates, guiding early infection control interventions while awaiting definitive WGS results. Within a tiered diagnostic framework, FTIR may therefore function as a rapid frontline tool for surveillance and preliminary clustering, with whole-genome sequencing reserved for high-resolution epidemiological reconstruction, resistance mechanism elucidation, and long-term surveillance. Such an integrated approach balances speed, resolution, and resource allocation, aligning well with the operational constraints of routine clinical microbiology laboratories.

The observed scale-dependent concordance between FTIR- and WGS-based distance structures suggests that FTIR spectra capture biologically meaningful phenotypic information, which may provide a conceptual basis for the future development of machine learning models aimed at predicting resistance or virulence phenotypes directly from spectral data. Importantly, this observation is intended as a conceptual future perspective and does not imply that predictive machine learning models were developed or validated within the scope of the present study. However, such applications would require dedicated model training and validation beyond the scope of the present study.

This study has several limitations that should be acknowledged. First, the number of analyzed isolates was limited, and all samples originated from a single tertiary-care center, which may restrict the generalizability of the observed resistance and phenotypic patterns. Second, genome assemblies were generated exclusively from short-read Illumina sequencing data, limiting the resolution of plasmid architecture and the precise structural context of mobile genetic elements. The use of short-read sequencing limits precise reconstruction of plasmid structures and may lead to ambiguous localization of mobile resistance determinants; however, the conservative classification strategy applied here minimizes systematic bias in chromosomal versus plasmid-level trends. Third, while FTIR measurements were performed under strictly standardized conditions, spectroscopic phenotypes remain inherently sensitive to growth conditions and metabolic state, which may influence inter-laboratory reproducibility. FTIR spectroscopy reflects potential phenotypic states of bacterial cells rather than directly measured functional phenotypes, and targeted experimental validation will be required in future studies. The inclusion of Pseudomonas mosselii primarily serves as a species-level control and does not allow a systematic evaluation of FTIR robustness beyond interspecies discrimination. Finally, the cross-sectional design of the study precludes direct inference on temporal evolutionary trajectories. Despite these limitations, the integrated multi-layered approach employed here provides robust insight into genotype–phenotype relationships in multidrug-resistant P. aeruginosa.

Taken together, these results support a complementary rather than a substitutive role for FTIR and WGS. While WGS provides a stable framework for lineage assignment and resistance mechanism identification, FTIR appears particularly sensitive to phenotypic modulation arising from regulatory changes, metabolic state, or cell envelope composition.

5. Conclusions

This study demonstrates that multidrug resistance in clinical P. aeruginosa isolates is organized around a conserved, efflux-centered intrinsic backbone that is variably supplemented by chromosomally integrated, high-impact resistance determinants. Resistance diversification is driven predominantly by genome-embedded mobile genetic elements rather than plasmid-mediated gene flux, highlighting a species-specific evolutionary strategy that promotes both stability and adaptability under sustained antimicrobial pressure. Despite a largely conserved core virulome, functionally relevant heterogeneity—particularly in type III secretion system effector profiles—indicates distinct pathogenic strategies within the cohort. The weak but scale-dependent concordance between FTIR- and WGS-based relatedness further underscores that phenotypic and genomic similarity capture complementary aspects of bacterial diversity.

These findings support the integration of rapid phenotypic tools with genome-based analyses to improve surveillance and characterization of multidrug-resistant P. aeruginosa in clinical settings. Future studies incorporating longitudinal sampling and multi-center isolate collections will be essential to further refine the role of FTIR spectroscopy within integrated genomic surveillance frameworks for multidrug-resistant P. aeruginosa.