Low-quality sequences and adapter sequences were removed using Trimmomatic. Transcriptome sequencing reads were compared to the reference genome (NCTC10258: https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_900475405.1/) using Bowtie2 software. Gene expression levels were assessed by the number of sequencing reads aligned per kilobase base length per million sequencing reads (FPKM). RNA-seq counts were modeled under a negative binomial distribution, with significance determined by an exact test (p-values adjusted for false discovery rate [FDR] via Benjamini–Hochberg). Genes with |log₂(fold change)| > 2 and FDR-adjusted p-value (padj) < 0.05 were classified as significant.

Low-quality sequences and adapter sequences were removed using Trimmomatic. The assembled genome was subjected to gene prediction and exhaustive annotation work, including functional resolution and taxonomic calibration, with the help of the Prokka tool. A circos plot of the genome was created using the Proksee online tool to visualize the structural features and functionality distribution in the genome in a clear and intuitive manner. GATK software was used to compare the sequencing results with the reference genome, and single nucleotide polymorphisms were identified for subsequent analysis. SNPs between different samples were compared, SNPs that differed within samples were identified, and key SNPs were labeled. The copy numbers of the efflux pump genes were estimated using CCNE-acc(https://github.com/biojiang/ccne), a component of the Carbapenemase-Encoding Gene Copy Number Estimator (CCNE) tool. We added related efflux pump genes to ccne/config/CARD_AMR_clustered.csv to enable targeted analysis. CCNE-acc assumes that read coverage for each position of the whole genome and the target gene-containing fragment follow a Poisson distribution, and it estimates the target gene copy number by this equation: copy number(CCNE-acc)_{target gene} = (λ_{target gene})/(λ_Genome), where λ_{target gene} and λ_Genome are the parameters of two Poisson distributions and are approximated by the medians of read coverages of target gene-containing fragments and the whole genome.

The MIC values of eravacycline for the three strains WJ_4, WJ_14, and WJ_18 were 4 mg/L, 0.5 mg/L, and 1 mg/L, respectively, prior to the induction treatment and increased to ≥16 mg/L after in vitro induction with eravacycline and subsequent drug-free subculture for 10 generations ( Table 1 ). The time required to induce a satisfactory resistant strain was approximately 40 days. This result indicates that we successfully and stably induced high MICs for eravacycline in three strains that originally had low MICs. To identify these changes, we named them WJ_4a, WJ_14a, and WJ_18a, respectively.

To identify the genes that may be responsible for eravacycline resistance, we performed differential gene expression analysis of the induced eravacycline-resistant strains and original strains. In the WJ_4 group and its resistant strain WJ_4a, we identified 121 upregulated genes and 438 downregulated genes ( Figure 1A ); in the WJ_14 group and its resistant strain WJ_14a, we identified 1,910 upregulated genes and 755 downregulated genes ( Figure 1B ); in the WJ_18 group and its resistant strain WJ_18a, we identified 500 upregulated genes and 413 downregulated genes ( Figure 1C ). Despite the different origins of these three strains, all of them were induced to acquire eravacycline resistance. Therefore, the differentially expressed genes may be the factors that contributed to the development of this resistance; specifically, it is possible that genes showing the same change trend in all three groups are the main cause. Therefore, we identified the genes that were simultaneously upregulated or downregulated in all three groups. We detected a total of 31 genes whose expression was upregulated in all three groups ( Figure 1D , Table 2 ), of which four were drug efflux transporter-related genes (smrA, smeD, smeE and smeF). We also detected a total of 21 genes whose expression was downregulated in all three groups ( Figure 1E , Table 3 ); however, no resistance- or drug efflux transporter-related genes were found among these genes. Despite this, the global regulator cyclic AMP receptor protein (CRP) expression was co- downregulated in all three groups, suggesting that resistance may be synergistically enhanced through multidimensional mechanisms. For example, to release the inhibition of efflux pump, to promote the formation of biofilm, and to affect other drug-resistant pathways (Liu et al., 2025).

To investigate the gene mutations associated with eravacycline resistance, we sequenced the genomes of the three resistant strains and the corresponding sensitive strains and focused on genes significantly upregulated in the resistant strains at the transcriptome level, especially drug efflux transporter-related genes. Interestingly, SNPs at the genome level were found only in the three significantly upregulated drug efflux transporter-related genes; the remaining genes did not exhibit mutations. These three genes, smeD, smeE and smeF, were identified as having 403 single-nucleotide polymorphisms in WJ_14a compared with the sensitive control WJ_14, including 359 synonymous and 44 nonsynonymous mutations, the latter of which may have a greater impact on protein structure. Specifically, the smeD gene had 8 nonsynonymous mutations, smeE had 16, and smeF had 20. Twenty-five single-nucleotide variants were identified in WJ_18a, including 22 synonymous and 3 nonsynonymous mutations, with smeD, smeE, and smeF having one nonsynonymous mutation each. Notably, no SNPs were found in WJ_4a ( Supplementary Tables S1 , S2 ). Moreover, no screened SNPs were shared between WJ_14a and WJ_18a.

To elucidate whether genomic amplification contributes to eravacycline resistance, we quantified the copy number (CN) of the smeDEF efflux operon in induced-resistant strains versus their susceptible progenitors using the CCNE-acc algorithm. Significant CN amplification of smeDEF was observed in all resistant strains (P < 0.001, paired t-test), with distinct magnitudes across isolates ( Table 4 ): WJ_4a exhibited a 4.2-fold increase, WJ_14a showed a 2.4-fold increase, WJ_18a displayed a 7.2-fold increase. This amplification aligns with: (1) Transcriptional upregulation of smeDEF; (2) the absence of nonsynonymous SNPs in WJ_4a, suggesting CNV-driven overexpression as a primary resistance mechanism in this strain; (3) the heterogeneity of resistance evolution pathways, where WJ_14a and WJ_18a combined CNV with coding mutations.

In addition to resistance genes, we studied the following SNPs. (1) Regulatory regions (Promoters/UTRs): No mutations were found in conserved regulatory motifs (e.g., smeT-binding sites) within 1.5 kb upstream of smeDEF. (2) Accessory genes (Transcriptional Regulators): A synonymous mutation in smeT (WJ_14a: c.213G>A) showed no impact on smeDEF promoter activity.

As to the potential SNPs in CRP binding sites. Sm exhibits distinct regulatory mechanisms compared to model organisms like E. coli. No observed mutations: Our whole-genome SNP analysis (Snippy pipeline) of resistant vs. susceptible strains revealed no nonsynonymous SNPs or regulatory variants in predicted CRP binding motifs upstream of efflux pump genes (e.g., smeDEF) or other CRP-associated loci (cAMP metabolism genes).

WJ_97 is resistant to minocycline, levofloxacin, and TMP-SMZ, has intermediate resistance to ceftazidime, and has MIC₉₀ values of 8 mg/L for tigecycline, 64 mg/L for doxycycline, 16 mg/L for eravacycline, 4 mg/L for moxifloxacin, and 64/2 mg/L for ticarcillin/clavulanic acid.

We performed an exhaustive annotation of the WJ_97 genome after completing its assembly. This included gene annotation using Prokka, annotation of CRISPR sequences, and annotation of CARD resistance genes. These annotation efforts provided a comprehensive and in-depth understanding of the genes in the WJ_97 genome and their functions, confirmed the presence of CRISPR sequences, and reveal possible drug resistance genes. WJ_97 possesses eight genes that are associated with drug resistance. Specifically, the efflux pump -associated genes carried by WJ_97 include smeD, smeE, smeF, smeR, and adeF in the RND family; tet(A) in the major facilitator superfamily (MFS); qacJ in the SMR family; and vanY, which is associated with structural changes in antibiotic target sites. Figure 2 shows the ring structure of the WJ_97 genome. We labeled important features such as genes, GC content, CDS regions, and resistance genes in detail to better represent the genome architecture. The scale on the genome ring clearly indicates the size of the genome, and the rich labeling of genomic features in the figure further deepens our understanding of the WJ_97 genome. This figure not only displays the structural features of the WJ_97 genome but also provides a valuable reference for exploring its biological properties and potential drug resistance. Through this comprehensive analysis of the WJ_97 genome, we were able to gain a deeper understanding of its biological mechanisms and functions, which provided an important direction for further research.

We compared the predicted protein sequences of all the genes derived from the CARD database and analyzed the drug resistance genes in the genomes of the sequenced strains ( Figure 3 ). We categorized these genes according to the mechanism of resistance and found that the highest percentage of them were genes encoding antibiotic-inactivating enzymes (41.4%). Efflux pump accounted for 34.4%, target alteration accounted for 13.8% and other mechanisms accounted for 10.4%, such as permeability reduction, biofilm formation and undefined pathway.

The WJ_97 genome was compared with the reference genome NCTC10258. NCTC 10258, as an international standard strain, possessed genetic stability and publicly available antibiotic resistance phenotype data, providing a baseline reference for understanding the acquisition mechanisms of clinical antibiotic-resistant strains. 549 SNPs were found in the coexpressed up-regulated genes and 28 SNPs were found in the coexpressed down-regulated genes. Sixty-five nonsynonymous mutations were found in the SNPs of the upregulated genes. We focused on the SNPs in smeD, smeE, and smeF, which are related to the resistance of the mutant strain to eravacycline, among which 16 nonsynonymous mutations were found in smeD, 28 nonsynonymous mutations were found in smeE, and 21 nonsynonymous mutations were found in smeF. In addition, there were 4 single-base deletion mutations in smeF. As deletion mutations that are not multiples of 3 result in amino acid frameshifts, thus having a greater impact on protein function, we believe that these mutations in the smeF protein in WJ_97 may be responsible for the emergence of resistance, although this needs to be verified by further experiments. The remaining 480 SNPs were synonymous mutations. Among the 28 SNPs in the coexpressed down-regulated genes, only two nonsynonymous mutations were found in the CRP gene, while the rest were synonymous mutations ( Supplementary Table S3 ).

A comparison of all resistance genes in WJ_97 with those in the mutant strains WJ4a, WJ14a, and WJ18a revealed a total of 313 SNPs. These variants included 8 deletions, 48 nonsynonymous mutations (21 in smeR, 7 in smeD, 8 in smeE, 6 in smeF, and 6 in adeF), 256 synonymous mutations, and 1 termination mutation (in the smeE gene) ( Supplementary Table S4 ).

There were synergistic activations of RND family pump smeDEF in the induced strains: Overexpression of smeD (1.4-3.2-fold upregulation), smeE (2.1-3.5-fold), and smeF (2.0-3.0-fold) was positively correlated with elevated MIC values. This finding is consistent with previous models of minocycline and tigecycline resistance (Chang et al., 2015; Cornely OA et al., 2018), further supporting the centrality of the RND family of pumps in tetracycline efflux. Our study provides the first evidence that genomic amplification of the smeDEF operon is a core mechanism conferring eravacycline resistance in Sm. Key findings support this: The amplification of smeDEF directly correlated with eravacycline MIC elevation, indicating a gene dosage effect where increased pump density enhances drug efflux capacity (Yaghoubi et al., 2022; Nicoloff et al., 2024). WJ_4a achieved high-level resistance (MIC=32 mg/L) without coding mutations in smeDEF, demonstrating that CNV alone can drive clinically relevant resistance. This challenges mutation-centric models and aligns with preadaption of intrinsic resistance genes under antibiotic pressure (García-León et al., 2014). Heterogeneity in CNV magnitude (e.g., WJ_18a: 7.2-fold vs. WJ_14a: 2.4-fold) reflects genetic background-dependent adaptation. Strains with lower CNV gains compensated via coding mutations (e.g., 44 nonsynonymous SNPs in WJ_14a), while WJ_18a combined moderate CNV with minimal mutations, indicating functional redundancy in resistance optimization (Andersson and Hughes, 2010).

CNV-mediated resistance may be reversible upon antibiotic withdrawal due to potential gene copy number instability, informing “drug holiday” strategies to restore baseline susceptibility—unlike stable SNP-based resistance (Nicoloff et al., 2024).

Based on the results, we analyzed regulatory and auxiliary factors influencing eravacycline resistance. CRP downregulation was consistently observed across all induced-resistant strains. CRP is a pivotal transcriptional regulator in bacterial antimicrobial resistance (AMR), modulating resistance through diverse pathways. CRP-cAMP complexes bind DNA to influence gene expression, enhancing resistance by altering membrane permeability, regulating efflux systems, and reducing antibiotic uptake. Additionally, CRP represses SOS response pathways, increasing persistence to β-lactams. These mechanisms underscore CRP’s role in fine-tuning bacterial AMR profiles in response to environmental stresses (Liu et al., 2025). In addition to affecting the efflux pump, its effect on biofilm formation is also worth considering. We found that the expression of the pga operon ( Table 2 ), which affects biofilm formation, was elevated in the in vitro induced strains (Kishii et al., 2020; Lin et al., 2020). It is necessary to further clarify in subsequent studies whether the drug resistance phenotype of the in vitro induced strain is related to biofilm formation.

Regarding auxiliary factors, no mutations were detected in promoter regions (e.g., smeT binding sites) or CRP-binding motifs, suggesting that resistance arises primarily from transcriptional dysregulation, genomic amplifications (CNVs) and resistance gene SNPs rather than sequence-specific alterations in regulatory elements. This emphasizes the dominance of post-transcriptional and structural changes in efflux pump genes (e.g., smeDEF) over direct promoter-driven adaptations (Andersson and Hughes, 2010).

Based on the results, the integrative resistance architecture in clinical isolate WJ_97 involves synergistic mechanisms, dominated by core efflux dysregulation through mutations in the smeDEF RND efflux complex (including a critical frameshift mutation in smeF predicted to alter substrate efflux) and dysregulation of adeF, which aligns with findings in induced strains and confirms smeDEF’s centrality. This is complemented by accessory mechanisms such as target alteration via vanY and inactivating enzymes (accounting for 41.4% of resistance genes), mechanisms largely absent in induced strains, alongside evolutionary divergence characterized by smeR mutations (21 nonsynonymous SNPs) and additional resistance genes like tet(A) and qacJ, reflecting adaptation under complex clinical antibiotic pressure. Ultimately, this multi-layered architecture demonstrates greater complexity than in vitro-induced resistance.