Development of Resistance in Escherichia coli ATCC25922 under Exposure of Sub-Inhibitory Concentration of Olaquindox

Quinoxaline1,4-di-N-oxides (QdNOs) are a class of important antibacterial drugs of veterinary use, of which the drug resistance mechanism has not yet been clearly explained. This study investigated the molecular mechanism of development of resistance in Escherichia coli (E. coli) under the pressure of sub-inhibitory concentration (sub-MIC) of olaquindox (OLA), a representative QdNOs drug. In vitro challenge of E. coli with 1/100× MIC to 1/2× MIC of OLA showed that the bacteria needed a longer time to develop resistance and could only achieve low to moderate levels of resistance as well as form weak biofilms. The transcriptomic and genomic profiles of the resistant E. coli induced by sub-MIC of OLA demonstrated that genes involved in tricarboxylic acid cycle, oxidation-reduction process, biofilm formation, and efflux pumps were up-regulated, while genes involved in DNA repair and outer membrane porin were down-regulated. Mutation rates were significantly increased in the sub-MIC OLA-treated bacteria and the mutated genes were mainly involved in the oxidation-reduction process, DNA repair, and replication. The SNPs were found in degQ, ks71A, vgrG, bigA, cusA, and DR76-4702 genes, which were covered in both transcriptomic and genomic profiles. This study provides new insights into the resistance mechanism of QdNOs and increases the current data pertaining to the development of bacterial resistance under the stress of antibacterials at sub-MIC concentrations.


Introduction
In recent years, resistant bacteria infection is one of the most pressing concerns of global societies, which has substantially contributed to the clinic and economic burden of the healthcare system [1]. During the use of antibacterials, microorganisms are frequently exposed to drug concentrations below the minimal inhibitory concentration (MIC), which is called sub-inhibitory concentration (sub-MIC) [2]. The origins of sub-MIC concentration of the antibacterial environment are versatile. For example, the concentration of natural antibacterials produced by bacteria and fungi in the environment is often lower than the MIC of pathogens [3]. Furthermore, carcasses and excreta of animals that have used antibacterial drugs generate antibacterial gradients in the body and in the surrounding environment [4]. The poor therapeutic schedule of antibacterials treatment and preventive use of antibacterials, as well as resistance will provide new insights for assessing the risk of antibacterial resistance caused by the low concentration of antibacterial agents.

Phenotypic
Changes of E. coli ATCC25922 Induced by Sub-MIC Of OLA In Vitro 2.1.1. Resistance Development of E. coli ATCC25922 Exposed to Sub-MIC of OLA Under aerobic conditions, the MIC value of OLA against E. coli ATCC25922 was 8 µg/mL. E. coli ATCC25922 was induced by sub-MIC of OLA in four dosage groups, 1/2× MIC(4 µg/mL), 1/4× MIC(2 µg/mL), 1/10× MIC(0.8 µg/mL), and 1/100× MIC(0.08 µg/mL). As shown in Figure 1, with the number of passages increased, the resistance rates gradually increased. From the 200th generation, 1/2× MIC and 1/4× MIC OLA-induced bacteria appeared colonies on plates containing 2, 4, and 8 times MIC of OLA, and the highest resistance rate was above 3% in the 600th generation. Bacteria induced by 1/10× MIC and 1/100× MIC OLA survived against 2× MIC and 4× MIC of OLA from the 300th and 500th generations, respectively, and the highest resistance rate was below 1% in the 600th generation. There were almost no drug-resistant colonies on the 16× MIC OLA plate for all the four groups, and 1/100× MIC OLA-induced group even had no resistant-colonies to 8× MIC OLA. Differently, 1/2× MIC and 1/4× MIC Enrofloxacin (ENR)-induced E. coli appeared colonies on 2× MIC and 4× MIC ENR plates from 100 generations, and the resistance rates for 2× MIC ENR exceeded 20% by the 600th generation ( Figure S1). Bacteria induced by 1/10× MIC and 1/100× MIC ENR survived against 2× MIC and 4× MIC of ENR from the 200th and 400th generations, respectively, and the resistance rate of the 1/10× MIC ENR-induced group to 2× MIC was more than 1% at the 500th generation, and the resistance rate of the 1/100× MIC ENR-induced group against 2× MIC reached 0.1% at the 600th generation. The resistance rates of the four groups in plates containing 8× MIC ENR at the 600th passage were 1.9%, 0.9%, 0.02%, and 0.003%, respectively. Except 1/100× MIC ENR group, the other three groups could induce resistant colonies against 16× MIC ENR at the 600th passage eventually.

Effects of Sub-MIC of OLA on the Morphology and Biofilm Formation of E. coli ATCC25922
Scanning electron microscopy of 1/2× MIC OLA or ENR-induced bacteria showed that the sub-MIC of antibacterials had little effect on the growth or morphology of E. coli ATCC25922 ( Figure S2). As shown in Figure 2, the biofilm formation ability of E. coli increased slightly with the increase of number of passages. Compared with control group, most of them belong to weak biofilm formation strains (OD value was 1 to 2 times that of the control group), and strong biofilm formation strains (OD value was 2 times larger than that of the control group) are 500th and 600th bacteria induced by 1/10× MIC and 600th bacteria induced by 1/4× MIC OLA.  (Table S1). The reads were mapped on the reference genome of E. coli ATCC25922 (NCBI accession No. GCF_000743255.1) (Supplementary Table S2). The closer the value of the correlation coefficient between samples is to 1, the more similar the expression patterns indicated. The correlation coefficients for gene expression pattern between biological triplicates of untreated, 1/10× MIC and 1/2× MIC OLA-treated bacteria were all ≥0.90, indicating that the correlation between triplicate samples in the same group is good. The square of the Pearson correlation coefficient (R 2 ) varied from 0.80 to 0.97 in 1/2× MIC OLA-treated group and from 0.85 to 0.94 in 1/10× MIC OLA-treated samples, indicating good operational stability and reliability of sequencing library (Supplementary Figure S3). The functional annotation was performed based on the KEGG, NOG, GO, and NR databases. The annotation of unigenes in the Swiss-Prot, KO, eggNOG, GO, and NR databases accounted for 86.73%, 63.22%, 13.09%, 49.43%, and 99.76% of the total unigenes, respectively (Supplementary Table S3), indicating that the quality of transcriptome sequencing was good and the function of most genes was annotated. The RNAseq data were deposited in the NCBI Sequence Read Archive (SRA) database under the accession number of PRJNA525186.
DEGs between OLA-treated and untreated groups were analyzed to better understand the mechanism of bacterial resistance development against OLA (Supplementary Tables S4 and S5). Compared to the untreated group (group B), 156 genes were up-regulated and 130 were down-regulated in 1/2× MIC OLA treated group (group C). A comparison of gene expression between group B and 1/10× MIC treated group (group D) showed that 255 genes were differentially expressed, of which 101 genes were up-regulated and 154 genes were down-regulated. A total of 269 differentially expressed genes between the group C and group D with 85 genes up-regulated and 184 genes down-regulated. An investigation of the lists of DEGs by Venn diagram found that a number of the same genes were up-regulated or down-regulated in the OLA-treated groups compared with the untreated group (Supplementary Figure S4). There were 121 co-differentially expressed genes in the OLA-treated groups, 38 genes were up-regulated and 83 genes were down-regulated. Co-differentially expressed genes which have been annotated in the GO and KEGG databases, are shown in Table 1. Among them, three of the up-regulated genes are related to transcriptional regulators, fourteen genes are involved in the redox process of bacterial metabolism, six genes encode phosphoglucose transferase system (PTS) proteins, six genes are involved in sugar and amino acid metabolism, and six genes are associated with biofilm formation, ferrous iron transport, anti-adapter protein, formate C-acetyltransferase, sensor histidine kinase, and one uncharacterized protein, respectively. Only 21 down-regulated genes were annotated in the GO and KEGG databases. Nine genes are involved in DNA binding and repair, four genes encode the membrane component proteins, two genes encode chemotactic proteins, and two genes encode phage shock proteins. Other five genes encode pyrimidine utilization transport protein, putative electron transport protein, aminopeptidase N, toxin HigB-1 and a hypothetical protein, respectively. The proteins corresponding to the 121 co-differentially expressed genes were matched with the proteins of Escherichia coli K12 MG1655, and then a cluster analysis of protein interactions was performed with 15 significant KEGG metabolic pathways obtained in the STRING database ( Figure 3). They are phosphotransferase system (PTS), metabolic pathways, alanine, aspartate and glutamate metabolism, carbon metabolism, biosynthesis of antibiotics, microbial metabolism in diverse environments, biosynthesis of secondary metabolites, glyoxylate and dicarboxylate metabolism, amino sugar and nucleotide sugar metabolism, fructose and mannose metabolism, glycolysis/gluconeogenesis, methane metabolism, tyrosine metabolism, butanoate metabolism, oxidative phosphorylation. In this network, oxidoreductases, such as Pfo, AdhE, and NuoM were shown to play a role in the important processes of pyruvate-ferredoxin/flavodoxin redox, DNA replication, and NADH-quinone redox. Because these proteins interact with more proteins, they are the key nodes of the protein interaction network.
Compared with the reference genome sequence of E. coli ATCC25922 (NCBI accession number GCF_000743255.1), RNA-seq results showed that the number of SNPs in the genome were 193, 120, and 122 from the three duplicates of 1/10× MIC OLA treated group, 128, 146, and 129 from the three duplicates of 1/2× MIC OLA treated group, and 86, 92, and 85 from the three duplicates of untreated group, respectively. RNA-seq results also showed that the number of indel sites in the genome was 30, 21, and 22 from the three duplicates of 1/10× MIC OLA treated group, 19, 24, and 18 from the three duplicates of 1/2× MIC OLA treated group, and 13, 18, and 12 from the three duplicates of untreated group, respectively. This result demonstrated that the number of mutations and insertion/deletion sites in the OLA treated groups was significantly increased. Among these SNPs, 23 SNP sites would result in a protein sequence change, 15 non-synonymous mutations, and eight single-point deletions ( Table 2). Genes with these mutations encode regulatory proteins of transcription, nucleotide pyrophosphatase family proteins, efflux system membrane proteins, etc. The SNPs were compared with the CARD resistance database, and the results showed that there were no mutations of the known drug resistance genes. Protein-protein interactions of E. coli ATCC25922 induced by sub-MICs of OLA generated following analysis of significantly different peptides input into the STRING database. The figure shows the protein-protein interactions networks that resulted from co-differentially expressed genes. " " represent known interactions from curated databases, " " represent known interactions have been experimentally determined; " " represent predicted interactions of gene neighborhood, " " represent predicted interactions of gene fusions, " " represent predicted interactions of protein homology. " " represent interactions of gene textmining, " " represent predicted interactions of gene co-expression, " " represent interactions of gene co-occurrence.

Genomic Profiles of E. coli ATCC25922 Induced by Sub-MIC Of OLA
The genome of E. coli ATCC25922 induced by 1/2× MIC OLA and resistant to 8× MIC OLA showed a total bases of 529,982,985 bp and a number of subreads of 64,908, with average subreads length of 8165 bp. Using the Canu software, the genome was assembled and showed a 5,158,951 bp length (Supplementary Figure S5), which contained 5384 predicted genes with an average G+C content of 50.49%. The protein sequences predicted by Glimmer software were extracted, and the length of the protein sequences were statistically plotted (Supplementary Figure S6). There were 4045 and 3626 predicted genes annotated in COG (Supplementary Figure S7) and GO database, respectively (Supplementary Figure S8). Also, 3181 predicted genes were involved in 29 KEGG pathways which were classified into "Metabolism" (12 out of 29), "Genetic Information Processing" (four out of 29), "Cellular Processes" (four out of 29), "Environmental Information Processing" (two out of 29), and "Organismal Systems"(seven out of 29) (Supplementary Figure S9). The genomic sequencing data were deposited in the NCBI Sequence Read Archive (SRA) database under the accession number of PRJNA525188. Using the genomic multiple alignment tool Mauve, the genome of the OLA resistant E. coli was compared with the reference genome of E. coli ATCC25922 (NCBI accession No. GCF_000743255.1). There are two gene fragments of shifting and one gene fragment of deletion ( Figure 4).  These lines indicate which regions in each genome are homologous. When a block lies above the center line the aligned region is in the forward orientation relative to the first genome sequence. Blocks below the center line indicate regions that align in the reverse complement orientation. Inside each block Mauve draws a similarity profile of the genome sequence. The height of the similarity profile corresponds to the average level of conservation in that region of the genome sequence.
Using MUMmer(v3.23), the genome of the OLA-resistant bacteria was compared with the reference genome (NCBI accession No. GCF_000743255.1), and 715 SNP sites were obtained, among which 270 SNP sites involving 41 genes existed in the exon region, causing non-synonymous mutations or codon termination. In the meantime, 82 indel sites were obtained, among which 27 indel sites existed in the exon region (Table 3). There were 113 SNP sites and 19 indel sites associated with the DNA replication fork binding protein CrfC, which plays important functions such as cell replication, recombination, and repair. Only three SNP sites and one indel locus were involved in DNA binding, transcriptional regulation, and recombination genes. There was an SNP site observed in drug/metabolite transporter gene yddG and a SNP site efflux gene acrB. There were four SNP sites and one non-frameshift deletion found in NAD + oxidoreductase gene. These results indicated that the bacterial resistance to QdNO drugs could be related to NAD + redox and efflux pump. Comparing the SNPs between the transcriptomic sequencing ( Table 2) and genomic sequencing (Table 3), six identical SNP sites, degQ, ks71A, vgrG, bigA, cusA, and DR76-4702, were found. vgrG and bigA have different mutation sites, and the other 4 genes have only one mutation site in both of the sequencing profiles. Other SNPs included amino acid metabolism, integral component of membrane, ion transport, and nucleic acid metabolism.  Using MUMmer (v3.23), the genome of the OLA-resistant bacteria was compared with the reference genome (NCBI accession No. GCF_000743255.1), and 715 SNP sites were obtained, among which 270 SNP sites involving 41 genes existed in the exon region, causing non-synonymous mutations or codon termination. In the meantime, 82 indel sites were obtained, among which 27 indel sites existed in the exon region (Table 3). There were 113 SNP sites and 19 indel sites associated with the DNA replication fork binding protein CrfC, which plays important functions such as cell replication, recombination, and repair. Only three SNP sites and one indel locus were involved in DNA binding, transcriptional regulation, and recombination genes. There was an SNP site observed in drug/metabolite transporter gene yddG and a SNP site efflux gene acrB. There were four SNP sites and one non-frameshift deletion found in NAD + oxidoreductase gene. These results indicated that the bacterial resistance to QdNO drugs could be related to NAD + redox and efflux pump. Comparing the SNPs between the transcriptomic sequencing ( Table 2) and genomic sequencing (Table 3), six identical SNP sites, degQ, ks71A, vgrG, bigA, cusA, and DR76 -4702, were found. vgrG and bigA have different mutation sites, and the other 4 genes have only one mutation site in both of the sequencing profiles. Other SNPs included amino acid metabolism, integral component of membrane, ion transport, and nucleic acid metabolism. bigA putative surface-exposed virulence protein 10 NS

Phenotype Characteristics of Sub-MIC Antimicrobial-Induced Bacterial Resistance
The previous study has demonstrated that two-fold ascended concentrations of OLA(from 1/2× MIC to 4× MIC) can contribute to the selection and maintenance of OLA-resistant bacteria against 16× MIC OLA in 20 days [13]. However, in this study, under sub-MIC concentrations of OLA (from 1/100× MIC to 1/2× MIC), bacteria took a longer time (60 days) to enrich resistants, with maximum resistant level of only 8× MIC (Figure 1). The step-wise selection by an antibacterial at concentration above MIC will enrich the existing resistant sub-population that acquire resistance by mutation or horizontal gene transfer, while sub-MIC selection can also drive the emergence of de novo resistant mutants [21], by increasing the mutation rate or stimulating the horizontal gene transfer in susceptive bacteria which do not die in sub-MIC concentrations [22]. Moreover, for the resistant mutants that already exist in the population, the sub-MIC concentration would enrich those with less fitness cost than the susceptive counterpart [10]. So, the exposure level of drug in clinical settings will influence the speed and the level of bacterial resistance development.
It has been demonstrated that stress responses gene soxS and type I fimbriae gene (fimG), metabolism gene (metK) were differentially expressed in all E. coli biofilm studies [23]. bssR, encoding a biofilm regulator, was discovered for the induction of biofilms in E. coli [24]. In this study, soxS and bssR as well as many carbon catabolism related genes were up-regulated in the OLA induced resistant bacteria (Table 1). Especially, the bssR gene was up-regulated by 2 times in the 1/2× MIC treated group and five times in the 1/10× MIC group, which was consistent with the amount of biofilm formation (Figure 2).
A phenomenon termed collateral sensitivity, whereby resistance to one antimicrobial simultaneously increases the susceptibility to another, was reported [25]. Collateral susceptibilities of mutants resistant to relevant antimicrobials against 16 antibiotics have been demonstrated, and it was shown that resistance mechanisms, in particular efflux-related mutations, as well as the relative fitness of resistant strains, are principal contributors to collateral responses [26]. Csaba et al also reported that various mechanisms of resistance, including target mutations and those mutations affecting drug uptake and efflux, are prone to induce collateral sensitivity [27]. Our data showed that resistance to OLA resulted in collateral responses to mequindox, ampicillin, compound sulfamethoxazole, florfenicol. However, the susceptibility to colistin of the OLA resistants decreased. The mechanism of action of mequindox is targeted to DNA and compound sulfamethoxazole affect the synthesis of folate. Florfenicol acts on the 50S ribosome of bacteria and amoxicillin inhibits the formation of bacterial cell walls. These results suggest that collateral susceptibility may be associated with antimicrobial resistance mechanism. The above results show that evolution of resistance towards a given antimicrobial agent frequently increases resistance to several other drugs. but it has remained unclear how frequently evolution of resistance increases sensitivity to other drugs. Cross-resistance profiles may be strongly correlated with drug resistance trajectories and the antibacterial mechanism of antibiotics. In addition, fitness costs may slow the evolution of the collateral responses.

Molecular Mechanism of Sub-MIC of QdNO-Induced Resistance
The PTS system located in the inner membrane is a multiprotein phosphorylation system that couples the transport of carbohydrates across the cytoplasmic membrane with their simultaneous phosphorylation. The global repressor Mlc regulates the expression of phosphoenolpyruvate: sugar phosphotransferase system (PTS) operon which is related to the uptake and utilization of sugars. It has been reported that mlc inhibits manXYZ, malT, ptsG, and pts operon, encoding enzyme II of the mannose PTS, the activator of maltose operon [28], enzyme IICB of the glucose PTS, and other PTS proteins, respectively. Our results showed that mutation appears among the mlc sequence, so the mlc may lose the ability of inhibiting manXYZ, malT, ptsG, and pts operon. MalT is a transcriptional regulator of the maltose regulon, which is involved in regulating the transport of maltose [29]. Under sub-MIC of OLA, it was shown that the expression of PTS system genes including mlc, manX, manY, manZ, ptsG, ptsH, ptsI, and non-PTS system gene malT were up-regulated. It can be speculated that the expression of the mannose phosphate transfer system is up-regulated, and the maltose transport system in the non-PTS is also more active. This indicates an increase in the uptake of carbon sources in resistant bacteria. Compared to the wild-type E. coli ATCC25922, some genes involved in glycolysis were up-regulated in the resistant bacteria induced by 1/2× MIC OLA, including fbaA (fructose-bisphosphate aldolase), fbaB (fructose-bisphosphate aldolase), and gap(glyceraldehyde 3-phosphate dehydrogenase). In 1/10× MIC OLA induced bacteria, there are also some genes involved in glycolysis up-regulated, including pgk (phosphoglycerate kinase), fbaA, pfkB (6-phosphofructokinase), ppsA (pyruvate, water dikinase), and gldA (glycerol dehydrogenase). Genes involved in pyruvate metabolic pathways, such as pfo (nifJ), were up-regulated. So, the pyruvate metabolism was up-regulated and produces more acetyl-CoA for the energy system. Moreover, many other metabolic pathways, such as glycine, histidine, and aspartate metabolism, were involved in the response to OLA as well. Genes involved in amino acid metabolism were also up-regulated, such as gcvP, zraS, and asnB. The metabolism of the major carbon and nitrogen substrates are afforded by multiple specific uptake systems and pathways. The genes involved in these pathways are controlled by both global and specific transcriptional regulators that may directly or indirectly sense key metabolites that are generated in the central carbon and nitrogen metabolic pathways [30]. Our study indicates that the metabolism of sugars and amino acids is enhanced in E. coli resistant to OLA, probably to remove oxidative damage and accelerate the synthesis of new substances to maintain the normal physiological function of the cells.
The adhE (aldehyde dehydrogenase) and icdA (isocitrate dehydrogenase) were up-regulated in this study, so the yield of NADPH increased. In the meantime, the over-expression of adhE increased the yield of acetaldehyde, thus reducing the amount of acetyl-CoA in the TCA cycle, resulting in the decrease of NADH/NAD ratio [31], playing a protective role against oxidative stress [32]. NifJ, a pyruvate ferredoxin (flavodoxin) oxidoreductase, catalyzes the oxidation of pyruvate to acetyl-coenzyme A and carbon dioxide in bacteria. It was revealed that NifJ can balance NADH-mediated process by regulating the NADH/NAD ratio [33], so the overexpressed of the nifJ gene also cause the reduction of NADH/NAD ratio. wrbA is a tryptophan repressor-binding protein that binds one molecule of flavin mononucleotide, which is involved in oxidative defense [34]. It appears to be a big electron acceptor pool to accommodate the electron of the electron donor, NADH or benzoquinone [35]. The CydAB quinol oxidases catalyze the reduction of oxygen by ubiquinol. cydB gene encodes a subunit of CydAB quinol oxidases, and the accumulation of NADH was prevented when the cydB mutated [36]. CydAB also play an important role in protection against oxidative stress by reducing H 2 O 2 content [37]. nuoM and nuoI encode NADH:ubiquinone oxidoreductase, whose hydrophilic domain catalyzes transferring two electrons from NADH to quinone [38]. In the resistant bacteria induced by the sub-MIC of OLA in our study, adhE, nifJ, wrbA, cydB, nuoM, and nuoI were significantly up-regulated (Table 1). It was reported that QdNOs was metabolized by a reductase to produce free radicals in bacteria and this process also involved electron transfer [39]. Therefore, reducing the accumulation of NADH by decreasing the NADH/NAD ratio and increasing the NADPH/NADP ratio might be one of the mechanisms of QdNO resistance.
MarA, SoxS and Rob are three homologous transcriptional regulators of E. coli. These homologous regulators are part of multiple regulatory mechanisms necessary for the adaptive response, including changing pH, the presence of antibiotics, oxidative stressors, and organic solvents, all of which threaten survival [40]. The soxRS regulon is involved in oxidative stress response and also responsible for the resistance of E. coli to antibacterials [41,42]. Only the upregulation of SoxS was found in our study, but not MarA and Rob. The possible reason for this is that cross talk between the mar and sox systems is limited under physiological conditions. Previous study showed that the expression of SoxS is not expected to activate marRAB transcription in the absence of salicylate, when MarR still represses the promoter [42]. A previous study also showed that SoxR is regulated by NADPH-dependent reductases, an increased NADPH demand of the cell counteracts SoxR reduction and increases soxS expression [43]. Several genes encoding NADPH-dependent reductases were up-regulated, such as nuoG, nuoH, and nuoI. So, the up-regulation of these genes may cause the up-regulation of soxS. Previous study showed that the alkyl hydroperoxide reductase (AhpC) has the function of scavenging ROS and peroxides in variety of aerobic and anaerobic bacteria [44,45]. Previous research has suggested that the knockout of grxB gene, which encodes glutaredoxins ubiquitous protein, cause lower levels of oxidoreductase activity, resulting in the increased sensitivity to oxidative stress in E. coli [46]. Overproduction of the electron transfer subunits of formate dehydrogenases, encoded by the fdoH, fdoI, and fdnG, can enhance oxidative stress tolerance [47,48]. In our study, soxS, ahpC, grxB, fdoH, fdoI, and fdnG gene were up-regulated. Our previous study shows that QdNO-induced DNA damage is related to the yield of free radicals, and free radical scavengers can alleviate the bactericidal activity of cyadox [49]. Therefore, free radical scavenging and oxidative stress tolerance may be one of the mechanisms of QdNOs resistance.
umuC encodes polymerase V, a low-fidelity DNA polymerase, which has the ability to trigger a repair pathway during the SOS response to DNA damage [50]. The RusA protein is a DNA structure-specific endonuclease that resolves holliday junction intermediates formed during DNA repair [51]. RutG, a broad-specificity pyrimidine permease, is predicted to be a pyrimidine transporter [52]. Previous study demonstrated that the genes involved in DNA metabolic process and SOS response were downregulated, including dinB, umuC, recA, etc. [53]. In this study, DNA damage repair genes, umuC, rusA, and rutG were down-regulated. The resistant bacteria may resist DNA damage, but the susceptive bacteria can be damaged by antibiotics. Therefore, the expression of genes related to DNA repair in drug-resistant bacteria is down-regulated.
A total of 50 non-synonymous SNPs and indels were observed between the genomes of 8× MIC OLA resistant bacteria and the reference strain (NCBI accession No. GCF_000743255.1). IraD is a factor aiding the survival of proliferating E. coli cells against multiple forms of DNA damage, including oxidative stress [54]. CdiA secretion protein provides immunity and inhibits the bacterial growth by a significant down-regulation of metabolic parameters, such as aerobic respiration, steady-state ATP levels, etc. [55]. A previous study suggests that the Type VI Secretion System core component VgrG contributes to the antimicrobial resistance in A. baumannii ATCC19606 [56]. The gene ydeE encodes a multidrug-efflux transporter of E. coli which confers resistance to dipeptides [57]. Previous studies reported that the cytotoxicity of bactericidal antibacterials predominantly results from lethal double-strand DNA breaks caused by incomplete repair of 8-oxo-deoxyguanosine (8-oxo-dGTPase), which is encoded by mutT gene to prevent incorporation of base analogs into DNA [58]. iraD, cdiA, vgrG, ydeE, and mutT genes involved in DNA damage, bacteria growth and antimicrobial resistance were mutated in our study, which may be related to the resistance of QdNOs. Previous study demonstrates that CrfC protein, distributed in the nucleoid poles, helps sustain the colocalization of nascent DNA regions of sister replisomes and promote chromosome equipartitioning [59]. So, the CrfC might be involved in reactions necessary for chromosome stability and protect DNA damage.
Comparing the SNPs in the transcriptomic data ( Table 2) with those in the genomic data (Table 3), six identical SNP genes were found. They are degQ, ks71A, vgrG, bigA, cusA, and DR76-4702. DegQ, a HtrA protease of the E. coli, is involved in the biosynthesis of outer membrane proteins [60,61]. ks71A encodes the fimbrial subunits of E. coli Ks71A, which is a hypothetical protein [62]. vgrG is a valine-glycine repeat G gene, previous study showed that the deletion of the vgrG gene, increased antimicrobial resistance to ampicillin/sulbactam, but reduced resistance to chloramphenicol [56]. The vgrG gene was mutated in the genome profiles. So, the vgrG gene might be related to OLA resistance. bigA encodes a putative surface-exposed virulence proteinin E. coli, and it also has been reported that the BigA protein acted as an adhesin protein of Brucella that mediates adhesion to epithelial cells [63]. CusA is an inner membrane proton-substrate carrier, responsible for transferring Cu(I) and Ag(I) ions [64]. The DR76-4702 gene has not been named, encoding a putative integral component of membrane. The outer membrane (OM) of Gram-negative bacteria performs the crucial role of providing an extra layer of protection to the organism, so the mutation of these genes provided a formidable barrier for bacteria to sustaining life.
A previous study has showed that the mechanism of an oxidative damage cellular death pathway involving the tricarboxylic acid cycle, a transient depletion of NADH, destabilization of iron-sulfur clusters, and stimulation of the Fenton reaction [65]. It was demonstrated that QdNOs genes and proteins involved in the bacterial metabolism, cellular structure maintenance, resistance and virulence were found to be changed, accompanying the SOS response and oxidative stress induced [49]. A previous study revealed that E. coli adapt the stress generated by QdNOs or develop specific QdNOs-resistance by the activation of antioxidative agent biosynthesis, protein biosynthesis, glycolysis, and oxidative phosphorylation [13]. Based on the mechanism of cell death induced by QdNOs and oxidative stress-mediated mode of action proposed for QdNOs in the above literature, it is supposed that the QdNOs resistance mechanism involved in outer membrane porin, efflux pump, the tricarboxylic acid (TCA) cycle, the oxidation of NADH, damage DNA repair and DNA replication and mutation of E.coli ( Figure 5). The study only elaborated the resistance mechanism of E. coli to OLA at sub-inhibitory concentrations of OLA, and could not fully cover all the resistance mechanisms of OLA under concentrations above MIC. This study only studied the OLA resistance mechanism of QdNOs, which may not represent the resistance mechanism of all QdNOs. Therefore, it is necessary to study the resistance mechanism of other QdNOs.

Drugs, Bacteria and Growth Conditions
OLA (purity of 99.8%) was purchased from National Institutes for Food and Drug Control (Beijing, China). Enrofloxacin (ENR) was bought from China Institute of Pharmaceutical and Biological Products Inspection (Beijing, China). The standard solutions of drugs were prepared Figure 5. The hyperthetical resistance mechanism involved in the development of QdNOs resistance. The decrease in the expression of the outer membrane porin gene ompD leads to the loss of OmpD protein, which reduces the intake of olaquindox. The up-regulation of multidrug resistance efflux pump gene acrB increases olaquindox efflux. The PTS systerm initiate a metabolic feedback dependent on tricarboxylic acid (TCA) cycle, stimulate the oxidation of NADH through the electron transport chain, promote superoxide formation, damage the structure of Fe-S clusters, lead to the formation of hydroxyl radicals and finally damage DNA; The Changes in the expression of genes umuC, rusA, and rutG trigger DNA repair; Changes in the expression of genes such as iraD, cdtA, vgrG, and ydeE cause gene mutations and affect DNA replication.

Drugs, Bacteria and Growth Conditions
OLA (purity of 99.8%) was purchased from National Institutes for Food and Drug Control (Beijing, China). Enrofloxacin (ENR) was bought from China Institute of Pharmaceutical and Biological Products Inspection (Beijing, China). The standard solutions of drugs were prepared according to the manufacturer's instruction at concentration of 1280 mg/L. E. coli ATCC25922 was purchased from China Center for Type Culture Collection (CCTCC). Mueller-Hinton broth, Mueller-Hinton agar (Mueller-Hinton broth supplemented with 1.5% agar) and Luria-Bertani (LB) agar was bought from HOPEBIO (Tsingtao, China). All other chemical reagents were of analytical grade.

Antimicrobial Susceptibility Testing
The MICs of antibacterials for wild-type and OLA-resistant E. coli ATCC25922 were determined using the broth micro-dilution method, according to the guidelines of the Clinical and Laboratory Standards Institute (CLSI) [66]. Specifically, 0.1 mL of two-fold serial dilutions of drugs were prepared and inoculated with 0.1 mL of bacterial overnight culture diluted into 5×10 5 CFU/mL. The tubes were incubated at 37 • C with shaking for 16-18 hours, the OLA cultures protected from light to avoid degradation of the antibiotic. The MIC was defined as the lowest concentration of an antimicrobial able to completely inhibit the microbial growth, as described in the CLSI protocols.

In Vitro Selection of Resistant Bacteria under Sub-Inhibitory Concentrations of OLA
To select for de novo generated resistant mutants at sub-MICs of OLA, a single colony of E. coli ATCC25922 was transferred into 5 mL of Mueller-Hinton broth and incubated at 37 • C overnight. Fresh cultures of bacteria were serially passaged by 1000-fold dilution in 5 mL batch cultures every 24 hours for 600 generations (60 passages, 10 generations per passage) in Mueller-Hinton medium containing OLA at concentrations of 1/2× MIC (4 µg/mL), 1/4× MIC (2 µg/mL), 1/10× MIC (0.8 µg/mL), and 1/100× MIC (0.08 µg/mL) or ENR at concentrations of 1/2× MIC (0.0075 µg/mL), 1/4× MIC (0.00375 µg/mL), 1/10× MIC (0.0015 µg/mL), and 1/100× MIC (0.00015 µg/mL) respectively, and the untreated E. coli culture was used as control. The liquid cultures were grown at 37 • C under aerobic conditions without shaking [10]. For every 100 generations of the bacteria culture, the number of colonies was counted, and the percentage of resistant cells in each culture was monitored by plating approximately 10 5 cells onto LB agar containing different concentrations of OLA or ENR as control (2× MIC, 4× MIC, 8× MIC, and 16× MIC). A subset of these cells was restreaked on the same antibiotic concentration to confirm that they were resistants.

Bacterial Morphology Observation and Biofilm Detection
To investigate the morphological changes of bacteria under sub-MIC of antimicrobial drugs, E. coli ATCC25922 induced by 1/2× MIC OLA and 1/2× MIC ENR as well as untreated group were observed by scanning electron microscopy. Immobilization, embedding, sectioning, staining and observation of bacterial samples were carried out as previously described [67]. Bacteria (about 10 6 CFU/mL) were harvested by centrifugation, washed twice with 0.1 M PBS (pH 7.4), and fixed in 0.2% glutaraldehyde in 0.1 M cacodylate buffer at 4 • C for 4 h, and then fixed in the perfluorocarbon containing 1% osmium tetroxide for 1 h. After three rinses in pure perfluorocarbon, the samples were dehydrated, embedded, sectioned, stained, and imaged with a TEM (H-7650, Japan) as previously described [68].
The method of biofilm formation detection of bacteria at various stages induced by 1/2× MIC, 1/4× MIC, 1/10× MIC, and 1/100× MIC of OLA was as follows [69]. After overnight incubation to logarithmic growth phase, the cultures were diluted to 0.5 McFarland standards in broth medium, respectively. A volume of 200 µL bacterial cultures was transferred to each well of a 96-well plate, and bacterial biofilms were formed followed by incubation at 37 • C for 24 h with no movement. After discarding the medium and rinsing the wells of the 96-well plate 3 times with sterile PBS solution, the biofilms were placed in a 0.1% (wt/vol) crystal violet solution for 15 min, rinsed again, and dried for several hours. To solubilize the adsorbed crystal violet, wells with stained biofilms were incubated in 30% acetic acid for 15 min. The absorbance was read at 590 nm on a microplate reader (Thermo Fisher Scientific, Shanghai, China). were pooled for RNA-Seq library construction. The mRNA was enriched from qualified total RNA using Oligo(dT) magnetic beads and fragmented chemically into small pieces at high temperature. Then, the first cDNA strand was synthesized by using random hexamers for reverse transcription with cleaved RNA fragments serving as templates, followed by second-strand cDNA synthesis using DNA polymerase I and RNase H. The double-strand cDNA then went through purification, end repair process, the addition of a single "A" base and ligation of the adapters, and was finally enriched by PCR implication to be ready for sequencing via Illumina HiSeq™ 2000 (Personalbio technology Co. Ltd., Shanghai, China).

Transcriptome Assembly and Gene Functions Annotation
A total of 39.0 GB raw data was obtained from nine samples via Illumina system. Clean reads were obtained after removing adapter sequences, filtering low-quality reads, and removing reads with an N ratio of greater than 5% from the raw reads. Each sample had more than 89.67% of clean data, with Q30 values greater than 90.03%. After adaptor sequences, ambiguous reads, and low-quality reads were removed, the sequence data were assembled with Trinity software [70]. RNA sequence reads were independently aligned with the genome sequence of E. coli ATCC 25922 (NCBI accession number GCF_000743255.1). Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); and GO (Gene Ontology). Principal component analysis was performed using the R statistical computing package [71]. Gene ontology analysis was performed by Blast2GO through searching the above databases [72]. The GO terms assign biological process, molecular function, and cellular component to the query sequences. KEGG analysis was performed using the online KEGG Automatic Annotation Server (KAAS; http://www.genome.jp/kegg/kaas) to obtain an overview of gene pathways networks.

Detection of DEGs
Gene expression was analyzed by DESeq (version 1.18.0). The screening conditions for DEGs included expression fold difference of log2|fold change| > 1 and significance of p-value <0.05.

Identification of Antimicrobial Resistant Genes
The set of resistance-associated genes was downloaded from the Antibiotic Resistance Genes Database (http://ardb.cbcb.umd.edu/). DEGs were annotated using the tool retrieved from the same web site, which was performed based on a homology search by BLASTp. The cutoff parameters used in this study were E-value <1e −5 and nucleotide sequence identity >50 %.

DNA Preparation, DNA Library Construction and Sequencing
Genomic DNA of resistant E. coli ATCC 25922 induced by 1/2× MIC of OLA and exhibiting a MIC of 64 µg/mL(equal to 8× MIC for wild-type susceptive strain) was extracted using a bacterial genomic DNA purification kit (Tiangen Biotech Co., Ltd., Beijing, China) according to the manufacturer's instructions. DNA was quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The isolated DNA was stored at -80 • C until use. Whole-genome sequencing was carried out using the PacBio Sequel platform (Frasergen Co. Ltd., Wuhan, China). A quantity of 5 µg input genomic DNA was used for 10 kb fragment library preparation. DNA was sheared with g-TUBE1 microcentrifuge tubes (Covaris, Woburn, MA) and SMRTbell libraries prepared with DNA Template Kit 2.0 (Pacific Biosciences, Menlo Park, CA). Sequencing primers were annealed to the SMRTbell template and samples were sequenced using C3 chemistry, Polymerase (version 5), SMRTAnalysis software (version 2.3.0), and a single SMRT cell (Pacific Biosciences, Menlo Park, CA). Raw read quality and average reference consensus were determined using SMRT Analysis software (version 2.3.0).

Genome Assembly and Annotation
Filtered reads were obtained from PacBio raw data using SMRT Analysis and were assembled by applying the Canu protocol [73]. The assembled genome sequence was annotated by the Institute for Genomic Research database, using the programs Glimmer (version 3.02) for the identification of protein-coding genes [74], tRNAscan-SE for the identification of tRNA genes [75], and RNAmmer for the identification of rRNA genes [76]. The non-coding RNA were identified using the Rfam (version 12.0) and Infernal [77]. Short tandem repeat and simple sequence repeats were identified by trf407b.linux [78], and clustered regularly interspaced short palindromic repeats (CRISPRs) were identified using CRISPRs in Environmental Datasets [79].

Statistical Analysis
All statistical analyses were performed using GraphPadPrism software (GraphPad Software Inc., LaJolla, CA, USA) with all data represented as the mean ± standard deviation (SD) from at least three independent experiments. Data analyses were performed using the t-test or one-way analysis of variance (ANOVA), and a p-value less than 0.05 was considered statistically significant.

Conclusions
Our data show the evolution of resistance of E. coli de novo selected by sub-MIC concentrations of OLA, and provide a reference predicting collateral sensitivity informed therapy. QdNOs are redox-activated, hypoxia-selective DNA-cleaving agents [80]. The oxidative DNA-damaging effects of QdNOs will stimulate the SOS response, oxidative stress, and other protection strategies in bacteria. Combined with the previous study of the QdNOs resistance mechanisms [20], we further elaborated the molecular mechanism of QdNOs resistance under sub-MIC. The intake of QdNOs was first reduced by adjusting the expression of membrane porin and efflux pumps and mutations in efflux pump genes, mainly including ompD and acrB. When the QdNOs enters the bacterial cell, it stimulates the bacteria to feedback by regulating the electron transport chains, resulting in the production and consumption of NADH to reduce the production of free radicals. The main genes involved in this process are soxS, malT, nouM, nouN, and icdA. When an oxidative stress reaction occurs, the free radicals were cleared by regulating the expression of certain redox genes and genes involved in DNA repair, involving genes of soxS, umuC, iraD, ahpC, wrbA, and rusA. The bacteria cell division was halted when exposed to CYA [25], so the mutations in the crfC gene, which is involved in bacterial replication, may be another key gene of QdNOs resistance. The deeper molecular mechanism involved in the origin of resistance development selected by QdNOs should be explored. Further studies are required to investigate the signaling pathway of the origin of the QdNOs resistance. Targeting the key molecular actors involved in the gene expression and inhibiting the mutagenic pathways of the emergence of resistance could be a key strategy to slow down the emergence of antibacterial resistance.
Supplementary Materials: The following are available online at http://www.mdpi.com/2079-6382/9/11/791/s1, Figure S1: Resistance rates of E. coli ATCC25922 exposed to ENR at sub-MIC concentrations of 1/2× MIC (A), 1/4× MIC (B), 1/10× MIC (C) and 1/100× MIC (D), Figure S2. Scanning electron micrograph of E. coli ATCC25922 induced by none (A), 1/2× MIC ENR (B) and 1/2× MIC OLA(C)., Figure S3. Correlation tests for 9 samples (triplicates in each group) of the untreated (group B), 1/2× MIC (group C) and 1/10× MIC OLA (group D) -treated E. coli ATCC25922. The abscissa and ordinate in the figure are sample numbers. The closer the block value is to 1, the higher the similarity is., Figure S4. Venn diagrams of DEGs. The numbers in each circle represents the total number of DEGs in the comparison combination, and the overlapping part of the circles represents the DEGs shared between the comparison groups. B, untreated group; C, 1/2× MIC OLA-treated group; D, 1/10× MIC OLA-treated group., Figure S5. The genome of 1/2× MIC OLA induced E. coli ATCC25922 resistant to 8× MIC OLA. The circles from outermost to innermost indicated the scale, GC content (red: > average value, blue: < average value), GC skew (purple: >0, orange: <0), noncoding RNA (tRNA in black, rRNA in red), minus-stranded CDS maps, plus-stranded CDS maps, minus-stranded base-modified map (full red circle), plus-stranded base-modified map (full blue circle), gene map of the restriction modification system, respectively., Figure S6. Protein length distribution of the genome of OLA resistant E. coli ATCC25922., Figure S7. COG function classification of OLA resistant E. coli ATCC25922. The abscissa is the group of the COG, and the ordinate is the number of genes annotated to the group., Figure S8. GO function annotation of OLA resistant E. coli ATCC25922. The ordinate is the three major classes of the GO term including biological processes, cellular components, and molecular functions, the abscissa is the number of genes annotated to the term (including the term of the term)., Figure S9. KEGG functions annotation of OLA resistant E. coli ATCC25922. The ordinate is the name of the KEGG metabolic pathway, and the abscissa is the number of genes annotated to the pathway. The genes are divided into five branches according to the involved KEGG metabolic pathway including Cellular Processes, Environmental Information Processing, Genetic Information Processing, Metabolism, and Organic Systems., Table S1: Table S1. Data Filtering Statistics, Table S2. RNASeq Map Statistics., Table S3. Reference genome annotation information statistics., Table S4. Number of differentially expressed genes between groups B and C., Table S5. Differentially expressed genes between groups B and D.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.