Soil samples were collected in four ice-free sites on King George Island (Antarctic Peninsula) during the austral summer of 2007 (Figure 1A). The sampling sites and collected samples include the following:

soil under Sanionia uncinata from the North Peak (62°04′849”S, 58°24′024”W; Figure 1B),

The method described by da Silva et al. (2017) was used to isolate culturable bacterial fraction. A preliminary evaluation of the 16S gene rrs indicated that the isolates belong to the genus Pseudomonas. The isolates are part of the Antarctic culture collection at the Microbial Molecular Ecology Laboratory (Federal University of Rio de Janeiro, Brazil).

Twenty-five psychrotolerant isolates affiliated with the genus Pseudomonas were selected for this study. Of these, 13 (52%) were isolated from ornithogenic soil, followed by the rhizosphere of the native plants C. quitensis (n = 6; 24%) and D. antarctica (n = 4; 16%). Two isolates (8%) were isolated from soil covered by the moss S. uncinata. The highest temperature at which we observed growth was 28°C. For this reason, the optimal incubation was 28°C for 24 h (antimicrobial resistance phenotypic screening) to 48 h (DNA extraction).

The genus of each isolate was confirmed using matrix-assisted laser desorption/ionization – time of flight (MALDI-TOF; Microflex LT, Bruker GmbH, Berlin, Germany). This experiment represented the beginning of the trial for Pseudomonas isolates (Figure 2), which was performed using the algorithm provided by the manufacturer. The colonies were transferred in triplicate to the plate “MSP 96 Polished Steel BC,” provided by the manufacturer. The plate was cleaned in accordance with the manufacturer’s instructions, with 70% alcohol followed by 80% trifluoroacetic acid. We added 1 μL of 70% formic acid (Tedia, Fairfield, Ohio, United States) and allowed the plate to fully dry at room temperature. Afterwards, 1 μL of the matrix α-cyano-4-hydroxycinnamic acid (Bruker GmbH, Berlin, Germany) diluted to 10 mg/mL in organic solvent [50% acetonitrile and 2.5% trifluoroacetic acid (Tedia, Fairfield, Ohio, United States)] was applied and dried at room temperature. The calibration strain was Escherichia coli ATCC 25922, which was the reference strain for the peaks in a spectrum of proteins between 2 and 20 kDa, as provided by the manufacturer (software FlexControl v.3.4, Bruker GmbH, Berlin, Germany). The results were compared with the spectra in MALDI Biotyper v 3.1, using the MBT Compass software and the MALDI Biotyper^® CA library (Bruker GmbH, Berlin, Germany).

The isolates were incubated at 28°C with constant shaking at 150 rpm, until the OD₆₀₀ reached 1 (approximately 48 h). Bacterial genomic DNA was extracted using the Wizard Genomic DNA Purification Kit (Promega, Madison, Wisconsin, United States), following the manufacturer’s instructions. The DNA was quantified using a Qubit fluorometer (Invitrogen, Waltham, Massachusetts, United States) with the Qubit double-stranded DNA high sensitivity Assay Kit (Life Technologies, Carlsbad, California, United States).

First, the genetic diversity of the isolates was assessed via a repetitive extragenic palindromic polymerase chain reaction (BOX-PCR). The final concentration of each reagent was 1.0 μM of primer BOXA1-R CTACGGCAAGGCGACGCTGACG (Versalovic et al., 1994), 1.25 u of GoTaq^® G2 DNA polymerase, 1X Green GoTaq^® reaction buffer, 1.5 mM of MgCl₂, 0.2 mM of deoxynucleotide triphosphate (Promega, Madison, Wisconsin, United States), and 50 ng/μL of genomic DNA with a final volume of 25 μL. The complete amplification cycle was one cycle of 95°C for 7 min, 30 cycles of 94°C for 1 min, 53°C for 1 min, 65°C for 8 min, and a final extension of one cycle at 65°C for 16 min in a thermocycler (Eppendorf, Hamburg, Germany). The products were analyzed using electrophoresis 1.5% agarose using the 1 kb DNA ladder (Thermo Fisher Scientific, Waltham, Massachusetts, United States). The run took place at 100 V for 30 min. The dendrograms were constructed using the program BioNumerics v.7, with default parameters (Biomérieux, Marcy-l’Étoile, France).

The amplification of the gene rrs was performed with 5 pmol/μL of each primer (27f AGAGTTTGATCATGGCTCAG and 1492r GTTTACCTTGTTACGACT) and a final fragment size of 1,465 base pairs (bp) (Lane, 1991). The reaction had a final volume of 50 μL, using the same concentrations of Taq, reaction buffer, MgCl₂, deoxynucleotide triphosphate, and genomic DNA described above. The cycle was performed as follows: one cycle of 94°C for 3 min, 35 cycles of 94°C for 40 s, 55°C for 1 min, 72°C for 2 min, and a final extension of one cycle at 72°C for 10 min (Lane, 1991). The PCR products were analyzed as described above.

Additionally, each amplicon was purified with the enzyme ExoSAP (Exonuclease I, Shrimp Alkaline Phosphatase; Thermo Fisher Scientific, Waltham, Massachusetts, United States). The enzyme was diluted to a ratio of 1:9 in nuclease-free water (Qiagen, Hilden, Germany). Next, the following cycle was performed: one cycle of 37°C for 15 min (enzymatic optimal temperature) and one cycle of 80°C for 15 min (denaturation). In addition, 5 μL of pure amplicon and 5 μL of each primer (at a final concentration of 5 pmol/μL) were inoculated in a microplate and submitted for Sanger sequencing (Macrogen, Seoul, South Korea). The sequences were evaluated, trimmed and aligned using BioEdit v.7.2 (Hall, 1999). The species were defined using the tool Sequence Match, with the nonparametric k-nearest neighbors’ method, available at the Ribosomal Database Project (http://rdp.cme.msu.edu, accessed on March 15, 2019).

Antimicrobial susceptibility was assessed using the disk diffusion method, according to the protocol M02 established by the Clinical Laboratory Standards Institute (CLSI) (2012a). The tests were performed with piperacillin (PIP, 100 μg), aztreonam (AZM, 30 μg), piperacillin-tazobactam (TZP, 110 μg), ceftazidime (CAZ, 30 μg), cefepime (FEP, 30 μg), imipenem (IPM, 10 μg), gentamicin (GM, 10 μg), norfloxacin (10 μg) and ciprofloxacin (5 μg). To check the evolutionary aspects, we also tested antimicrobials known to be ineffective against P. aeruginosa. The list included tetracycline (TE, 30 μg), sulfamethoxazole-trimethoprim (25 μg), NEO (30 μg), chloramphenicol (C, 30 μg), ertapenem (ETP, 10 μg), ampicillin (AM, 10 μg), amoxicillin-clavulanate (AMC, 30 μg), cephalothin (CF, 30 μg), and cefotaxime (CTX, 30 μg). These antimicrobials were selected to evaluate whether the Antarctic P. fluorescens isolates exhibited similar resistance patterns to P. aeruginosa.

To evaluate extended-spectrum beta-lactamase (ESBL) phenotypes, PIP, AZM, CAZ, and FEP disks were positioned 2.5 cm from TZP, while AM, CF, and CTX disks were placed 2.0 cm from AMC. The antimicrobial disks represented the product “sensifar” and were commercially obtained from Cefar (São Paulo, Brazil), except PIP, which was prepared using the lyophilized drug from MilliporeSigma (Burlington, Massachusetts, United States). In addition, P. aeruginosa ATCC 27853 was the positive control strain, and the data were interpreted using CLSI M100 [Clinical Laboratory Standards Institute (CLSI), 2023].

The minimum inhibitory concentration (MIC) of colistin (CL; MilliporeSigma, Burlington, Massachusetts, United States) was accessed using broth microdilution. Each 0.5 McFarland suspension was diluted with a cation-adjusted medium (0.2 mL of Ca²⁺ and 0.1 mL of Mg²⁺ to each 100 mL of Mueller Hinton; Difco Laboratories Inc., Detroit, Michigan, United States) [Clinical Laboratory Standards Institute (CLSI), 2012b]. Serial dilutions of CL, with concentrations between 0.032 and 256 μg/mL, were evaluated. The strains P. aeruginosa ATCC 27853 (MIC 0.5–4 μg/mL) and E. coli ATCC 25922 (MIC 0.25–2 μg/mL) were used as susceptible controls, and E. coli C153 (carrying mcr-1, MIC 8 μg/mL) was employed as the CL-resistant control.

The MIC of the beta-lactams CAZ, FEP, AZM, and aminoglycoside NEO (MilliporeSigma, Brulington, Massachusetts, United States) were analyzed using agar dilution [Clinical Laboratory Standards Institute (CLSI), 2012b]. The beta-lactams were dissolved and diluted according to the method in previous work [Clinical Laboratory Standards Institute (CLSI), 2023]. Moreover, NEO was dissolved and diluted with sterile distilled water to a final concentration of 50 mg/mL, following the manufacturer’s instructions. Serial plates with concentrations varying between 1 and 256 μg/mL were analyzed, and P. aeruginosa ATCC 27853 was the control. The results were interpreted using the MIC breakpoints for other non-Enterobacterales [Clinical Laboratory Standards Institute (CLSI), 2023]. When the breakpoints for other non-Enterobacterales were unavailable (as in the case of CL), P. aeruginosa breakpoints were used.

Fourteen non-clonal isolates were selected for WGS. The selection criteria aimed to include isolates with high-level resistance phenotypes toward beta-lactams (specifically CAZ and FEP) and CL, and those with the most susceptible profiles.

Paired-end sequencing libraries (2×150 bp; 450 bp insert size) were constructed using 5 μg/μl of genomic DNA, following the NEBNext Fast DNA Fragmentation and Library Preparation Kit (New England Biolabs Inc., Ipswich, Massachusetts, United States). The quality control analysis of the final libraries was performed using the 2100 bioanalyzer (Agilent Technologies, Santa Clara, California, United States) and was visualized using electrophoresis 1.2% agarose. All samples were sequenced on the Illumina Hi-Seq 2500 platform (Illumina, San Diego, California, United States).

The quality of raw sequences was evaluated with FastQC v.0.11.5 (Andrews, 2010), and the reads and adaptors were trimmed using fastp v.0.23.4 (Chen et al., 2018) with default quality filter of >Q15. The genomes were assembled using Unicycler v.0.5.0 (Wick et al., 2017) with tested k-mer sizes of 27,53,71,87,99,111,119,127. The quality assessment of each assembly was checked with QUAST v.5.2.0 (Gurevich et al., 2013), CheckM2 v.1.0.2 (Chklovski et al., 2024) and GUNC v.1.0.2. As a quality filter, we considered N50 > 70 Kb (QUAST), completeness >90% (CheckM2), contamination <5% (CheckM2) and clade separation score > 0.45 (no chimeric contig) (GUNC). MOB-suite (Robertson and Nash, 2018) was used to identify plasmids in the draft genomes. The plasmid database used is available online at: https://zenodo.org/records/10304948/files/data.tar.gz, accessed on August 25, 2024. Default parameters were used for each software.

To determine the taxonomic classification of each Pseudomonas strain based on their genomes, we performed an analysis using the Genome Taxonomy Database Toolkit (GTDB-Tk) v.2.3.2 (Chaumeil et al., 2022) with the Classify workflow (“classify_wf”) and database r214. In the “ani_screen’ step it uses Mash v.2.3 (Ondov et al., 2016) to find the best hits among the representative genomes in the r214 database, then FastANI v.1.32 (Jain et al., 2018) to identify the species of the query genome using ≥95% as Average Nucleotide Identity (ANI) cutoff (Arahal, 2014). If the ANI analysis does not identify the query genome species the next steps are performed. In the “identify” step it employs Prodigal v.2.6.3 (Hyatt et al., 2010) and HMMER v.3.4 (Finn et al., 2011; Finn et al., 2015) for the identification of 120 bacterial phylogenetic marker genes and performs a multiple sequence alignment. The “align” step concatenates and filters the alignment. Finally, the “classify” step uses pplacer v1.1.alpha19-0-g807f6f3 (Matsen et al., 2010) to determine the place of the genome in the GTDB-Tk reference tree. The genomes were annotated using Prokka v.1.14.6 (Seemann, 2014).

The average nucleotide identity (ANIb) based on Basic Local Alignment Search Tool+ (BLAST+) and the correlation indices of tetra-nucleotide (TETRA) signatures of all analyzed genomes was run using default parameters on Pyani (Pritchard et al., 2016). The results were combined to assess the relationship between the genomes using the R package v.4.0.3, using the Euclidean distance and dist function from the statistics package for distance calculations. The hclust function from the statistics package was applied, using the average method for clustering calculations. Further details are provided in the online manuals (available at: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/dist and https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/hclust, accessed on March 10, 2024).

A scatterplot of the correlation of ANIb and TETRA values was generated using the ggplot2 package (Wickham, 2016), and the correlation was evaluated using Spearman’s correlation and the cor.test function from the statistics package. The cor.test calculates an exact p-value when using “cor.test (clusteredAni$height, clusteredTetra$height, method = “spearman”).” The shapiro.test function from the statistics package was applied to assess the normality of the data, revealing a Gaussian distribution. The shapiro.test calculates an approximate p-value when using “shapiro.test (clusteredAni$height) and shapiro.test (clusteredTetra$height).” Further details are provided in the online manuals (available at: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/shapiro.test and https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/cor.test, accessed on March 10, 2024). Additionally, the confidence of each clade was calculated with a bootstrap of 100 replicates using pvclust 2.2–0 (Suzuki and Shimodaira, 2006), and the factoextra package was employed to generate the dendrogram. For the visualization of the graphs,

We ran Benchmarking Universal Single-Copy Orthologs (BUSCO) (Simao et al., 2015) and BUSCO Phylogenomics (Waterhouse et al., 2018) to create the supermatrix. Besides, MAFFT (Kuraku et al., 2013) was employed to align the sequences of the supermatrix. We used the “--auto” option from MAFFT, which selects the appropriate alignment strategy amongst FFT-NS-2, FFT-NS-i and L-INS-I, according to the size of the input data. Gblocks (Talavera and Castresana, 2007) extracted the best-aligned blocks, using “sequence type equals protein (−t = p)” and “minimum length of a block equals 5 (−b4 = 5)” as the parameters. The output was converted to the Phylogeny Inference Package (PHYLIP) format using ClustalW2 (Larkin et al., 2007). The phylogenomic analysis was performed using RAxML (Stamatakis, 2014), using 100 bootstrap repetitions. The substitution model was defined by “-mPROTGAMMAWAG,” in which the model of heterogeneity is “GAMMA” and the substitution model is “LG.” Further details are provided in the online manual (available at: https://cme.h-its.org/exelixis/resource/download/NewManual.pdf, accessed on March 10, 2024). The tree was visualized and colored using iTOL v.7 (Letunic and Bork, 2024).

Genome Unclutterer (GUNC) v1.0.6 (Orakov et al., 2021) with clade separation score (CSS) > 0.45 was used to detect chimeras, contamination and the annotation of plasmids. For the taxonomy curation, type strain genome server (available at Type Strain Genome Server (dsmz.de); accessed on March 23, 2024) and GTDB-Tk v2 were employed (Chaumeil et al., 2022).

Antimicrobial resistance and further genes of interest were identified using the Genome Annotation tool available at the Pathosystems Resource Integration Center (PATRIC; available at Bacterial and Viral Bioinformatics Resource Center | BV-BRC, accessed on June 24, 2024) (Davis et al., 2020). For the annotation of ARGs, we applied the keywords “resistance,” “beta-lactam,” “penicillin,” “aminoglycoside,” “aph,” “aac,” “chloramphenicol,” “cat,” “polymyxin,” “arn,” “pmrK,” “tetracycline,” “tetR,” “efflux,” “ABC transporter,” and “ABC-type” using default parameters for Bacteria/Archaea. “Isopenicillin N epimerase” was detected for most of the isolates, but it was not included due to its importance on the antimicrobial biosynthesis pathway (not described in this work).

The keywords “heavy metal,” “arsenic” and “prophage” were also investigated. Relevant genes for the adaptation to extreme environments were searched and the function was manually assigned based on the annotation (Wattam et al., 2017).

The following databases were searched for ARGs using the ABRicate v.1.0.1 Pipeline: MEGARes v.3.0 (6635 sequences) (Lakin et al., 2017; Bonin et al., 2023), ResFinder v.4.1 (3077 sequences) (Bortolaia et al., 2020; Florensa et al., 2022), NCBI AMRFinderPlus v3.12.8 (5386 sequences) (Feldgarden et al., 2021; Feldgarden et al., 2022), ARG-ANNOT (2223 sequences) (Gupta et al., 2014), and CARD v.3.2.4 (2631 sequences) (Alcock et al., 2023; McArthur et al., 2013). We used ABRicate v.1.0.1 with a minimum identity threshold of 80% and a minimum coverage threshold of 80%. Databases were queried sequentially and overlapping predictions from multiple databases were considered a validation of results and retained in the final analysis. ABRicate was run on the Galaxy version 24.1.4.dev0 server, where its user-friendly interface allows for parameter configuration and database selection without requiring complex workflows. All databases (DbType nucI) were last updated November 4, 2023.

This study employed the Genomic Island Prediction Software (GIPSy) to check for genomic island availability (Soares et al., 2016). GIPSy uses an eight-step workflow for genomic island prediction, with each step incorporating specific default parameters to identify genomic features: Step 1 processes input files; Step 2 applies a G + C content deviation cutoff of 1.5 standard deviations; Step 3 uses a sensitivity setting of 0.95 for codon usage deviation (Colombo/SigiHMM); Step 4 predicts transposase genes with an HMMer e-value of 0.0001; Step 5 detects virulence or resistance factors using BLASTP with an e-value of 0.000001; Step 6 performs reciprocal BLAST with an e-value of 0.000001; Step 7 identifies tRNA flanking regions with an HMMer e-value of 0.0001; and Step 8 combines results from all previous steps to predict pathogenicity and resistance islands. All steps were executed using default settings provided by the software. In cases where regions were associated with more than one type of genomic island (pathogenic, metabolic, resistance, or symbiotic), PAIs and RIs were retained simultaneously, with overlapping regions plotted together at the same locus on the circular genomic graphs to reflect their dual classification. This research applied BLAST Ring Image Generation (Alikhan et al., 2011) to represent and evaluate the position of genomic islands in various strains of Pseudomonas, and the similarity between strains of the same species. The strain P. antarctica LMG 22709 (NZ_LT629704.1) was selected as the reference genome for GIPSy Island predictions. The strain P. carnis BML-PP010 (BQHE01000001.1) was applied as a reference for pathogenic strain. Genomic ring images were generated, and for groups with more than one sample, the largest genome was selected as the central ring for the circular plots.

The rrs sequences encoding the 16S gene are provided under the nucleotide accession numbers MK681799.1 to MK681824.1. The complete genome sequence data, including raw sequence reads, genome assemblies, and annotations of Pseudomonads applied in this study, were submitted to GenBank under the BioProject accession PRJNA1183857. Supplementary Table S2 displays the nucleotide accession numbers of the isolates evaluated using phylogenomic analysis.

A screening with six distinct classes of antimicrobials was performed using the disk diffusion method in accordance with CLSI standards [Clinical Laboratory Standards Institute (CLSI), 2012a; Clinical Laboratory Standards Institute (CLSI), 2023]. Similar to the control strain P. aeruginosa ATCC 27853, the 25 isolates presented resistance to representatives of four classes of antimicrobials: C, sulfamethoxazole-trimethoprim, TE, and beta-lactams (AM, AMC, CTX, CF, and ETP). The evaluated isolates in this work did not present a halo distortion, suggesting the absence of ESBL.

All isolates presented a putative phenotype susceptible to NEO. The diameter of the inhibition zones varied between 19 and 26 mm, in contrast to the control strain (no observed halo).

Regarding the beta-lactams with breakpoints available for the clinical treatment of other non-Enterobacterales, all isolates presented phenotypes intermediate or resistant to AZM. Moreover, 11 isolates (50%) were intermediate/resistant to CAZ, six (25%) to FEP and one (4%) to IPM (Supplementary Table S1).

The MIC by agar dilution confirmed the low susceptibility to AZM, CAZ, and FEP, with MICs varying from 32 to ≥256 μg/mL, whereas P. aeruginosa ATCC 27853 displayed an MIC of 4 μg/mL for CAZ and AZM and 8 μg/mL for FEP. This method also confirmed the susceptibility of the Antarctic isolates to NEO, with MICs ranging from 2 to 8 μg/mL, whereas the control strain P. aeruginosa ATCC 27853 presented an MIC ≥256 μg/mL.

As stated by Clinical Laboratory Standards Institute (CLSI) (2023), an MIC ≤2 μg/mL is considered intermediate to CL, and a value ≥4 μg/mL is considered resistant. Likewise, all 25 isolates described in this work were intermediate or resistant to this polymyxin. The MIC was 0.5 μg/mL for four isolates, 1 μg/mL for seven, 2 μg/mL for two, and 4 μg/mL for one. Eleven isolates (44%) displayed higher MICs than the positive control E. coli C153 (16 μg/mL): two were 128 μg/mL, and nine were ≥ 256 μg/mL. Table 1 presents the phenotypic data for 14 isolates selected for WGS. Supplementary Table S1 presents the complete data for the 25 phenotypically analyzed isolates.

The MALDI-TOF analysis (scores between 1,703 and 2,259) confirmed that the isolates belong to the genus Pseudomonas. According to BOX-PCR, seven clusters ranged from two to four isolates, represented by clones (Supplementary Figure S1). Isolates representative of such genetic diversity or displaying diverse antimicrobial susceptibility profiles were selected for WGS.

The average genome size of the Pseudomonas was 6.5 Mb (varying between 5.8 and 7.6 Mb) and presented an average GC ratio of 59.7% (between 58 and 60%). In addition, an average of 6,135,071 protein coding sequences were identified. Integrated prophages were found across each genome, in which the most common encoded protein was Gifsy-2, followed by an antirepressor, CP4-57 regulatory, and Lp2 protein 6. Table 2 summarizes the general features of the draft genomes of the Pseudomonas isolates.

According to the phylogenomic analyses based on the average nucleotide identity (ANI) derived from the complete genome distance matrix and the correlation of tetra-nucleotide frequency (TETRA), the P. fluorescens isolates from Antarctic samples exhibited high similarity to each other (ANIb >98%). The isolates belong to the P. fluorescens group, forming two distinct clusters (Supplementary Figure S2). The first and largest cluster comprised the isolates C291, O329, and O62 from a recent common ancestor, whereas C290, S191, S101, D118, and O160 were distant from each other. In contrast, the second cluster displayed less divergence between the isolates O230, O39, D277, O11, D47, and O64. Moreover, ANIb and TETRA varied between 0.0 and 0.3, with an outlier after 0.4.

The taxonomy was curated using GTDB-Tk (Chaumeil et al., 2022) and the Type Strain Genome Server, and all isolates belong to the P. fluorescens group. In general, three potential new species were identified (Figure 3). Apart from O62 and O64, presenting a completeness of 99.99%, all genomes presented a completeness of 100%. Furthermore, all presented a very low contamination ratio [between 0 (C290) and 1.49 (D118)], indicating that contaminant contigs from other genomes were not identified.

This study investigated the genes potentially involved in TE (tetR), aminoglycoside [aph(3′) and aac(6′)], polymyxin (arnC, D, E, F), fosfomycin (fosA) and, to a lesser extent, C (cat). Unknown quinolone and fosmidomycin resistance proteins were also detected. The quinolone resistance protein, which could possibly confer resistance to the first generation (i.e., nalidixic acid; not assayed in this work), was identified in nine isolates and presented 100% coverage and over 99% identity with other P. fluorescens isolates. The penicillinases penicillin-insensitive transglycosylase (EC 2.4.2.-), transpeptidase penicillin-binding protein (PBP-1C), other PBPs, class C beta-lactamases (EC 3.5.2.6) and metallo-beta-lactamases were frequently identified (Figure 3; Table 3).

Hundreds of copies of genes encoding efflux pumps were predominantly found in each sample. In order of frequency, all five well-known efflux pump families were detected: ABC transporters (including the gene encoding the macrolide-specific efflux protein MacA), resistance nodulation division, major facilitator superfamily, multidrug and toxic compound extrusion, and small multidrug resistance proteins (Table 3; Supplementary Table S3). The overexpression of efflux pumps, alongside the activity of beta-lactamases, could be responsible for the observed phenotype of resistance towards AZM, CAZ and FEP. Further experiments with efflux pump inhibitors could help in the understanding of the resistance mechanism.

The beta-lactamase bla_PFM-2 responsible for carbapenem resistance, which was identified in only one of the tested isolates (strain D47), was identified using NCBI AMRFinderPlus. In the databases CARD and MEGARes, only efflux pump genes from the resistance nodulation division family were detected (Supplementary Table S3). No ARGs were detected using the databases ResFinder and ARG-ANNOT. Moreover, only one plasmid was detected (isolate D118), which did not exhibit antimicrobial resistance determinants (Supplementary Figure S3). Overall, the data suggest that the variety of efflux pumps in the P. fluorescens genomes evaluated in this work plays a significant role in the observed antimicrobial resistance.

Apart from the isolate D118, PAIs and RIs were identified in all characterized genomes. The PAIs were more frequent, varying from six to 22, whereas up to nine RIs were identified. In general, the genomes described in this work presented higher numbers of PAIs and RIs than the pathogenic reference strain (P. carnis BML-PP010; BQHE01000001.1). The isolate D47 presented the highest total PAIs and RIs (n = 29), whereas S101 was even lower than the pathogenic reference strain (11 versus 13; Table 3). Supplementary Figure S4 presents the circular genome comparison plots displaying the islands identified against the genus Pseudomonas, reference strain P. antarctica LMG 22709 (NZ_LT629704.1) and pathogenic reference Pseudomonas carnis BML-PP010 (BQHE01000001.1).

Copies of genes encoding the heavy-metal response regulator, heavy-metal sensor histidine kinase and heavy-metal resistance transcriptional regulator (hmrR) were identified in the genomes. One copy of a membrane-bound cytochrome biogenesis cycz-like domain, annotated as a heavy-metal associated domain, was identified in each genome. In addition, genes encoding arsenic resistance proteins were also found, varying from zero to five copies.

The DNA repair system genes (recA, recX, mutL, and mutS) responsible for the resistance to ionizing radiation were identified in each genome. Likewise, lexA (signaling and regulation) and kdpD (osmotolerance) were also observed. Genes conferring resistance to high pressure or high temperature (dnaK, groEL, dnaJ, grpe, groES, htpX, and rpoH) were also detected, and two genes responsible for the resistance to low temperatures (cspD and cshA) were found. The genes pstA, pstB, and pstC involved in the optimization of phosphorus uptake were detected. The frequency of these bacterial adaptation genes was generally low (from one to three copies), except for dnaJ, where the genomes presented three to five copies each.