To analyze the drug resistance status of bacteria isolated from animals and the environment of the animal farms, 658 bacterial strains were obtained from the poultry and livestock anal fecal swabs, and sewage and soil of the animal farms in Wenzhou, China. An fosA-like gene, designated fosA13, was found in an isolate from poultry named DW0548. Bacterial species identification was performed via 16S rRNA gene homology and genome-wide average nucleotide identity (ANI) analyses. Table 1 lists the plasmids and strains used in this study.

The MICs of the antimicrobials were determined via the agar dilution method according to Clinical Laboratory Standards Institute (CLSI) guidelines (CLSI, 2024). Medium with glucose-6-phosphate (G6P) at a constant concentration of 25 μg/mL was used when the MIC of fosfomycin was tested. E. coli ATCC 25922 was used as a quality control. Table 2 shows the MIC data for the 25 antibiotics from the six antibacterial categories.

The bacteria were cultured in liquid LB medium and incubated at 37°C for 16 h. Total DNA was extracted from the bacteria with the Generay Genomic DNA Small Volume Preparation Kit (Shanghai Generay Biotech Co., Ltd., Shanghai, China). Whole-genome sequencing was completed on the Illumina NovaSeq and PacBio RS II platforms at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China). The Illumina short reads and the PacBio long reads were assembled using MEGAHIT v1.2.9 (Li et al., 2016) and Trycycler v0.5.1 (Wick et al., 2021), respectively. Using Pilon v1.24, the quality of the genome sequence from PacBio sequencing was corrected by mapping the Illumina short reads onto the PacBio read assembly (Walker et al., 2014). The open reading frames (ORFs) were predicted with Prokka v1.14.6 (Seemann, 2014). The functions of the predicted proteins were annotated by searching the ORFs against the NCBI nonredundant protein database with DIAMOND v2.0.11 (Buchfink et al., 2021). The promoter region of a gene was predicted using the BPROM tool (http://www.softberry.com/berry.phtml?topic=bprom&group=programs&subgroup=gfindb). Annotation of drug resistance genes was conducted by using Resistance Gene Identifier v5.2.0 (RGI) and the Comprehensive Antibiotic Resistance Database (CARD) database (https://github.com/arpcard/rgi) (McArthur et al., 2013). Homology analysis of the 16S rRNA gene from the target genome was performed by comparison with the 16S ribosomal RNA sequence database in NCBI (McArthur et al., 2013). The ANI was computed using FastANI v1.33 (Jain et al., 2018). Calculation of digital DNA−DNA hybridization (dDDH) values was performed on the basis of the Type strain Genome Server (TYGS) online database (Type Strain Genome Server) (Lian et al., 2021). Multiple sequence alignment of fosA13, fosA, fosA2, and other related genes was performed with MAFFT v7.487 (Katoh and Standley, 2013). A neighbor-joining (N-J) phylogenetic tree including FosA13 and other functionally characterized fosfomycin resistance enzymes was constructed with MEGA11 (Kumar et al., 2018). The phylogenetic tree of FosA13 with other Fos proteins was visualized with the online website iTol (iTOL: Interactive Tree Of Life) (Letunic and Bork, 2021).

Referring to a previous publication (Zhao et al., 2023), primers to clone the predicted resistance gene with its promoter region and the primers to clone the open reading frame (ORF) of the gene were designed. The PCR products were amplified and then inserted into the vectors pUCP19 and pCold I, respectively, via the DNA ligation kit Ver.2.1 (Takara Bio, Inc., Dalian, China). The recombinant plasmids pUCP19-fosA13 and pCold I-fosA13 were transformed into E. coli DH5α and E. coli BL21, respectively. The inserted sequences in the recombinants were verified by Sanger sequencing. Table 3 shows the primer sequences and details.

The methods used to express and purify the recombinant FosA13 protein were mainly based on a previous publication (Gao et al., 2022). In brief, the recombinant strain (pCold I-fosA13/BL21) was cultured in LB broth. When the OD value of the culture reached between 0.45 and 0.55, the expression of the protein was induced by 1 mM isopropyl-β-dithiogalactopyranoside (IPTG). Bacteria were harvested by centrifugation (4,000 × g, 10 min) at 4°C, resuspended in 5 mL of nondenaturing lysis buffer, and fragmented by ultrasonication for 10 min at 4°C. After centrifugation, the supernatant containing the recombinant protein was collected. The recombinant protein was purified with BeyoGold His-tag purification resin using a nondenaturing elution solution from a His-tag protein purification kit (Beyotime, Shanghai, China). The His-tag of the recombinant protein was removed by enterokinase (EK enzyme). The purity of the protein was determined via SDS−PAGE, and the protein concentration was subsequently determined using a BCA protein concentration assay kit (Beyotime, Shanghai, China).

The kinetic parameters of purified FosA13 with fosfomycin were analyzed via high-performance liquid chromatography (HPLC) using a Thermo Scientific AcceLA HPLC system (Thermo Fisher Scientific, Inc., China) with a 100 μL final reaction volume at 37°C. The 100 μL reaction system consisted of 10 mM GSH, 250 µm MnCl₂, 100 mM KCl, and 100 mM Tris-HCl, and gradient concentrations of fosfomycin (0, 25, 50, 100, 200, 400, and 800 µM) were added in a volume of 10 μL. After the reaction system was preheated for 20 minutes, the purified FosA13 protein was added. The reaction was conducted for 5 minutes, and the assay was carried out on the analysis system. The mobile phases were a mixture containing K₂HPO₄ (200 mM), KH₂PO₄ (200 mM), acetonitrile and methanol at a percentage volume ratio of 1.68: 78.32: 10: 10 (Arca et al., 1990; Bernat et al., 1997, Bernat et al., 1999). The steady-state kinetic parameters k _cat and K _m were nonlinearly regressed against the initial reaction rate via the Michaelis−Menten equation in Prism (version 8.0.2) software (GraphPad software, CA, United States). The value is the average of three independent measures.

The GenBank accession numbers for the novel fosfomycin fosA13 gene, the plasmid pMMDW0548 and the chromosome sequences of DW0548 are PQ600006, CP173708 and CP173707, respectively.

To clarify the mechanisms of antimicrobial resistance in bacteria isolated from animals and the environment, we sequenced 658 bacterial isolates isolated from poultry and livestock anal fecal swabs and environmental samples from animal farms in Wenzhou, China. Annotation of genomic data revealed resistance genes against different antibiotics. Notably, among the predicted genes were numerous putative fosfomycin resistance genes, including but not limited to homologs of fosA, fosB, fosC, and fosX. These genes shared less than 80.0% aa sequence identity with functionally characterized fosfomycin resistance genes. Some of these genes, including the fosA2, fosA5, fosB, fosLC2 and fosL1 homologs, were randomly selected and cloned, and their resistance functions were determined. Finally, a gene homologous to fosA2 (designated fosA13 in this work) that conferred resistance to fosfomycin was identified, and it was encoded in the chromosome of an isolate named DW0548.

To confirm the resistance function of this gene, the ORF with its promoter region was cloned into the vector pUCP19, and the recombinant plasmid containing fosA13 was transformed into E. coli DH5α cells. The transformant strain carrying fosA13 (pUCP19-fosA13/E. coli DH5α) showed a 4-fold increase in the MIC value to fosfomycin (2 μg/mL) compared with that of the control strain (pUCP19/DH5α, 0.5 μg/mL) ( Table 2 ). The MICs of the previously cloned fosfomycin resistance genes varied. Compared with the recipients, the cloned fosY (Chen et al., 2022) and fosA6 (Guo et al., 2016a) genes increased the MIC levels by 16- and 32-fold, respectively, whereas fosI (Pelegrino et al., 2016) and fosA7 (Rehman et al., 2017) increased the MIC levels by 128- and >256-fold, respectively ( Supplementary Table S1 ).

16S rRNA gene homology analysis revealed that the 16S rRNA gene of the isolate DW0548 shared the highest similarity (96.0% coverage and 99.0% identity) with that of M. morganii M11 (NR_028938.1). In addition, ANI analysis of all 759 Morganella genomes in the NCBI database revealed that 78 of these genomes had ≥ 95.0% ANI (the threshold value to define a bacterial species) with the isolate DW0548 genome. Seventy-three of these genomes were from the M. morganii genomes. The result of the digital DNA–DNA hybridization (dDDH) analysis revealed that the isolate DW0548 presented the highest dDDH value (75.7%) with M. morganii NBRC 3848, which was greater than the threshold (70.0%) for classifying a bacterial species. Therefore, on the basis of the results above, the isolate DW0548 was ultimately included in the species M. morganii and was thus designated M. morganii DW0548.

The whole genome of M. morganii DW0548 consists of a chromosome and a plasmid named pMM548. The chromosome size was approximately 4.31 Mb, and the average GC content was 50.35%, with 4,380 coding sequences (CDSs). The plasmid was 96,349 bp in length with the average GC content of 50.45%, and it encoded 109 CDSs ( Table 4 ).

The result of the susceptibility test demonstrated that M. morganii DW0548 had high MICs (≥ 8 μg/mL) for 11 of the 25 antimicrobials tested, which included members of 5 classes of antimicrobials, with 4 aminoglycosides (streptomycin, kanamycin, spectinomycin and paromomycin), 4 β-lactams (cefazolin, cefothiophene, cefoxitin and cefuroxime), 1 quinolone (naphridixic acid), tetracycline, and chloramphenicol ( Table 2 ).

In examining the relationship between the drug resistance phenotype and genotype, we found that even though the bacterium presented high MICs for antimicrobials from the 5 classes tested, only two genes (a β-lactam resistance gene, bla _DHA-16, and a phenicol resistance gene, catII), which shared ≥ 80% aa similarity with functionally characterized antimicrobial resistance genes (ARGs), were identified. No aminoglycoside, tetracycline or quinolone resistance genes were found in the whole-genome sequence ( Supplementary Table S2 ). Similar to the isolate DW0548, other M. morganii strains have been reported to be resistant to β-lactams, aminoglycosides, tetracyclines, fluoroquinolones, fosfomycin, and other types of antibiotics (Liu et al., 2016).

A comparison of FosA13 with these functionally characterized proteins in the CARD revealed that it had the greatest aa sequence similarities with FosA (55.6%), followed by FosA2 (55.2%), FosA3 (55.2%), FosA4 (55.2%), FosA5 (55.2%), FosA6 (55.2%), FosA7 (52.9%), FosA7.5 (52.9%), and FosA8 (52.9%). Sixteen function and structure essential residues of FosA (Beharry and Palzkill, 2005a) are conserved in FosA13. Half of them act as ligands for Mn²⁺ (His7, His64 and Glu110) and K⁺ (Ser94, Ser98 and Glu95) and ligands (Arg119 and Tyr100, within the hydrogen-bonding site of fosfomycin) involved in FosA binding to fosfomycin, and the other half of the residues (W34, Y39, W46, C48, Y65, D103, H107, Y128) were located in the putative fosfomycin/GSH binding channel. Although FosA13 shares only about 50–60% identity with these functionally characterized FosA proteins, the active site residues essential for FosA function remain unchanged ( Figure 1 ). Further analysis of the evolutionary relationship between FosA13 and different glutathione transferases revealed that FosA13 was most similar to FosA2 and formed a new branch ( Figure 2 ).

The fosA13 gene is 432 bp in length and encodes a 143 aa protein (FosA13). The predicted molecular weight of the mature glutathione-S-transferase is 16.48 kDa, with a pI of 6.08. The purified FosA13 has the ability to hydrolyze fosfomycin, with a K _m of 0.427 ± 0.007 µM, k _cat of 6.43 ± 0.04 s^-1 and k _cat /K _m of (1.50 ± 0.02) × 10⁴ M^-1. s^-1. In terms of enzyme kinetics, compared with the other two FosA proteins, FosA13 showed 6-fold lower hydrolytic activity against fosfomycin than that of FosA (Beharry and Palzkill, 2005b) (k _cat /K _m of 1.5 × 10⁴ M^-1). s^-1 vs. 9.0 × 10⁴ M^-1. s^-1) and approximately 5-fold lower than that of FosA6 (Guo et al., 2016) (k _cat /K _m of 1.5 × 10⁴ M^-1. s^-1 vs. 0.3 × 10⁴ M^-1. s^-1).

To investigate the distribution of fosA13-homologous genes, the nucleotide sequence of fosA13 was used as a query to search for similar genes in the NCBI nucleotide database. As a result, a total of 82 similar genes with similarities between 87.3% and 100.0% were obtained, all of which were derived from the M. morganii genomes ( Supplementary Table S3 ). Among these 82 genes, only one had a similarity of 87.3%, and all the others had similarities of ≥ 92.13%. The gene context of the fosA13-homologous genes was further analyzed. The 20-kb sequences with the fosA13-homologous genes at the center were intercepted and clustered. Finally, these 83 20-kb fragments (including the gene identified in this study) were grouped into 13 clusters with a similarity threshold of 90.0%. The gene context of 13 sequences consisting of one from each of the 13 clusters was compared ( Figure 3 , Supplementary Table S3 ).

As mentioned above and illustrated in Figure 3 , the fosA13 and fosA13-homologous genes and their surrounding sequences were relatively conserved in the M. morganii genomes. Further structural analysis revealed that no mobile genetic element (GME) was present in their flanking regions. Upstream of fosA13 are genes encoding proteins related to nifj [ferredoxin (flavodoxin) oxidoreductase], fldI [phenyllactate dehydratase activator] and fld [(R)-phenyllactyl-CoA dehydratase beta subunit], whereas downstream of fosA13 are genes related to pcpR (PCP degradation transcriptional activation protein), yijE_1 (putative cystine transporter YijE_1) and perR (HTH-type transcriptional regulator PerR). However, many other fos genes, such as fosC2 (Wachino et al., 2010), fosA3 (Wachino et al., 2010), fosA5 (Ma et al., 2015), fosA6 (Guo et al., 2016) and fosA8 (Poirel et al., 2019), are related to MGEs and are encoded on plasmids.