The Effect of Iodine-Containing Nano-Micelles, FS-1, on Antibiotic Resistance, Gene Expression and Epigenetic Modifications in the Genome of Multidrug Resistant MRSA Strain Staphylococcus aureus ATCC BAA-39

Application of supplementary drugs which increase susceptibility of pathogenic bacteria to antibiotics is a promising yet unexplored approach to overcome the global problem of multidrug-resistant infections. The discovery of a new drug, an iodine-containing nano-molecular complex FS-1, which has proven to improve susceptibility to antibiotics in various pathogens, including MRSA strain Staphylococcus aureus ATCC BAA-39TM, allowed studying this phenomenon. Chromosomal DNA and total RNA samples extracted from the FS-1 treated strain (FS) and from the negative control (NC) cultures were sequenced by PacBio SMRT and Ion Torrent technologies, respectively. PacBio DNA reads were used to assemble chromosomal DNA of the NC and FS variants of S. aureus BAA-39 and to perform profiling of epigenetically modified nucleotides. Results of transcriptional profiling, variant calling and detection of epigenetic modifications in the FS variant were compared to the NC variant. Additionally, the genetic alterations caused by the treatment of S. aureus BAA-39 with FS-1 were compared to the results of a similar experiment conducted with another model organism, E. coli ATCC BAA-196. Several commonalities in responses of these phylogenetically distant microorganisms to the treatment with FS-1 were discovered, which included metabolic transition toward anaerobiosis and oxidative/osmotic stress response. S. aureus culture appeared to be more sensitive to FS-1 due to a higher penetrability of cells by iodine bound compounds, which caused carbonyl stress associated with nucleotide damaging by FS-1, abnormal epigenetic modifications and an increased rate of mutations. It was hypothesized that the disrupted pattern of adenine methylated loci within methicillin-resistance chromosome cassettes (SCCmec) may promote excision of this antibiotic resistance determinant from chromosomes while the altered pattern of cytosine methylation was behind the adaptive gene regulation in the culture FS. The selection against the antibiotic resistance in bacterial populations caused by abnormal epigenetic modifications exemplifies possible mechanisms of antibiotic resistance reversion induced by iodine-containing compounds. These finding will facilitate development of therapeutic agents against multidrug-resistant infections.


INTRODUCTION
The misuse of antibiotics, as well as inappropriate prescription and overuse, has led to strong selective pressure, resulting in the survival and wide distribution of drug-resistant pathogens that threaten the public health system globally (Barbosa and Levy, 2000;Ventola, 2015). Acquired drug resistance to antimicrobial agents has been a widely recognized problem (Tenover, 2001;Spellberg et al., 2004;Wall, 2019;Rocha-Granados et al., 2020). Bacteria can respond rapidly to environmental changes owing to their short growth cycles, which enable them to evolve and adapt rapidly and survive under unfavorable conditions (Clatworthy et al., 2007). Currently, the rate of resistant bacteria is increasing, whereas the development of new antibiotics has dramatically declined over the past two decades (Spellberg et al., 2004).
Various factors are involved in the emergence and spread of drug resistant bacteria. These factors include mutations, which can modify target proteins; the transfer of genetic material, known as horizontal gene transfer; selective pressure in healthcare and community settings, which facilitates the development and distribution of multiple resistant bacteria; and in some cases, inability to detect emerging resistance phenotypes (Wang et al., 2017;Kumar et al., 2019). Resistance is generally associated with a decrease in bacterial fitness. It was expected that the physiological cost for bacteria to maintain resistance genes in the absence of antibiotics would be large enough to select substantially against the drug resistance. However, this has not been generally observed for various reasons. Firstly, the fitness cost for bacteria to maintain resistance is not always large enough to select for loss of the resistance alleles; therefore, even after removal of the drug, resistance may remain in the population for an extended period (Sjolund et al., 2003;Andersson and Hughes, 2010). Compensatory mutations and/or genetic regulatory mechanisms can also compensate for the large fitness cost of resistance by only activating resistance in the presence of the drug (Durão et al., 2018;Lin et al., 2018). Lastly, resistance mutations may provide the resistant strains with a fitness advantage by conferring increased virulence (Roux et al., 2015).
With the increasing occurrence of multidrug-resistant bacteria, monotherapy treatment is gradually becoming less adequate, necessitating the use of drug combination therapies Abbreviations: BM, base modification; BMF, bipartite methylated fragment; m4C, N 4 -methylcytosine; m6A, N 6 -methyladenine; MBC, minimal bactericidal concentration; modA, epigenetically modified adenine; modC, epigenetically modified cytosine; modG, epigenetically modified guanine; modT, epigenetically modified thymine; MRSA, methicillin-resistant Staphylococcus aureus; PacBio, Pacific Bioscience technology; SCCmec, staphylococcal methicillin-resistance chromosome cassette; SMRT, single molecule real-time sequencing. (Nigam et al., 2014;Zhu et al., 2014). Therefore, it is of great importance to devise strategies that focus on the application of supplementary drugs to increase susceptibility to regular antibiotics by inhibiting bacterial growth while reversing the selection for resistance. Strategies to combat antibiotic resistance with combination drug therapy have been proven promising since the late 1940s. Co-administration of streptomycin and para-aminosalicylic acid showed reduced evolution of resistant Mycobacterium tuberculosis strains (Dunner et al., 1949). Drug combinations are also currently being used in most treatments of infectious diseases (Poulikakos et al., 2014). Drug-induced reversion of antibiotic-resistant pathogens into sensitive phenotypes is a prospective approach to target the mechanisms and evolution of bacterial resistance (Baym et al., 2016;Ilin et al., 2017;van Niekerk et al., 2018).
Various mechanisms that reduce or invert the selective advantage of antibiotic resistance have previously been studied. Administering antibiotics with chemicals that inhibit specific resistance mechanisms is generally accepted as the bestestablished approach to neutralize the evolutionary advantages of resistant strains. However, this does not necessarily create a competitive disadvantage and the relative prevalence of resistance within a bacterial population is not reduced. Negative selective pressure is required to reduce resistant strains in the population, even when the antibiotic is present (Sjolund et al., 2003(Sjolund et al., , 2005. More recent studies focused on using the evolutionary and physiological interactions between drugs not only to neutralize the selective advantage of resistant strains, but also to select against them actively by imposing a direct cost on resistance. According to Baym et al. (2016), combinations of drugs can interact to impose a direct cost on resistance in three ways. Firstly, co-administrated drugs can cause suppressive interactions, where when drug A suppresses drug B, bacteria that became resistant to drug A can be inhibited more strongly by drug B owing to loss of their protective effect. This causes the concentration of bacterial growth to be non-monotonic, and sensitive bacteria grow better than resistant bacteria at high concentrations of drug B when drug A is present. Secondly, co-administered drugs can interact synergistically, creating a concentration regime where sensitive bacteria can continue to grow, while the growth of resistant strains is inhibited. Thirdly, resistance to drug A can increase the efficacy of drug B owing to evolutionary trade-offs, allowing selection against resistant strains even when drug A is present (Baym et al., 2016).
Key drivers of evolution are mutations of an organism's DNA, which generate genetic variation (Hershberg, 2015). Bacteria are constantly faced with the challenge to maintain fitness in changing environments, subjecting them to various stresses, which include oxidative stress and DNA damage. In response to these stresses, bacteria alter their phenotypes by modulating gene expression levels. For a long time, it was generally believed that adaptable variations in bacterial populations are possible only through mutational processes. The role of mutations in heritable phase variations of bacterial populations is a certainty. Bacteria have phase-variable genes, referred to as contingency genes, that are highly mutable compared to housekeeping genes (Phillips et al., 2019a,b). Slipped-strand mispairing (SSM) is described as a mechanism that promotes phase variation by reversible mispairing of short repeat sequences between mother and daughter strands during DNA replication (Henderson et al., 1999;Low and Casadesús, 2008;Harhay et al., 2019). SSM is not the only mechanism of phase variation. There are other mechanisms involving homologous recombination, insertion/excision of transposons and site-specific recombination. Recently, the role of differential DNA methylation and other epigenetic mechanisms of phase variations enabling the adaptation of bacterial populations to harsh or changing environments including antibiotic treatment regiments were pointed out (van der Woude and Bäumler, 2004;Sánchez-Romero and Casadesús, 2020). DNA methylation enables bacteria to control the reversible ON/OFF switching of important genes epigenetically. An example of this includes the pyelonephritis-associated pili (pap) operon in uropathogenic E. coli controlled by DNA adenine methylase (Dam). Switching between the ON/OFF states of the papBA genes determines the binding of two proteins at two GATC sites, before and after the promoter. The operon is turned to the ON state when methylation occurs proximal to the promoter, and vice versa (Hernday et al., 2002;van der Woude and Bäumler, 2004). Another example includes the gene flu, encoding for the outer membrane protein Ag43 in E. coli. This gene is regulated at three GATC sites and its expression is repressed by the oxidative stress response protein, OxyR. Binding of OxyR to a GATC site masks this site, therefore blocking methylation by Dam, and turning expression to the OFF state (Haagmans and van der Woude, 2000). Other study showed that Dam methylases can modulate pathogenicity of Klebsiella pneumoniae (Fang et al., 2017) and regulate the expression of drug resistance related genes including multidrug efflux pumps that leads to a fast acquisition by E. coli of a resistance to various antibiotics through epigenetic DNA modifications (Adam et al., 2008). Application of PacBio SMRT sequencing demonstrated a crucial role of directed DNA methylation in micro-evolution of the outbreak clones of methicillin-resistant Staphylococcus aureus (MRSA) leading to attenuation of agr transcriptional regulator expression and upregulation of genes involved in stress response and biofilm formation (Sullivan et al., 2019).
As epigenetic modifications are so important for the virulence and antibiotic resistance of pathogenic bacteria, the possibility of using the phase variation processes and epigenetic modifications to induce and direct antibiotic resistance reversion in populations of multidrug resistant pathogens may be assumed; however, it has never been explored. Recently, a new medicine, FS-1, was introduced into clinical practice as a supplement to antibiotic therapy for drug-resistant tuberculosis (Ilin and Kulmanov, 2014;Ilin et al., 2017;Korotetskiy et al., 2017;Islamov et al., 2018;van Niekerk et al., 2018). FS-1 is an iodine-containing nanomicelle, which reverses the susceptibility of resistant pathogens to conventional antibiotics. Sequencing of extensively drugresistant strains of M. tuberculosis (XDR TB), which were reverted to the sensitive phenotype by treatment with FS-1, showed no mutations to explain this phenomenon (Ilin et al., 2015(Ilin et al., , 2017. Epigenetic mechanisms of antibiotic resistance reversion were hypothesized. Later, the antibiotic resistance reversion was reproduced on beta-lactam antibiotic resistant strain Escherichia coli ATCC BAA-196 TM (Korotetskiy et al., 2019a) and MRSA Staphylococcus aureus ATCC BAA-39 TM (Joubert et al., 2019), which are more suitable for laboratory experiments than XDR TB.
Staphylococcus aureus is a representative of Gram-positive bacterial human pathogens, responsible for a wide variety of clinical manifestations. The occurrence of antibiotic resistance in S. aureus has been studied widely in view of their ability to develop resistance to nearly any antibiotic (Chambers and Deleo, 2009), exemplifying the adaptive evolution of bacteria during the antibiotic era. S. aureus is a fast-growing bacterium and proper growth can easily be maintained in a laboratory, making it simple and inexpensive to work with. This makes S. aureus a desirable model for studying the evolution of drug resistance and establishing new approaches to combat the resistance by induced antibiotic resistance reversion. The aim of this research was to study antibiotic resistance phenotype, transcriptome and epigenetic modifications in the genome of multidrug-resistant Staphylococcus aureus ATCC BAA-39 cultivated with a sublethal concentration of the drug FS-1. Additionally, the differential gene expression and epigenetic modifications in FS-1 treated S. aureus ATCC BAA-39 were compared to the results obtained on the phylogenetically distant model organism, E. coli ATCC BAA-196 (Korotetskiy et al., 2020b).

Bacterial Cultures
The model multidrug-resistant microorganism Staphylococcus aureus ATCC BAA-39 TM was obtained from the ATCC collection and kept in a freezer at −80 • C. Bacteria were cultivated on Mueller-Hinton (MH) liquid or solid media (Himedia, India) (Vitko and Richardson, 2013).

Antibiotic Resistance Reversion
Bacteria were inoculated into test-tubes with 10 ml of liquid MH medium supplemented with FS-1 (450 µg/ml) that corresponds to 1/2 minimal bactericidal concentration (MBC) of the drug estimated for the strain S. aureus BAA-39. As the negative control, the culture was cultivated in the same medium without the drug. Test-tubes were incubated at 37 • C for 24 h and then 0.1 ml aliquots of the cultures (in average 1.2 × 10 9 CFU/ml, no significant difference in the optical density of the overnight NC and FS cultures was observed) were transferred to fresh tubes with the corresponding media. After 10 passages, the experimental and control bacteria were cross-inoculated in three repetitions into tubes with drug-containing and drugfree media for further overnight incubation (24 h), followed by DNA and RNA extraction. The scheme of the experiment is shown in Figure 1. The susceptibility of bacterial cultures was evaluated by measuring the growth inhibition zones around disks (Himedia, India) with the following antibiotics: amoxicillin, clindamycin, gentamicin, methicillin, oxacillin, imipenem, tobramycin, cefamandole and erythromycin. The bacterial cultures were denoted as susceptible (S), intermediate (I) or resistant (R) by comparing the experimental results with the threshold inhibition zones recommended for the given antibiotics in the CLSI (2010) and by Sarker et al. (2014).

Radiological Study of FS-1 Complexing With Chromosomal DNA
To carry out the radiological study, isotope 131 I was used for the synthesis of FS-1. The documented radioactivity of the isotope 20 MBq/ml was controlled at the beginning of the FIGURE 1 | The scheme of the antibiotic resistance reversion experiment. Different experiments were denoted as (FS − ) -FS-1 treated culture grown for RNA extraction on MH medium without the drug FS-1; (FS + ) -FS-1 treated culture grown on MH supplemented with 450 µg/ml; (NC + ) -negative control culture grown on MH supplemented with 450 µg/ml; (NC − ) -negative control culture grown on MH without the drug FS-1.
Bacteria were incubated in a thermo-shaker overnight (24 h) on MH medium at 37 • C. The culture growth was centrifuged at 5,000 g and the cells were resuspended in saline to achieve 1 × 10 9 CFU/ml. Aliquots of 0.4 ml of the marked FS-1 solution (total radioactivity of an aliquot was 270 Kbq) were added to 7.6 ml of the bacterial suspension and rigorously shaken. Samples were incubated for 1 h at 37 • C in a thermoshaker, Thermomixer comfort (Germany). After incubation, bacterial cells were washed from FS-1 by three rounds of centrifugation at 5,000 g and re-suspended in sterile saline. DNA samples were extracted from bacterial cells using PureLink Genomic DNA Kits (Publication Number: MAN0000601, Revision 2.0), following the manufacturer's recommendations. DNA purity and concentration were controlled by NanoDrop 2000C (United States) at wavelengths 260 and 280 nm. Then the DNA samples were resolved in scintillation solute (ULTIMA GOLD LLT, PerkinElmer) to a final volume of 5 ml. The residual radiation was measured by β-spectrometer Hidex 300 SL (Finland).

DNA Extraction and PacBio Sequencing
DNA samples were extracted from bacterial cells using PureLink Genomic DNA Kits (Thermo Fisher, MAN0000601, Revision 2.0) following the manufacturer's recommendations.
Samples were prepared in three repetitions according to a guide of preparation of SMRT bell templates for sequencing on the PacBio RS II System. The samples were sequenced in Macrogene (South Korea) on SMRT Cell 8Pac V3 cells using the DNA Polymerase Binding Kit P6 and DNA Sequencing Reagent 4.0 v2 following the SMR Tbell 20-kb library preparation protocol. The raw PacBio reads are available at NCBI SRA database under accession numbers SRX5852332 (run SRR9077059) and SRX5843551 (run SRR9067444) for NC and FS genomes, respectively.

Total RNA Extraction and Ion Torrent Sequencing
Isolation of total RNA was performed with the RiboPure Bacteria Kit (Ambion, Lithuania) according to the developer's guidelines. The quality and quantity of the resulting RNA were determined using the Qubit 2.0 Fluorometer (Thermo Scientific, United States) and Qubit RNA HS Assay Kit (Invitrogen, United States). Purification of total RNA from ribosomal RNA was carried out using the MICROBExpress Bacterial mRNA Purification Kit (Ambion, Lithuania), as recommended by the developer. The efficiency of the template-RNA purification was determined on the Bioanalyzer 2100 (Agilent, Germany) with the RNA 6000 Nano LabChip Kit (Agilent Technologies, Lithuania).
The RNA fragment library was prepared by enzymatic fragmentation with the Ion Total RNA Seq Kit V2 (Life Technologies, United States). Barcoding of the formed library was carried out with the Ion Xpress RNA-Seq Barcode 01-16 Kit (Life Technologies, United States), according to the manufacturer's instructions. RNA sequencing was performed in three repetitions for each condition using the Ion Torrent PGM sequencer (Life Technologies, United States) with the Ion 318 Chip Kit V2. Resulted 12 sets of the raw RNA reads (NC and FS variants cultivated on MH and MH with 450 µg/ml FS-1 in three repetitions, see Figure 1) were deposited at NCBI SRA database under accession numbers SRR10728551-SRR10728554.
Reads obtained were demultiplexed with FASTQCMCF prior to being trimmed on quality (Phred score = 21) using the Raw RNA-Seq Data Processing pipeline implemented in Unipro UGENE v35.1 (Golosova et al., 2014). Reads shorter than 30 bp were filtered out.

Genome Assembly and Annotation
The raw PacBio reads were converted into XML subreads by the function HdfSubreadSets of the SMRT Link v5.0.1 software package using the parameters set by default. The complete genome assembly was performed using polished_falcon_fat pipeline, available from the SMRT Link v5.0.1 software (Chin et al., 2013;Roberts et al., 2013). Two genome variants, NC and FS, were obtained as single contig sequences without gaps and ambiguities with the respective coverage of 966 and 906 overlapping reads per nucleotide. The completeness of the final assemblies was evaluated using the benchmarking universal single-copy orthologous (BUSCO) software (Simão et al., 2015). Genome annotation was performed using the RAST Server 1 (Aziz et al., 2008) and then manually corrected. The locations of horizontally transferred genomic islands were identified by the program SeqWord Genome Island Sniffer (Bezuidt et al., 2009). The same program was used to identify the replication origin and terminus on bacterial chromosomes by GC-skew between the leading and lagging strands (Freeman et al., 1998). Antibiotic resistance genes were predicted using CARD-RGI Web-based tool 2 (Alcock et al., 2020). This research project was registered in the BioProject database at NCBI under the accession number PRJNA480363. PacBio reads generated for this study are available from the BioProject Web-site. The assembled and annotated genomes of the variants NC and FS of Staphylococcus aureus ATCC BAA-39 TM were published in NCBI GenBank under the accession numbers CP033505 and CP033506, respectively.

Gene Expression Profiling
The differential expression was done using the R-3.4.4 software. Firstly, a reference index was built for each reference genome using the buildindex function available in the Rsubreads package (Bioconducter). The obtained RNA fragments were aligned to the relevant reference genome with the use of the function align. The aligned BAM files and relevant GFF annotation files were then used as input for the featureCounts function to obtain gene counts. The R packages DESeq2 (Bioconducter) and GenomicFeatures was then used for the differential expression analyses (Love et al., 2014). The DESeq2 algorithm normalizes read counts by gene lengths and by total number of reads in samples. Grouping of co-expressed genes was performed by using Principle Coordinates Analysis (PCoA) algorithm implemented in the program PAleontological STatistics (PAST) 4.02 3 (Hammer et al., 2001). Downstream networks of regulated genes were constructed using the Web based tool PheNetic (De Maeyer et al., 2013) based on the regulation network available from the PheNetic Web site, which was designed for the strain S. aureus XN108 (CP007447.1) (Zhang et al., 2014). Pairs of homologous genes in the genomes S. aureus BAA-39 and E. coli BAA196 were identified using the program GET_HOMOLOGUES (Contreras-Moreira and Vinuesa, 2013).

Profiling of Epigenetic Modifications
PacBio sequencing allows the prediction of epigenetically modified bases by using the base call kinetics analysis, a useful technique for base modification profiling of bacterial genomes (Liu et al., 2019). Processing of modified bases by SMRT sequencing technologies requires a longer time than the average. The time of base calling may vary sporadically in a rather wide range. The program accounts for repeated base call delays termed interpulse duration (IPD) in overlapping reads to score the likelihood of the epigenetic modification at a given site that is estimated as base modification (BM) score. Tools available in the SMRT Link v5.0.1 software were used with an in-house Python script to generate a pipeline to perform base call kinetic analysis on the PacBio reads generated from chromosomal DNA (Joubert et al., 2019). The pipeline consists of the following steps: (i) the complete genome consensus sequences in FASTA format, obtained by assembling PacBio reads, were indexed by the program samtools to be used as the reference sequence for PacBio read alignment; (ii) PacBio reads were converted from the original BAX.H5 format to BAM format by the tool bax2bam; (iii) reads stored in BAM files were aligned against the indexed reference sequence by the tool blasr (Chaisson and Tesler, 2012); (iv) aligned reads in BAM format were sorted by locations and indexed by samtools sort and index functions; and (v) sorted and indexed BAM files were analyzed by the tool ipdSummary to evaluate the base call kinetics for every nucleotide in the reference genome (the output file * _kinetics.csv). The program stores all the estimated parameters together with context sequences into an output file * _basemods.gff. In step (vi), contextual motifs of base modifications were searched by the tool motifMaker. Thereafter, the epigenetic profiles of the studied genomes were visualized using an in-house Python script, which uses the * _kinetics.csv and * _basemods.gff output files. This in-house pipeline was created to allow a detailed analysis of intermediate files.
To check the correctness of predictions by the in-house pipeline, the standard SMRT Link DNA modification prediction protocol, ds_modification_motif_analysis described in the paper by O'Dell et al. (2018), was used also.
In addition to scoring each nucleotide, the program analyses the context information around the modified bases and identifies contextual motifs. The most common type of epigenetic modification is the methylation of adenine (m6A) and cytosine (m4C) residues. However, this analysis identifies many other nucleotides delaying SMRT sequencing owing to modifications 3 http://folk.uio.no/ohammer/past/ of unknown nature. Thereafter, in the text, these unknown modifications are denoted modA, modG, modC and modT.

Statistical Analysis
All measurements were performed at least 3 times. Average values and standard deviations were used for comparison. To evaluate patterns of gene expression at an experimental condition compared to the control condition, RNA samples were isolated in three repetitions. DESeq2 algorithm estimates fold change and p-values for every gene. Genes with two-fold or above differential expression and p-value smaller than 0.05 were considered as regulated.
SMRT Link v5.0.1 DNA modification pipeline identifies nucleotide positions with a significant base call delay repeated in multiple overlapped reads. The program calculates several statistical parameters such as the interpulse duration (IPD) ratio of the average base call time to the expectation; and the base modification (BM) scores. BM scores are phred-transformed p-values to validate that a kinetic deviation existing at this base position is statistically reliable. BM scores 14 and 21 correspond respectively to p-values 0.05 and 0.01. Only those bases that gained BM scores above 20 in all three repetitions were taken into consideration. As in this study, NC and FS variants were sequenced in three repetitions, an additional level of validation of epigenetic modifications was implemented using the inhouse Python 2.7 script, which selected only those modifications predicted in all three sets of PacBio reads. Average BM scores were calculated for the selected sites.
Pearson correlation coefficients (C p ) of gene co-expression were calculated by Eq. 1 implemented in an in-house Python 2.7 script.
where x i and y i are fold-change values estimated for every i's homologous gene, N -total number of compared homologous genes shared by two genomes.

Antibiotic Resistance Reversion
Antibiotic resistance reversion induced by FS-1 was evaluated after cultivation of the multidrug-resistant bacterium for 10 days on the medium containing a sub-lethal 1/2 MBC of FS-1 with daily re-inoculation into fresh medium. The negative control cultures were cultivated in the same medium without FS-1. The scheme of the experiment is shown in Figure 1. Susceptibility of S. aureus to antibiotics before the experiment and after cultivation at experimental (FS) and negative control (NC) conditions was detected by measuring the growth inhibition zones around antibiotic-containing disks in three repetitions ( Table 1). The initial S. aureus BAA-39 strain was resistant to many tested antibiotics, except for imipenem and cefamandole ( Table 1). The cultivation of multidrug-resistant bacteria on the medium without antibiotics was favorable for selecting more sensitive variants. As a result, after 10 passages on the regular medium, the culture gained susceptibility to methicillin and oxacillin. Rapid decline of bacterial drug resistance on the antibiotic-free medium due to the fitness cost associated with drug resistance mechanisms was described previously. Mechanisms of drug resistance decline involves accumulation of mutations reducing the fitness cost, which, due to absence of antibiotics in the medium, reduced drug resistance (Dunai et al., 2019). Other authors suggested epistatic mechanisms of regulation of the level of antibiotic resistance and the associated fitness cost (Durão et al., 2018). The significant methicillin resistance reduction observed in the culture NC may be explained also by excision of staphylococcal methicillinresistance chromosome cassettes (SCCmec) from chromosomes that occurs naturally on the methicillin-free media with a rate from 10 −8 to 10 −5 depending on the tested strain and the presence of other genomic islands (Stojanov et al., 2015;Almebairik et al., 2020). Cultivation with 1/2 MBC FS-1 induced antibiotic resistance reversion in the treated bacterium compared to the initial strain and the variant NC. After 10 passages on the medium with FS-1, the initially resistant S. aureus BAA-39 strain became susceptible to all the tested antibiotics except for cefamandole ( Table 1).

Complete Genome Assembly
DNA reads, generated from genomes of the NC and FS variants of S. aureus ATCC BAA-39, were assembled into single contig sequences of bacterial chromosomes with no gaps or ambiguities.
No plasmids were identified in these genomes. The lengths of the NC and FS genomes were 2,791,218 bp and 2,792,888 bp, respectively. The completeness of the final assemblies was considered satisfactory, as revealed by the BUSCO software. The obtained genomes were shorter than the previous assembly, GCA_000146385.1, which was 2,865,318 bp, excluding gaps between 83 contigs. The difference in the genome lengths could be due to selection for reduced genome size during the cultivation of the strain in laboratory conditions starting from the date of its isolation in 2010, or due to redundancy of multiple contigs of the previous assembly. The genomes obtained were annotated using the RAST Server and published at GenBank NCBI with the accession numbers CP033505.1 (assembly GCA_003827735.1) and CP033506.1 (assembly GCA_003827835.1) for the NC and FS variants, respectively. The staphylococcal methicillinresistance chromosome cassette (SCCmec), comprising mecA and mecR genes (Holden et al., 2004;Jani et al., 2017), was found in both genomes ∼50 kbp downstream from the replication origin. However, the coverage of PacBio reads overlapping the SCCmec region decreased significantly in both NC and FS genomes, indicating that this mobile element could be unstable in the culture and that there was a tendency to drop this element when S. aureus was cultivated without the selective pressure of the antibiotic that is in agreement with previously published reports (Stojanov et al., 2015;Almebairik et al., 2020). This may explain the fact that both NC and FS cultures gained sensitivity to methicillin after 10 passages of cultivation on the media without the antibiotic (Table 1) Grouping of genes of S. aureus ATCC BAA-39 by the Principle Coordinates Analysis (PCoA) algorithm based on their fold change expression data at 5 tested conditions is shown in Figure 2. Clusters of co-regulated genes are organized in tree branches linked by minimal span tree graphs and distributed between sectors A, B, C and D created by principle coordinates 1 and 2. The gene nodes lying outside of the 95% confidence ellipse are characterized by statistically reliable co-expression patterns. These genes show counter-regulation in opposite sectors. For example, the genes of sector A were characterized by differential levels of expression NC − < FS + < FS − < NC + , while the expression levels of the genes of sector D was NC + < FS − < FS + < NC − . These genes were strongly up-or down-regulated in the culture NC in response to the addition of FS-1 to the medium. In the culture FS adapted to the presence of FS-1 in the medium, these genes were continuedly expressed on an average level and did not significantly respond to the addition of FS-1. Unexpectedly, the expression of these genes in the culture FS was stronger when it was transferred to the medium without FS-1; however, the gene expression variations between the conditions FS − and FS + in many cases were statistically unreliable (p-values > 0.05). The most regulated genes of sector A were tenA aminopyrimidine aminohydrolase of the thiamine salvage pathway; genes of the thiamine biosynthetic pathway thiDEM; potassium-transporting ATPase kdbAB; extracellular polysaccharide biosynthetic genes cap5B and supH, and riboflavin biosynthetic genes ribE, ribH and ribD. The intensively regulated genes of sector D, which were strongly inhibited by FS-1, were Larginine biosynthetic genes argG and argH; ald2-ilvA1-steT-norB operon associated with efflux pump synthesis and activation; lactose and maltose uptake and utilization genes mtlA and ldhA; alpha-phosphotrehalase treC; PTS threhalose transporter treP and trehalose operon transcriptional repressor treR.
The expression pattern of the genes of sector B is characterized by the following order of expression levels: FS − < NC − < NC + < FS + ; and the genes of sector C were regulated oppositely: FS + < NC + < NC − < FS − . These genes were sensitized by 10 passages of cultivation with FS-1 and responded to the addition of FS-1 by a strong positive or negative regulation. In the culture NC, the differential expression of these genes at the conditions NC − and NC + in many cases was statistically unreliable. One of the most regulated genes of sector B is carbamoylphosphate synthase carB. Carbamoyl-phosphate synthesized by this enzyme may be further used either for biosynthesis of L-arginine or pyrimidines. Because the next steps of Larginine biosynthesis by argininosuccinate synthase argG and argininosuccinate lyase argH were strongly downregulated at this condition, it may be assumed that the activation of carB was dictated by the need for DNA synthesis or repair. The activation of gene carB was co-regulated at these conditions with upregulation of protein and nucleotide deglycase hchA. Other co-regulated genes of sector B were genes for DNA lesion repair proteins (uvrABC); chaperon dnaJ; non-heme ferritin ftnA; spermidine/putrescine import protein potA; cobalamin biosynthesis ferrochelatase cbiX; respiratory nitrate reductase narIJ; iron/manganese ABC transporter sitBC; suppressor of fatty acid biosynthesis fapR; mercury reductase merA2; components of glucose and glucosamine PTS uptake systems ptsG and glvC; and ncRNA teg25as (HMPRNC0000_2674).
Oppositely regulated genes of sector C include biotin carboxylase accC; cytochrome d ubiquinol oxidase subunit appC; nickel/cobalt efflux protein nixA; peptidoglycan polymerase mgt; betaine biosynthetic gene betB; glycine/betaine transporter opuD; ribulosamine/erythrulosamine 3-kinase yniA involved in protein deglycation and 10 genes encoding regulatory ncRNA. Half of these ncRNA genes, HMPRNC0000_0900, 1335, 2206 (RNAIII), 2747 and 2902, were expressed only in the culture FS when it was cultivated on the normal medium without FS-1. Other five genes, HMPRNC0000_0052, 1456, 1494 (sRNA233), 1661 and 2390, were expressed also in the culture NC on the medium without FS-1. Addition of FS-1 to the medium completely abolished or significantly inhibited the expression of these regulatory elements. The presence of FS-1 in the medium halted the expression of three protein coding genes: ferrous iron transporter HMPREF0783_0807, periplasmic oligopeptide-binding protein oppA and ATP-binding subunit clpL of the Clp protease.
The effect of FS-1 on expression of several important genes involved in cell division and replication, and the genes of the methicillin resistance cassette (SCCmec) is shown in Table 2. As expected, the addition of FS-1 to the culture NC does inhibit the cell growth/division associated genes and antibiotic resistance genes. An unexpected observation was that the incubation of FS-1 adapted culture FS on the medium without the drug FS-1 caused an even stronger suppression of bacterial growth and the expression of the antibiotic resistance genes. Comparison of the gene expression in pairs of conditions, NC − /NC + and FS − /FS + , confirmed that the removal of FS-1 from the medium stressed culture FS in a very much similar way as the addition of FS-1 to the medium stressed the culture NC. In total, 127 genes were regulated at both experimental conditions, NC − /NC + and FS − /FS + , and 85 of these genes (67%) were counter-regulated meaning that the addition of FS-1 to the culture NC caused the same effect on these genes as the removal of FS-1 affected these genes in the culture FS. Such gene expression regulation may be associated with the response to the medium replacement stress.
Remaining 42 genes, which were co-regulated at the experimental conditions NC − /NC + and FS − /FS + , were of interest as they represented an effect caused by FS-1 independently whether the bacterial culture was adapted to the presence of the drug or not. There are 27 genes activated by FS-1 in both cultures, NC and FS, which include genes of Nar respiratory nitrate reductase operon; genes involved in biosynthesis of riboflavin (rib), molybdopterin (moaC), cobalamin (cbiX) and extracellular polysaccharides (supH); α-amylase malA; glucose and maltose uptake genes ptsG, ptbA and glvC; ribonucleoside diphosphate reductase subunits nrdIE; flavin-utilizing monooxygenase yhbW; oxidoreductase nfrA; thiol peroxidase tpx; mercuric reductase merA2; ribosomal protein S4 and SSU maturase ylqF important for translational accuracy and several other metabolic enzymes. There were 15 genes strongly inhibited by FS-1 in both cultures, NC and FS. They include de novo purine biosynthesis repressor purR; sugar phosphate antiporter uhpT; NADH-ubiquinone oxidoreductase mnhA2; alkaline shock gene asp; transmembrane transporter gudP and transporter of small amino acids and other nitrogen containing compounds opuD; sucrose uptake gene scrA and Clp protease gene involved in thermotolerance and intracellular multiplication.  counter-regulated genes ( Figure 3A). However, more significant positive correlation, 0.231, was obtained when the gene expression pattern in E. coli at condition NC − /FS + was compared to the pattern of gene expression in S. aureus at condition FS − /FS + (Figure 3B). Gene expression comparison at different conditions allowed the identification of commonly regulated genes in these distant model microorganisms under the effect of FS-1. Only eight protein coding genes were positively regulated in FS + compared to NC − in both organisms ( Figure 3A). These genes were heavy metal ion exporting P-type ATPase copA (Sitthisak et al., 2007), kdpBC K + -transporting ATPase subunits and DNA-binding response regulator kdpE expressed in bacteria in response to the osmotic stress (Epstein, 1986;Gründling, 2013); chaperons dnaJ and dnaK; capsular biosynthesis gene epsD and alkyl hydroperoxide reductase ahpF that protects cells against DNA damage by alkyl hydroperoxides. FS − /FS + comparison confirmed that genes copA and ahpF are positively regulated by FS-1. Osmotic stress response genes kdpBC and epsD appeared counter-regulated at this condition as they were stronger activated in the S. aureus culture FS when it was transferred to the medium without the drug FS-1. It may be concluded that the presence of FS-1 in the medium and removal of this drug after culture adaptation to the presence of the drug, both invoke an osmotic stress to the bacterial culture. Genes betAB catalyzing biosynthesis of osmoprotectant glycine betaine from choline (Falkenberg and Strøm, 1990) were strongly downregulated by FS-1 in both microorganisms. The reason for this inhibition of betaine synthesis is unknown but may be associated with complexation of glycine-betaine with iodine (Blackwell et al., 2001).

Comparison of Gene Expression
Genes, which were co-regulated by FS-1 in S. aureus and E. coli, are summarized in Figure 4 by a PheNetic network of functionally associated genes. In both genomes, the presence FIGURE 4 | PheNetic networks of regulated genes grouped by functional or regulatory interaction between these genes or their products using the downstream regulation algorithm of the program PheNetic. Important functional groups of genes are labeled. Upregulated genes are depicted by yellow nodes and downregulated genes -by blue nodes (vertices). Color intensity indicates the level of transcriptional regulation (fold-change values). Gray nodes are genes or transcriptional regulators involved in the network, which expression was not reliably changed at the aforementioned conditions. Green edges show activation relations; blue edges -inhibitory relations; and brown edges -ambivalent or neutral relations. Direct regulation by transcriptional regulators are indicated by arrowheads.
of FS-1 in the medium had a strong inhibitory effect on enzymes of tricarboxylic acid cycle (TCA), sucABCD; fatty acids β-oxidation pathway (acetyl-CoA acetyltransferase atoB) and glyoxylate pathway associated with the fatty acid catabolism (citrate synthase gltA). In contrast, fatty acid biosynthesis was activated by upregulation of acetyl-CoA carboxyltransferase accD. The common response of S. aureus and E. coli to the presence of FS-1 in the medium was a switch to anaerobiosis that could be associated with the need to reduce the oxidative stress caused by iodine. Another hypothesis was that iodine may damage cytochrome molecules. Indeed, the transcription of all cytochrome o ubiquinol oxidase subunits, cyoABCDE, was brought to an almost complete stop in E. coli; however, in S. aureus these genes were constantly expressed. TCA cycle enzymes were strongly inhibited in S. aureus while fermentation genes, alpha-acetolactate decarboxylase badA, Llactate dehydrogenase ldhA; formate acetyltransferase activating enzyme pflA and formate acetyltransferase pflB, were activated. This switch to the anaerobic lifestyle most likely was regulated by arcR that is the main regulator of anaerobiosis in S. aureus (Kohler et al., 2008). This transcriptional regulator was upregulated in NC and FS cultures when cultivated with FS-1. Fermentation related pathways were activated in E. coli by FS-1 in a similar way. Both these microorganisms used the nitrate/nitrite respiration pathway that was indicated by activation of nitrite reductase nirB and heme biosynthesis gene hemE (in both microorganisms), and the respiratory nitrate reduction chain, narH, narJ, narH and narG, in S. aureus. The analysis of the levels of expression of two NAD-dependent glyceraldehyde-3phosphate dehydrogenases, gapA1 (HMPRNC0000_0860) and gapA2 (HMPRNC0000_1815), the former one is involved in glycolysis and the latter one controls gluconeogenesis in S. aureus (Chaffin et al., 2012), showed that the glycolysis was approximately 14-fold more active than the gluconeogenesis in S. aureus at the given conditions. Both these genes were upregulated in S. aureus FS that reflects a general activation of glycolysis and gluconeogenesis in the culture adapted to the presence of FS-1 in the medium. The transcriptional comparison at the conditions FS − /FS + showed 5-fold upregulation of gapA1 with no significant change in the level of expression of gapA2. Moreover, cggR repressor of gluconeogenesis (Kohler et al., 2008;Chaffin et al., 2012) was activated by FS-1. Such activation of glycolysis is typical for S. aureus growing at anaerobic conditions (Kohler et al., 2008); however, it was opposite to what was observed in E. coli where FS-1 treatment caused upregulation of the genes involved in gluconeogenesis and downregulation of their glycolytic counterparts (Korotetskiy et al., 2020b). It may be explained by the fact that E. coli, in contrast to S. aureus, has several alternative glycolytic pathways bypassing glycolysis, which are Entner Doudoroff, Embden-Meyerhof-Parnas and oxidative pentose phosphate pathways (Hollinshead et al., 2016). The genes of the Entner Doudoroff pathway, phosphogluconate dehydratase edd and keto-hydroxyglutaratealdolase eda, were 1.8 and 1.6-fold upregulated in the FS-1 treated E. coli  Another commonality in the response of S. aureus and E. coli to FS-1 was activation of multiple enzymes involved in the nucleotide and protein biosynthesis, translational accuracy control and DNA repair. The latter group is represented by methyl-directed mismatch repair protein MutS, nucleotide excision repair proteins UvrA and UvrC, and NAD(+)dependent DNA ligase LigA. Two chaperons, DnaJ and DnaK, also were activated in both microorganisms induced by FS-1.
There were only a few genes, which generally were counterregulated by FS-1 in S. aureus and E. coli (Figure 3B). Two genes of the mixed acid fermentation pathway, fumarate reductase frdA and acetate kinase ackA, were strongly upregulated in E. coli but downregulated in S. aureus. Purine salvage hypoxanthineguanine phosphoribosyltransferase hpt was expressed in E. coli and in the S. aureus culture NC, but it was transcriptionally silent in the S. aureus culture FS at all tested conditions. Moreover, this gene was 2.4-fold upregulated in E. coli at condition NC − /FS + , but in S. aureus this gene was 5-fold downregulated at condition NC − /NC + . This opposite response of one of the most important gene of purine turnover in E. coli and S. aureus remained unclear. Other oppositely regulated genes were cyclic pyranopterin phosphate synthase moaA involved in the molybdenum cofactor biosynthesis pathway, ferric uptake regulator perR, NADP-dependent acetaldehyde dehydrogenase aldB and carbonyl stress DNA and protein repair deglycase hchA. All these genes were strongly upregulated by FS-1 in S. aureus and downregulated in E. coli. The gene lmrB2 encoding the multidrug efflux protein was transcribed only in the S. aureus culture NC and FS, but its homolog emrY in E. coli was transcriptionally silent at the given conditions. This gene was constantly expressed in S. aureus under condition NC − /FS + , but it was 3.3-fold upregulated under condition FS − /FS + .

General Patterns of Modified Bases
Base modifications were predicted in three repeats of sequencing of DNA samples generated from NC and FS genomes. Epigenetically modified nucleotides, e.g., methylated, partially oxidized or halogenated residues, are characterized with significantly increased BM scores. However, BM scores of epigenetic modified nucleotides may be lower, if the modification occurred only in a fraction of the bacterial population.
The distribution of BM scores calculated for all bases in NC and FS genomes of S. aureus is shown in Figure 5 by lines representing frequencies of adenine (A), thymine (T), guanine (G) and cytosine (C) residues along the BM score ranks. Axis Y represents decimal logarithms of the numbers of nucleotides per rank. In both genomes, 95% bases had BM scores below 14 and 99% bases had scores below 21. In the genome NC (Figure 5A), numbers of nucleotides with higher BM scores decreased gradually except for adenine residues having one peak in the range of BM scores from 50 to 65, and a larger escalation of adenine residues with BM scores 170-260. The latter peak corresponds to adenine residues N-methylated at sixth carbon atoms (m6A methylation). Similar increase in number of higher scored cytosine residues corresponds to methylation of these residues either at fourth or fifth carbon atoms (m4C and m5C methylation). Higher BM scores determined for guanine and thymine residues, but also for several adenine and cytosine residues, were associated with other types of epigenetic modifications; which cannot be recognized by the current version of the program ipdSummary. They were denoted as modA, modG, modT and modC.
In the genome FS, a drastic increase of all nucleotides with BM scores in the range from 50 to 95 was observed. In Figure 5B, these characteristic peaks are pointed out by an arrow. The area of escalation of high scored nucleotides (BM ranks 170-290) became wider in the genome FS.

Patterns of Adenine Methylation in S. aureus
Adenine methylation and its role in gene expression regulation and adaptation of prokaryotes were discussed in many publications (Casadesús and Low, 2006;Bierne et al., 2012;Sánchez-Romero et al., 2015). Therefore, comparison of methylation patterns was in the focus of this study.
Detailed analysis of context sequences adjacent to methylated adenine residues revealed that the majority of m6A sites were bipartite methylated fragments (BMF) of DNA with the first methylated adenine on the direct strand and the second one on the reverse complement strand with six or eight nucleotides in between ( Figure 6C). According to REBASE classification (Roberts et al., 2015), this is type-I methylation that is often associated with DNA restriction-modification systems. Type-I methyltransferases are frequent in bacterial genomes and usually are comprised of three genes encoding a methyltransferase, motif recognition subunit and cognate restriction endonuclease. The best candidate for this role in the genome of S. aureus BAA-39 is type-I methyltransferase hsdM associated with the specificity unit gene HMPRNC0000_1947; however, there were no associated restriction enzymes around these genes. The addition of FS-1 to the medium has an opposite effect on the expression of these two genes (Supplementary Table 1) that implies their transcriptional control by different promoters. Comparison of the context sequences showed no significant sequence conservation of nucleotides enclosed by the bipartite m6A residues, except that they were AT-rich. Contrary, the triplets upstream and downstream of the m6A sites were semiconserved (shown by capital Ns in Figure 6C). The frequencies of combinations of different triplets at 5 -and 3 -ends of BMF elements are shown in Figures 6A,B. Two types of BMF elements were found. The sequences of type 1 were asymmetric ( Figure 6A). On one side, the flanking triplet was NCC-m6A (or m6A-CCN in the reverse-complement variants), where N is A, T, G or C in order of the frequency of occurrence. The opposite flanking triplets were loosely conserved. They could be CAA, CGA, TTG, ATG, TGG, AGG, AAG, TAG, ACG, CAG, CAG, GTG, CTG, CTG, GAG, TCG, GCG, CGG, GGG, CGG, CCG, or CCG in order of the frequency of occurrence. Spacer sequences were 8 bp long. In total, there were 1,105 BMF elements of this type (82%). Pairs with loosely conserved sequences at the 5 -end and conserved 3 -end, NNNA-(n) 8 -TGGN, were denoted as direct type 1 BMF elements (T1-direct), as they were found predominantly in the leading replichore (the clockwise chromosomal sequence directed from the replication origin to the terminus). The most frequent motif of this type found on the chromosome 54 times was CAAA-(n) 8 -TGGT. The complement type 1 BMF elements, CCNA-(n) 8 -TNNN, were more frequent in the lagging replichore (counter-clockwise sequence from the replication origin to terminus) and thus denoted as T1-complement BMF elements ( Figure 6A). The most frequent motif of this type found on the chromosome 44 times was ACCA-(n) 8 -TTTG. In the BMF motifs shown above, nucleotides in bold are methylated sites, and the nucleotides opposing the bipartite methylation on the complement DNA strand are denoted with an italic font.
BMF elements of type 2 (they also may be classified as REBASE type-I methylation) have equally conserved flanking triplets at the 5 -and 3 -ends. These elements were shorter with 6 bp in the spacer regions. Four groups of these sequences were found. In Figure 6B, these groups are separated by diagonal lines and titled T2-1, -2, -3 and -4. Sequences of the groups T2-1 and T2-4 were defined as direct sequences (T2-direct) and sequences of the groups T2-2 and T2-3 were defined as complement sequences (T2-complement). Direct sequences were more frequent in the leading replichore and complement ones were frequent in the lagging replichore. In total, there were 240 instances of type 2 BMF elements (18% of all instances). The most frequent motif of the direct type is TGAA-(n) 6 -TACA (23 instances), and TGTA-(n) 6 -TTCT (15 instances) is the most common motif of the complement type. The disproportional density of the direct and complement BMF elements found on the leading and lagging replichores are shown in Table 3 and illustrated in Figure 7A. This methylation should be controlled by an alternative methyltransferase. The most likely candidate is the orphan (not associated with any restriction enzymes) adeninespecific methyltransferase HMPRNC0000_1848 (upregulated in the medium with FS-1).
A common property of BMF elements of both types, direct and complement, is that the upstream m6A residues are almost always located on the direct DNA strand (clockwise from the replication origin), and the downstream counterpart m6A residues are on the reverse complement strand. Of 1,355 BMF elements found in NC and FS genomes, only two in NC and four in FS were in reversed order, with the upstream m6A on the reverse-complement strand. All these facts demonstrate that the location and orientation of BMF elements on the chromosome of S. aureus is of importance for this microorganism.
The total number of BMF elements in NC and FS genomes remained the same: 1,355. The complete list of identified BMF elements with respective locations, coverage and BM scores is in Supplementary Table 2. Among them, 1,321 BMF elements were at the same locations on the chromosomes NC and FS, but in 34 cases the BMF elements in NC were replaced with alternative elements in FS. Ten such replacements were in the body of the methicillin resistance cassette SCCmec. BMF were infrequent in SCCmec compared to other parts of the chromosome including identified prophages. It may be concluded that the SCCmec insert was subjected to alternative methylation depending on the growth condition and stresses experienced by the microorganism.
A striking difference in profiles of the m6A methylated sites in the NC and FS variants of S. aureus was in the distribution of BM scores of the m6A nucleotides along the chromosome ( Figure 7B). In NC, the BM scores were significantly higher near the chromosomal origin of replication, but they decreased progressively in the direction toward the replication terminus on both replichores. The culture cultivated with FS-1 showed a different distribution of BM scores along the chromosome with a smaller difference between the BM values of the m6A sites around the replication origin and the terminus. BM scores indicate either the efficacy of methylation of the site (proportion of methylated nucleotides at this location in the population) or the local coverage. Analysis of coverage values of PacBio reads aligned against the genome consensus sequences showed the same parabolic distribution with lower values around the replication terminus (Supplementary Figure 1). Due to the continuous replication, the areas around the replication origin in bacterial cells are semi-diploid that increases the coverage of randomly generated reads compared to the area around the replication terminus (Marczynski and Shapiro, 1993). It may be assumed that in bacterial populations with different distribution of dividing and stationary cells, the shape of the coverage curve will be different. The treatment of S. aureus with FS-1 caused an obvious phase shift in the population illustrated in Figure 7B by the average BM score lines. Strangely, the treatment of E. coli BAA-196 cells with FS-1 affected the population oppositely by making the coverage/BM score curves more concave in the FS-1 treated culture ( Figure 7C). The chromosome of the latter organism also was rich with m6A methylated sites; however, the major methylation motif was typical for Enterobacteriaceae REBASE type-II Dam methylation GATC (Marinus and Løbner-Olesen, 2014;Blow et al., 2016). Only a small fraction of bipartite m6A methylated sites were associated with type-I methylation motifs AAC-(n) 6 -GTGC and GCAC-(n) 6 -GTT, which resemble the corresponding BMF elements of the genome of S. aureus BAA-39.

Patterns of Other Epigenetically Modified Nucleotides
Methylated cytosine residues identified by the program ipdSummary as m4C methylated sited were infrequent and they were not associated with any sequence motifs. However, several m4C sites showed high BM scores in all three PacBio SMRT sequencing repeats. In contrast to the m6A methylated sites discussed above with preserved locations in both the genomes, the m4C methylated nucleotides, even those with the highest BM scores, were distributed differently in the genomes NC and FS. Table 4 shows the distribution of the high scored m4C sites with BM scores above 160. This cutoff value was selected by analyzing the BM score distribution shown in Figure 5. Only 5 high scored m4C sites out of 15 in NC and 19 in FS genomes were located at the same positions in both genomes: two sites in the genes hyaluronate lyase precursor hysA and DNA-methyltransferase hsdM, and other three sites in non-coding regions. Other genome specific high scored m4C sites were found in several coding sequences and frequently in transcribed regions not associated with protein coding genes, which could be unidentified genes for regulatory RNA molecules. It may be assumed that the alternative m4C methylation potentially can interfere with these regulatory elements. Two FS specific high scored m4C sites were found in the methicillin resistance cassette SCCmec. No cytosine specific methyltransferases were found in the genome of S. aureus BAA-39 except for two uncharacterized orphan methyltransferases HMPRNC0000_1724 and HMPRNC0000_1746 (both are significantly down-regulated at the FS − condition). There was an increase of modified nucleotides of all types with moderate BM scores in the range from 70 to 95 in the FS-1 treated S. aureus (Figure 5). The distribution of these modified nucleotides was further analyzed as shown in Figure 8. An increased density of modified nucleotides was observed in the second half of the lagging replichore (downstream of the replication terminus) and in the first half of the leading replichore (downstream of the replication origin). This disproportional distribution of the modified bases was observed in both genomes; however, there are more modified bases in the genome FS compared to NC.
A possible effect of base modifications on the gene expression level was analyzed. Taking the average level of gene expression in the culture FS as 1, the level of expression of genes with additional modG modifications compared to NC was 0.732 on average, while in those genes with unchanged number of modG, the expression level was 1.134 on average. There was a clear tendency of avoidance or counter-selection of modG modifications from highly expressed genes. However, it should be admitted that because of the significant variability in gene expression values within the groups, the difference between average values was statistically unreliable, meaning that there were many other factors affecting gene expression besides the modG modifications. When the numbers of modA, modG, motT and modC modifications were considered, the difference in the level of gene expression was smoothed out to 0.836 for the genes with novel nucleotide modifications and 1.096 for the unmodified genes. It shows that the other types of modifications did not contribute to the effect of avoidance of highly expressed genes, or their effects were ambivalent.

Complexing of FS-1 With Bacterial DNA
Epigenetic changes in bacterial chromosomes under the effect of the iodine-containing nano-molecular complex FS-1 suggested possible complexing of DNA with nano-micelles and/or with the released iodine, which can potentially halogenate chromosomal nucleotides, predominantly thymine residues (Ilin et al., 2013).
To check this hypothesis, complexing of FS-1 with the chromosomal DNA was investigated experimentally.
To study the ability of FS-1 to penetrate bacterial cell wall and create complexes with the chromosomal DNA, the drug was synthesized with the radioactive isotope 131 I (20 MBq/ml). The residual radioactivity was measured in DNA samples extracted from the treated-and-washed cells. It was found in three repetitions that the radioactivity of the extracted DNA constituted 12.76% ± 9.04 of the residual radioactivity of the treated cells that corresponded to 0.14 ± 0.015 Bq/ng of the extracted DNA. The residual radioactivity can be associated with both complexing of FS-1 with DNA and with a direct halogenation of nucleotides by iodine isotopes. In a similar study with E. coli ATCC BAA-196, the residual radioactivity of the extracted DNA constituted only 0.46% ± 0.15 of the residual radioactivity of the treated bacterial cells measured after washing the cells (Korotetskiy et al., 2020b).

DISCUSSION
Drug-induced reversion of antibiotic-resistant bacterial pathogens into antibiotic-susceptible phenotypes is a promising way to control outbreaks of drug-resistant infections. There are two basic approaches to antibiotic resistance reversion: (i) direct inhibition of the resistance mechanisms; and (ii) active selection against drug resistance in bacterial populations using evolutionary and physiological interactions between drugs and bacteria. The latter approach was discussed and modeled in the review by Baym et al. (2016); however, these authors did not gain deep insight into possible molecular mechanisms of counter-selecting against drug resistance in bacterial populations because of absence of workable models. The introduction of FS-1 into clinical practice (Ilin and Kulmanov, 2014;Ilin et al., 2017;Korotetskiy et al., 2017;Islamov et al., 2018;van Niekerk et al., 2018) allowed detail investigation of this phenomenon. A series of experiments on selected multidrug resistant model microorganisms was performed, first on E. coli ATCC BAA-196 (Korotetskiy et al., 2020b) and currently on S. aureus ATCC BAA-39. FS-1 steadily converts the treated bacteria to antibioticsensitive phenotype. Bacteria remain sensitive even after the removal of FS-1 from the medium at least for several generations. While the FS-1 treated bacterial populations regain the antibiotic resistance at in vitro experimental conditions under the selective pressure of antibiotics (data on the dynamics of antibiotic resistance regaining were not shown due to significant variations between the repeated experiments), the in vivo experiment published before showed a promising therapeutic effect of co-administration of FS-1 with antibiotics (Ilin et al., 2017). This drug can be applied prior to antibiotics or in combination with them. Supplementary drug may induce phase variations in bacterial populations and favor the selection of drug-sensitive variants. Mechanisms of phase variations can be genetic, such as mutations in genes and promoter regions (Henderson et al., 1999;Gründling, 2013), or epigenetic, including altering methylation/modification patterns (Casadesús and Low, 2006;Low and Casadesús, 2008;Phillips et al., 2019a), nucleotide oxidation and DNA phosphorothioate modifications (Liu et al., 2019). Recent publications demonstrated the important role of phase variations in host-pathogen and host-environment interactions, and in bacterial adaptation (Haagmans and van der Woude, 2000;Hernday et al., 2002;van der Woude and Bäumler, 2004;Phillips et al., 2019a,b). The role of epigenetic modifications in antibiotic resistance development has also been demonstrated (Nyce et al., 1993;Adam et al., 2008). Nevertheless, the possibility to use induced phase variations to revert antibiotic resistance has never been explored. In this work, a third-generation sequencing approach was used to investigate possible genetic and epigenetic mechanisms of antibiotic resistance reversion induced by FS-1 in S. aureus ATCC BAA-39.
In the conducted experiment, S. aureus culture was cultivated in test-tubes for 10 days with daily re-inoculations to fresh MH medium containing 1/2 MBC FS-1 (450 µg/ml). This culture was denoted FS. Negative control (NC) culture was cultivated in parallel on the MH medium without FS-1.
An expectation was that antibiotic resistance reversion may be induced by altering DNA methylation profiles. Methyltransferases identified in the genome of S. aureus BAA-39 were not associated with DNA restriction-modification systems and played other roles rather than the protection of DNA against cognate restriction enzymes. It was hypothesized that type-I methylation in this genome may be associated with architecture-imparting sequences (AIMS). The m6A methylation pattern of S. aureus BAA-39 showed well-organized global replichor-oriented distribution of BMF elements that resembles the distribution of AIMS on bacterial chromosomes (Hendrickson et al., 2018). AIMS are semi-conserved 6-8 bp sequences preferentially abundant on the leading strand. AIMS may be synonyms of FtsK orienting polar sequences (KOPS) assisting the directional loading of the FtsK translocases (Bigot et al., 2006), which are important for proper chromosomal replication and DNA repair. Spacer regions of BMF elements in S. aureus, as well as in E. coli, were 6-8 bp long, the same as AIMS. Further study should be performed to investigate possible relationships between AIMS, KOPS and BMF elements as architectural units of bacterial chromosomes. The absence of AIMS in horizontally acquired genomic islands suggests a possibility that these sequences may help to prevent insertion of mobile genetic elements by allowing a higher frequency of mutations in these inserts (Hendrickson et al., 2018). The global m6A methylation pattern persisted the treatment with FS-1 except for the SCCmec genomic island. The treatment with FS-1 increased instability of this regions and possibly promoted its excision from the chromosome of S. aureus due to improper methylation of internal BMF elements. Indeed, a drop of coverage of PacBio DNA reads generated from this region was observed in both, NC and FS cultures, suggesting an excision of the genomic islands at least in a fraction of cells; however, the coverage difference between NC and FS genomes was insignificant. A similar destabilization of the virulence plasmid in the FS-1 treated E. coli BAA-196 was reported in a previous publication (Korotetskiy et al., 2020b). The varying rate of excision of SCCmec cassettes from Staphylococcal populations was reported in another publication (Almebairik et al., 2020). However, it should be noted that in the current work, the rate of SCCmec excision in NC and FS cultures was not measured. The hypothesis of destabilization of antibiotic resistance genetic inserts by the treatment with FS-1 must be checked in future studies.
The altered adenine and cytosine methylation landscape in the FS-1 treated S. aureus was clearly demonstrated; however, it should be accepted that the role of the differential nucleotide methylation must be proved in further studies as the altered patterns of methylation equally likely may be either gene regulation factors or the result of differential gene regulation at different growth conditions. Nevertheless, the repeated culturespecific methylation of the NC and FS genomes of S. aureus confirmed in three repetitions of PacBio genome sequencing demonstrates a biological importance of this phenomenon and possibility to use the methylation landscape as a genetic marker of the initial and FS-1 treated culture variants.
The treatment of bacterial cultures with FS-1 caused an increase in numbers of modified nucleotides in the chromosome of S. aureus (Figures 5, 8). These global epigenetic modifications, which chemical nature remained disputable, could be an integral part of induction of antibiotic resistance reversion. Relatively low BM scores (65-80) of identified modified sites and the absence of contextual motifs imply a random pattern of nucleotide modifications caused by FS-1. Analysis of the distribution of modG sites showed a tendency to avoid highly expressed genes that may indirectly indicate an interference of the modG modifications with the gene transcription or an increased rate of mutation of the affected genes. The modG modifications may correspond to O-6-guanine methylation resulted from an abnormal activity of DNA methylases. O-6-methylguanine is a highly mutagenic adduct provoking single nucleotide mismatches during replication (Gu et al., 2017). Guanine residues may be oxidized to 7,8-dihydro-8-oxoguanine that is another strong mutagen due to its complementation with both cytosine and adenine (Neeley and Essigmann, 2006). It may be expected that the rate of mutations should increase in the FS-1 treated S. aureus. Frequent re-inoculations and the growth on the simple nutrient rich medium were favorable for accumulation of genetic mutations in genes, which may be out of importance at this laboratory condition. In total, 104 protein coding genes in NC and 271 genes in FS were truncated owing to frameshift mutations. The distribution of truncated genes in both genomes identified by comparison to the initial whole genome sequence of this microorganism (AEEK00000000.1, published in 2013), is shown in Figure 9. Particularly, gene mecA encoding for an alternative penicillin-binding methicillin-resistant protein PBP2a has a frameshift mutation at its 5 -end in the genome FS that makes the encoded protein shorter by 61 amino acids. Alignment of Ion Torrent reads generated from the total RNA samples showed that the mutated gene was transcribed possibly from an alternative promoter. The level of expression of this gene was sixfold downregulated in FS compared to NC grown on the medium without the drug FS-1. It may be a result of either this mutation, or an alternative gene regulation, or excision of the genomic island from the chromosome. The neighbor gene mecR encoding for methicillin resistance regulatory sensor-transducer was 3-fold down-regulated in FS at this condition ( Table 2).
The gene expression patterns in S. aureus and E. coli induced by the treatment with FS-1 showed some level of similarity (Figure 3). While the calculated Pearson correlation coefficient, 0.231, is moderate, it could not be higher as the homologous genes of these phylogenetically distant microorganisms were expressed differently even at the same condition (data not shown). An interesting finding was that the general pattern of gene expression in the negative control E. coli was more similar to that in the FS-1 treated S. aureus rather than to NC culture ( Figure 3). Nevertheless, there was a commonality in the general trend of gene regulation in both microorganisms (Figure 4): strong down-regulation of the TCA pathway; switch from aerobic respiration to anaerobic respiration or fermentation; increased fatty acid synthesis with prohibited β-oxidation of fatty acids; activated synthesis of proteins and nucleotides; activation of genes associated with osmotic stress response but strong downregulation of synthesis of betaine osmoprotectants; and selective activation or inhibition of transmembrane transporters possibly caused by an attempt to diminish the uptake of iodinecontaining compounds. Noteworthily, while S. aureus may survive anaerobic or microaerophilic conditions for long time, anaerobiosis affects strictly the virulence and antibiotic resistance potential of this bacterium (Fuchs et al., 2007). The transition to anaerobiosis caused by the treatment with FS-1 may explain the observed increase in susceptibility to antibiotics ( Table 1).
One striking difference was the strong upregulation of nucleotide deglycase hchA in the FS-1 treated S. aureus. HchA is a component of the major nucleotide glycation repair system activated under carbonyl stress that reverses methylglyoxal and glyoxal damage mostly in guanine residues via nucleotide sanitization and direct nucleic acid repair (Richarme et al., 2017). Carbonyl stress in bacteria is associated with strong oxidative stress and intracellular accumulation of various reactive carbonyl compounds, which may be halogens and their products (Semchyshyn and Lushchak, 2012;Kosmachevskaya et al., 2015;Ezraty et al., 2017). Nitrogen and oxygen atoms in purine and pyrimidine bases are the most attractive targets for reactive carbonyl compounds, which lead to their chemical modification resulting in cytotoxicity and mutagenicity Frontiers in Microbiology | www.frontiersin.org (Ellis, 2007;Liu et al., 2010;Semchyshyn and Lushchak, 2012). In E. coli treated with FS-1, the homologous gene hchA was downregulated. It suggests that the transmembrane flux of iodine-bound compounds into FS-1 treated E. coli cells was significantly lower even that in both cases FS-1 was applied in 1/2 MBC concentrations calculated for respective microorganisms (450 µg/ml for S. aureus and 500 µg/ml for E. coli). In contrast to S. aureus, no increase in the number of modified nucleotides and no increase in the rate of frameshift mutations was observed in E. coli treated with FS-1 (Korotetskiy et al., 2020b). In concordance with this was the detection of only traces of the radioactive iodine in the DNA extracted from the cells of E. coli treated with 131 I marked FS-1, while in the respective DNA sample extracted from the treated S. aureus cells, the residual radioactivity of DNA constituted 12.76% ± 9.04 of the total radioactivity of washed bacterial cells. One commonality in the response of S. aureus and E. coli to the treatment with FS-1 in terms of epigenetic modifications lies in an alternative distribution of high scored methylated cytosine residues that could be a common mechanism of population phase variation in response to environmental stresses.
A recent study confirmed an activation of genes associated with the oxidative stress response in S. aureus in 5 min after injection of FS-1 into growth medium . Changing of the redox potential in proximity with the chromosomal DNA may promote the oxidation of nucleotides. Moreover, close co-location of two oxidized guanine residues on the chromosome can lead to double-strand DNA breaks and cell death (Foti et al., 2012). These authors claimed also that closely spaced 8-oxo-deoxyguanine lesions leading to lethal double-strand DNA breaks and RNA mistranslation contribute significantly to the cytotoxicity of many broadly used antibiotics including beta-lactams, quinolones and aminoglycosides. In this case, the DNA oxidation by FS-1, which may be inherited by bacterial daughter cells in several generations, may increase the sensitivity to antibiotics even when FS-1 is removed from the medium.

CONCLUSION
This study demonstrated reversion of antibiotic resistance obtained on the reference multidrug resistant strain S. aureus ATCC BAA-39 treated with the iodine-containing nano-micelle drug FS-1. Comparison of complete genome sequences, transcriptional profiles and landscapes of epigenetic modifications suggested several possible molecular mechanisms of the antibiotic resistance reversion, which include (i) profound transition of bacterial metabolism toward anaerobiosis due to damaging of oxygen-depended terminal cytochrome:ubiquinol electron transfer complexes by iodine that eventually caused unfavorable conditions for bacteria to maintain the initial level of antibiotic resistance; (ii) increased rate of mutations in the FS-1 treated cells due to oxidation and/or halogenation of chromosomal DNA nucleotides that eventually may affect many important cellular systems including antibiotic resistance genes; (iii) osmotic, oxidative and carbonyl stresses experienced by the FS-1 treated cells may aggravate the antibiotic resistance fitness cost; (iv) the treatment with FS-1 may cause destabilization and loss of horizontally acquired antibiotic resistance genomic islands (e.g., SCCmec) and virulence plasmids that was reported for the FS-1 treated E. coli BAA-196 (Korotetskiy et al., 2020b).
In total, the expression of 2,952 protein or regulatory RNA coding genes of S. aureus BAA-39 participating basically in all biological processes were affected during 10 passages of the cultivation with FS-1 (Supplementary Table S1). Further study will aim at identification and experimental validation of common mechanisms of action of iodine-containing nanomicelles on taxonomically distant microorganisms, S. aureus, E. coli (Korotetskiy et al., 2019a(Korotetskiy et al., , 2020b, A. baumannii (Korotetskiy et al., 2020a) and M. tuberculosis (Ilin et al., 2015(Ilin et al., , 2017. The drug FS-1 has successful passed clinical trials in Kazakhstan as a supplement of the antibiotic therapy against multidrug resistant tuberculosis 4 (acc. NCT02607449). Future studies will identify specific aspects of the action of the drug on other multidrug resistant pathogens causing nosocomial infection outbreaks to allow better formulation of novel drugs against each individual pathogen. The comparison of the gene expression patterns of S. aureus and E. coli carried out in this work has showed several important particularities of the responses of these bacteria to FS-1. While S. aureus suffered mostly from DNA damaging, oxidative and carbonyl stresses, E. coli seemed to be more protected from the influxes of iodine into the cells possibly due to much stronger inhibition of almost all transmembrane nutrient uptake system and cytochrome:ubiquinol electron transfer complexes. However, this strategy led the FS-1 treated E. coli culture to a nutrient stringency and stronger osmotic stress. Further studies will be performed to ensure that these responses including the antibiotic resistance reversion are inherent also in clinical isolates of nosocomial infections. This work will be performed on a collection of multidrug resistant strains isolated recently from clinics in Kazakhstan, sequenced and deposited at NCBI (Korotetskiy et al., 2019b).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/nuccore/CP033505.1 and https://www.ncbi.nlm.nih.gov/ nuccore/CP033506.1.

AUTHOR CONTRIBUTIONS
OR did the data processing and visualization, funding acquisition, manuscript writing and edition. IK did the sequencing, data processing and visualization, manuscript writing and edition. MJ did the data processing and visualization, manuscript writing and edition. SS did the sequencing, data processing and visualization. AJ did the microbiological procedures, data processing and visualization, manuscript writing and edition. NS did the microbiological procedures, data processing. AI did the project management and conceptualization, funding acquisition, manuscript writing and edition. All authors contributed to the article and approved the submitted version.

FUNDING
Sequencing was funded by the grant (O.0776) from the program, "Study on the reversion of antibiotic resistance in pathogenic microorganisms" provided by the Industrial development and industrial safety committee of the Ministry of Industry and Infrastructural Development of the Republic of Kazakhstan. Genome assembly, annotation, bioinformatic analysis and student support were funded by the South African National Research Foundation (NRF) grant (105996).