Comparative transcriptomic analysis reveals the molecular mechanisms related to oxytetracycline- resistance in strains of Aeromonas hydrophila Aquaculture

in the 4-fold minimal inhibitory concentration (MIC) and 8-fold MIC resistant strains, respectively. Furthermore, the bioin- formatics analysis revealed that the sulfur metabolism-related genes were down-regulated. The genes responsible for mannitol metabolism and the efflux pump system were up-regulated in resistant strains, compared to the susceptible strain. Therefore, it suggests that these three pathways may be involved in the OXY resistance evolution in A. hydrophila . The outcome of the transcriptomic data was further validated through quantitative reverse transcription-PCR (qRT-PCR) and Western blot analysis. Overall, the obtained data provides a deeper insight into the intrinsic molecular mechanism of OXY resistance evolution in A. hydrophila .


Introduction
Antibiotic resistance among bacterial pathogens has been on the rise since the discovery of antibiotics, posing significant challenges for infectious disease treatment (Hasannejad-Bibalan et al., 2019). During their prolonged evolution, bacteria have developed resistance against several antibiotics using a range of tactics such as enhancing the efflux drug system, reducing the outer or inner membrane permeability, altering the structure of antibiotics and producing antibiotic hydrolytic enzymes (D'Costa et al., 2012;Kumar and Schweizer, 2005). However, the evolution of drug resistance in bacterial pathogens has been evident even before the discovery of antibiotics. Furthermore, the evolution of drug resistance is continuously established when the bacteria are confronted with harmful toxins and low doses of antibiotics. For example, genomic sequence analysis revealed that a staphylococci strain isolated from permafrost carried several antibiotic resistance genes (Kashuba et al., 2017), which indicates that some kinds of antibiotic resistance are naturally occurring due to spontaneous resistance evolution.
In recent years, research has focused on the antibiotic resistance mechanisms in various bacterial pathogens using whole-genome sequencing and transcriptomic approaches. For instance, Zhu et al. (2017) investigated the mechanism of enrofloxacin resistance in Aeromonas hydrophila using transcriptomic methods. The authors identified transcriptomic differences between enrofloxacin susceptible and resistant strains of A. hydrophila, revealing that enrofloxacin affects multiple biological functions in the resistant strain. Moreover, proteomic analysis has also been used to understand the antibiotic resistance mechanisms in several bacterial pathogens by identifying the antibiotic resistance related proteins and pathways. In our recent studies, we have identified several differentially expressed proteins (DEPs) in resistant strains of A. hydrophila under oxytetracycline (OXY) treatment by quantitative proteomics technologies Yao et al., 2018). Further to this, the outcome of these proteomic studies have revealed that various biological processes such as fatty acid biosynthesis, chemotaxis and central metabolic pathways were involved in antibiotic resistance in bacterial pathogens, other than the more classical antibiotics resistance mechanisms previously reported Li et al., 2016Li et al., , 2018bLi et al., 2018aLi et al., , 2019Su et al., 2018).
A. hydrophila is a Gram-negative bacterium belonging to the Aeromonadaceae family. It is an important fish bacterial pathogen responsible for global economic losses in the aquaculture sector. A. hydrophila is known to cause deleterious septic outbursts with high mortality rates in various fish species and poses a serious hazard to public health (Baldissera et al., 2019). In fish farms, A. hydrophila affects several types of fish species such as the Ctenopharyngodon idella (Grass carp), Oreochromis niloticus (Nile tilapia), Piaractus mesopotamicus (Small scaled pacu) and Cyprinus carpio (Common carp) (Claudiano et al., 2019;Yang et al., 2016;Yun et al., 2019). Farmers continuously use different classes of antibiotics, including fluoroquinolones, trimethoprim, β-lactams, ciprofloxacin and OXY to control A. hydrophila infections in aquaculture (Rodgers and Furones, 2009;Suzuki and Hoa, 2012). However, in the past few decades, A. hydrophila has developed resistance to most antibiotics . This prompted us to investigate the differentially expressed genes (DEGs) in A. hydrophila under a sub-lethal concentration of OXY to reveal genes that may have an important involvement in antibiotic resistance. Therefore, in this study, we compared the DEGs among OXY susceptible and resistant strains of A. hydrophila using RNA sequencing to investigate the molecular mechanisms of OXY resistance evolution in A. hydrophila.

Bacterial strains and culture conditions
A. hydrophila ATCC 7966 was cultivated in Luria Bertani (LB) broth and incubated at 30 • C in a shaker at 200 rpm. The minimal inhibitory concentration (MIC) of A. hydrophila ATCC 7966 to OXY was 2.5 μg/mL. Further, the OXY resistant strains were induced from a susceptible strain by antibiotics stress assay. Briefly, the OXY resistant strain was obtained from A. hydrophila ATCC 7966 by ten sequential subcultures in 1/2 MIC concentration of the OXY in LB medium at 30 • C until the MIC increased from the susceptible strain MIC to 4 fold (10 μg/mL) and 8 fold (20 μg/ mL) MIC, respectively as previously described .

RNA extraction, RNA-Seq library construction and sequencing
A single colony from 1 fold, 4 fold and 8 fold MIC bacterial strains were cultured in 5 mL LB medium overnight at 30 • C. Then the bacterial cells were diluted to 5 mL fresh LB medium at a ratio of 1:100 and incubated at 30 • C at 200 rpm until they reached the mid-logarithmic phase (OD 600 reached 1.0). Generally, we have isolated ten members of each of the evolved communities. Among these, we have taken three isolated members for RNA-sequencing. Then, the total RNA was extracted from all samples by TRIzol® reagent according to the manufacturer's instruction (Invitrogen). The contaminating genomic DNA was removed by RNase-free DNase I (Takara Biotechnology, Dalian, China). Further, the RNA quality and quantity were determined by 2100 Bioanalyzer (Agilent) and ND-2000 (NanoDrop Technologies), respectively. After that, the RNA-seq transcriptome library was constructed with the TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA) using 2 μg of total RNA and then the RNA-seq sequencing library was sequenced with the Illumina HiSeq×10 (2 × 150 bp read length) as previously described by Parkhomchuk et al. (2009). The RNA-sequencing analysis was done for three independent experimental samples.

Bioinformatics analysis
The data generated from the Illumina platform were used for bioinformatics analysis using I-Sanger Cloud Platform (www.i-sanger.com) from Shanghai Majorbio Bio-pharm Technology Co., Ltd. Briefly, DESeq2 software was applied for statistical analysis of the three biological replicates. Then, the DEGs were selected based on the ratios between OXY susceptible and resistant strains with |log2FC|≥1 (Fold Change (FC) in log2 scale) and P-adjust value (Student's t-test) at <0.05. Further, the overlapping of the DEGs was analyzed through the Venn diagram by the online tool Venny (http://bioinfogp.cnb.csic.es/too ls/venny/index.html). More, the GO terms of the DEGs were analyzed with the Goatools tool (EMBL, European Molecular Biology Laboratory) and the metabolic pathway of the DEGs were analyzed by the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Li et al., 2018a).

Western blot analysis
To validate the outcome of the transcriptome data, Western blot analysis was performed for the related DEGs as the method previously described (Cheng et al., 2019;Sun et al., 2019). Briefly, the protein samples were separated using 10 % Sodium Dodecyl Sulfate Polyacrylamide Gel Electrophoresis (SDS-PAGE) gels and then transferred to polyvinylidene fluoride (PVDF) membrane for 10 min at 2.5v in a 10 × Tris/Glycine buffer using a Trans-Blot Turbo Transfer System (Bio-Rad, Hercules, CA, USA). After blocking in 5% skim milk with phosphate-buffered saline pH-8.0 containing 0.1 % Tween-20 for 1 h (PBST), the membranes were incubated in a 1:2000 diluted primary antibody at 4 • C for overnight. Then, the membranes were washed five times with PBST and incubated with horseradish peroxidase conjugated secondary antibody in PBST at room temperature for 1 h. More, the expression of the targeted proteins (A0KP37, A0KNL0, A0KNL1, A0KFW5, A0KP35, A0KPP0, A0KGL1, A0KHC0, A0KGL1, A0KHC0, A0KLX3 and A0KP36) were detected by Clarity™ Western ECL Substrate (Bio-Rad) and scanned with the ChemiDoc MP imaging system using Image Lab software (Bio-Rad).

Quantitative reverse transcription-PCR (qRT-PCR) analysis
The DEGs obtained from transcriptome analysis were further validated by qRT-PCR. The primers used in this study were designed by the Primer Premier 5 software and the sequences of the primers were tabulated in Table 1. Moreover, the 16 s rRNA gene was used as the internal control. The total RNA was isolated from OXY susceptible and resistant strains of A. hydrophila using RNAiso Plus (TaKaRa Bio, Tokyo, Japan) based on the manufacturer's protocol. Then, 1 μg of cDNA was synthesized using Prime-Script TM RT reagent kits (Takara Shuzo, Otsu, Japan). Furthermore, the 2× Real Star Green Fast Mixture (GenStar, China) was used as a fluorescent reporter in qRT-PCR analysis. The fluorescence quantitative amplification program consisted of the following conditions with pre-incubation for 2 min at 95 • C, followed by 40 cycles of 15 s at 95 • C, 15 s at 60 • C and the final extended step with the 30 s at 72 • C. Then, the normalized data were calculated by the 2 − ΔΔct method (Livak and Schmittgen, 2001). A total of three independent experimental samples were taken for the qRT-PCR analysis.

Illumina sequencing and quality assessment
A total of three independent biological replicates of susceptible, 4 fold and 8 fold MIC OXY resistant A. hydrophila samples were used for RNA sequencing. Further, a total of 23,826,025, 23,619,297 and 20,191,193 bp of raw reads were obtained for susceptible, 4 fold and 8 fold MIC resistant strains, respectively through RNA sequencing analysis. More, the obtained total genome size and GC content were 4,744,448 bp and 61.55 %, respectively. After filtering and trimming, we have obtained 23,632,948, 23,469,585 and 20,042,380 bp of clean reads for susceptible, 4 fold and 8 fold MIC resistant strains with trim rates of 97.23, 98.22 and 97.25 % (clean Q20), respectively ( Table 2). The obtained data have indicated that a successful transcriptome sequencing of A. hydrophila under a gradient resistance level. Furthermore, the clean reads of DEGs were used for subsequent analysis.

Comparative analysis with the reference genome
Trimmed and clean reads of the A. hydrophila transcriptome were compared with the reference genome of A. hydrophila ATCC 7966 from the National Center for Biotechnology Information (NCBI) GCA_000014805.1 by the Burrows-Wheels software.  (Table 3).

Analysis of DEGs
Each colour in Fig. 1A represents a different sample and the bulk of the expansion represents the region in which the expressions of the nine samples were most concentrated. Further, it is also indicated that the gene expressions in all the samples were at the same level (Fig. 1A). Furthermore, the principal component analysis (PCA) has shown that the three samples are significantly different over the antibiotic resistance evolution and the 4 fold and 8 fold MIC samples have a cross-link ( Fig. 1B). Additionally, the correlation coefficient varied from 0.94 to 1.0, which suggests good data repeatability (Fig. 1C). In this study we identified a total of 4725 genes (Supplement Table 1) which include 4285 mRNA and 440 sRNA. Among a total obtained mRNA, the differentially expressed mRNA genes were 188, in which 69 and 119 genes were up-regulated and down-regulated, respectively. Further, a total of 28 differentially expressed sRNA genes were identified, including 11 upregulated and 17 down-regulated genes (Supplement Table 2). The Venn diagram showed that 6 up-regulated genes such as mtlR, AHA_0550, AHA_0551, AHA_2940, AHA_1316, AHA_2744 and 13 downregulated genes such as AHA_1590, AHA_1589, AHA_3324, AHA_1588, AHA_1586, AHA_1587, AHA_1595, AHA_2856, AHA_3322, AHA_2713, AHA_0810, AHA_2867, AHA_0762 were overlapped in the 4 fold and 8 fold MIC resistant strains. Compared to the 4 fold MIC resistant strain, the 8 fold MIC resistant strain has shown a greater number of DEGs, which suggested that an increasing level of resistance evolution (4 fold MIC to 8 fold MIC) has augmented the number of DEGs (Fig. 1D).
The volcano graph has indicated that the mass of DEGs was extremely significant (|log2FC|≥1, adjusted P < 0.05). The 8 fold MIC resistant samples have shown a total of 185 DEGs when compared to 4 fold MIC resistant samples, which has 22 DEGs with significant differences ( Fig. 2A&B). To further explore the functionally related genes with similar expression patterns, we have selected 26 genes that may participate in the OXY resistance evolution. The 8 fold MIC sample showed the extent of down-regulated genes, mainly involved in sulfur metabolism (metQ, modA, AHA_2577, AHA_3371, AHA_3565,cysA,cysI,cysH,cysC,cysN,cysD,and cysG). On the contrary, the genes related to mannitol metabolism (AHA_0549, AHA_0550, and AHA_0551), sulfur metabolism (hmp and hcp) and efflux pump system (AHA_0021, AHA_0022, and AHA_0023) were down-regulated in the susceptible strain. In contrast, at the high resistance levels (4 fold MIC to 8 fold MIC), all these genes were up-regulated. Overall the obtained data revealed that the gene expressions have significant differences among    the 4 fold and 8 fold MIC resistant strains (Fig. 2C).
We then compared the reference genome with the six databases (NR, Swiss-Prot, Pfam, COG, GO and KEGG) to comprehensively obtain the functional information about DEGs, which is used for functional annotations (Table 4). The unigene annotations in the COG, GO and KEGG databases were about 92.02, 75.68 and 59.91 %, respectively. The transcriptome data were further analyzed by COG classification, which showed 3943 unigenes clustered into 21 functional categories. Among these, the amino acid transport and metabolism (7.15 %), transcription (6.06 %) and inorganic ion transport and metabolism clusters (5.76 %) were identified as major pathways.

GO annotation and KEGG pathway analysis of DEGs
We analyzed the GO enrichment of the DEGs among the susceptible and resistant strains. When compared with the susceptible strain, the 4 fold MIC and 8 fold MIC resistant strains have shown a total of 22 and 185 DEGs, respectively. Furthermore, the several biological processes in the GO terms were over represented compared to the other GO terms. Among the different GO terms, 1883 (58.06 %) genes were related to the biological process. The top 10 biological processes with the highest degree of enrichment are shown in Fig. 3. At 4 fold MIC sample, phosphoenolpyruvate dependent sugar phosphotransferase system, nitrogen cycle metabolic process and carbohydrate transmembrane transport were enriched in two genes. Meanwhile, the other seven biological processes, including polyol transport, organic hydroxy compound transport, molybdate ion transport, mannitol transport, mannitol metabolic process, hexitol metabolic process and denitrification pathway were enriched in one gene (Fig. 3A). Furthermore, as the degree of resistance increased to 8-fold MIC, the three processes, such as oxidation-reduction process (43 genes), cellular respiration (8 genes) and ATP metabolic process (7 genes), were enriched with most of the genes. Moreover, the three processes of sulfur metabolism (sulfate assimilation, hydrogen sulfide metabolic process and hydrogen sulfide biosynthetic), respiratory electron transport chain, oxidative phosphorylation and ATP synthesis coupled electron transport process were enriched in six genes. Additionally, the acetyl-CoA metabolic process was enriched in three genes (ack-2A, acsA and sucB) (Fig. 3B).
We further analyzed these DEGs by KEGG pathways analysis. A total of 1494 genes were mapped to 64 pathways. Many genes related to multiple pathways were enriched in both 4 fold and 8 fold MIC resistant strains. The five pathways, including nitrogen metabolism, pyruvate metabolism, microbial metabolism in diverse environments, citrate cycle (TCA cycle) and propanoate metabolism were enriched in both the 4 fold and 8 fold MIC resistant strains. Interestingly, over the increasing resistance evolution (4 fold MIC to 8 fold MIC), the genes related to sulfur metabolism (cysD, cysN, cysC, cysH, cysI, cysA, AHA_3371 and AHA_2577) pathway were significantly enriched. Further, the results of KEGG are consistent with the results of GO enrichment analysis. More, the results of the GO and KEGG analyses have revealed that the OXY resistance affects multiple biological functions in A. hydrophila, mainly in energy biogenesis such as sulfur metabolism and mannitol metabolism. So, the overall obtained transcriptomic data clearly have indicated that the genes related to sulfur metabolism, mannitol metabolism and efflux pump system may play an important role in OXY resistance in A. hydrophila when the concentration is increased (Fig. 4A&B).

Validation of transcriptome data at the protein level
We have selected several proteins related to DEGs such as A0KP37, A0KNL0, A0KNL1, A0KFW5, A0KP35, A0KPP0, A0KGL1, A0KHC0, A0KLX3 and A0KP36 for Western blot analysis to validate the outcome of the RNA sequencing data. The obtained results showed that the cysG (A0KP37), cysH (A0KNL1), cysA (A0KFW5), cysN (A0KP35), AHA_3793 (A0KPP0), groL (A0KGL1), AHA_1130 (A0KHC0), AHA_2766 (A0KLX3) and cysD (A0KP36) proteins were substantially down-regulated in both the 4 fold and 8 fold MIC samples compared to the susceptible strain (Fig. 5). In general, many of the protein expression tendencies were consistent with the transcriptome data. However, not all the protein levels are consistent with the RNA levels. For example, some proteins such as A0KP35, A0KPP0 and A0KLX3 were decreased in the 8 fold MIC sample but increased slightly in the 4 fold MIC sample. On the other hand, the protein A0KNL0 was decreased in the 4 fold MIC sample but increased in the 8 fold MIC sample. However, the genes related to all of these selected proteins were substantially decreased at the RNA level in Table 4 Statistical results of the gene functional annotation.

Validation of transcriptome data at the mRNA level
Finally, to further validate the transcriptome data, the qRT-PCR analysis was performed for a total of 20 candidate genes related to the efflux pump system, mannitol and sulfur metabolism. Among 14 sulfur metabolism genes, 2 genes such as hcp and hmp were up-regulated and 12 genes such as cysA, cysC, cysD, cysH, cysI, cysN, cysG, AHA_3371, modA, metQ, AHA_2577 and AHA_3565 were down-regulated in the 8 fold MIC resistant strain. Further, the mannitol metabolism related genes (AHA_0549, AHA_0550, and AHA_0551) and efflux pump system related genes (AHA_0021, AHA_0022, and AHA_0023) were up-regulated in the 8 fold MIC resistant strain. More, the correlation analysis showed a strong positive relationship between the experiment and data (R 2 = 0.803), which made the validation of transcriptome data by qRT-PCR reliable and repeatable (Fig. 6).

Discussion
The transcriptome is the comprehensive collection of expressed RNA transcripts in a cell. Its portrayal is vital for decoding the functional intricacy of the genome and understanding the cellular activities in organisms, including development, growth, pathogenicity and immune defence (Xiang et al., 2010;Zhang et al., 2018;Tang et al., 2020). Microarray and expressed sequence tag analyses have long been used to understand the molecular mechanisms underlying antibiotic resistance in bacterial pathogens (Dally et al., 2013). However, Next-Generation Sequencing (NGS) platforms such as Illumina HiSeq×10 are better suited to quantifying transcripts' expression at low levels than microarrays and EST analyses (Rao et al., 2015) because this revolutionary sequencing technique confirms the direct transcript profiling without any compromise.
This study determined the comparative transcriptome differences between susceptible and resistant strains of A. hydrophila using the Illumina Hiseq sequencing platform. The obtained data revealed the   (Zhu et al., 2017) have reported the successful transcriptome sequencing of enrofloxacin susceptible and enrofloxacin resistant A. hydrophila strains with a high percentage of clean reads. Similarly, in this study, the obtained transcriptomic data have shown a high percentage of clean reads in OXY susceptible and resistant A. hydrophila strains, compared to the reference genome. Therefore, it represents the quality of transcriptomic sequencing data that met the demand for subsequent studies.
GO annotation and KEGG pathway enrichment analysis were done to reveal the biological function of the significant DEGs between the OXY susceptible and resistant A. hydrophila strains. The outcome of GO analysis revealed that at the 4 fold MIC degree of resistance, polyol transport, nitrogen cycle metabolic process, mannitol transport, mannitol metabolic process and denitrification pathway were highly enriched. In 8 fold MIC, the processes such as sulfate assimilation, oxidation-reduction process, hydrogen sulfide metabolic process and hydrogen sulfide biosynthetic process were highly enriched. The outcomes of the KEGG pathway analysis have exposed that a significant percentage of genes were highly enriched in the phosphotransferase system (PTS), nitrogen metabolism and fructose and mannose metabolism in the 4 fold MIC OXY resistant A. hydrophila strain. Also, a substantial percentage of genes were highly enriched in beta-lactam resistance, sulfur metabolism, microbial metabolism in diverse environments and citrate cycle (TCA cycle) in 8 fold MIC OXY resistant A. hydrophila strain. Overall, the obtained data from GO and KEGG analyses have revealed that most of the genes were highly enriched in sulfur and mannitol metabolisms, which also supported the outcome of Heat map data. We have selected some DEGs related to sulfur metabolism, mannitol metabolism, and efflux pump system for further analyses based on the obtained data.
To validate the outcome of transcriptome sequencing results, we have performed Western blot and qRT-PCR analyses for the candidate proteins or genes related to the sulfur metabolism, mannitol metabolism and efflux pump system. (Rusniok et al., 2009) have reported that the five proteins encoded by the cysD, cysH, cysI, cysJ and cysN genes have the ability to reduce the sulfate (SO 4 2− ) to hydrogen sulfide (H 2 S) in Neisseria meningitidis for their sulfur source. In line with this experimental result, we have observed the down-regulation of cysD, cysH, cysI and cysN in 4 fold and 8 fold MIC OXY resistant A. hydrophila strains in qRT-PCR analysis. This indicates that A. hydrophila may utilize the sulfate as a sole sulfur source by this cysDHIN gene system. Therefore, the down-regulation of these genes in A. hydrophila on OXY treatment may affect the energy metabolism pathway. Furthermore, the expressions of 6 proteins [Siroheme synthase 2 (A0KP37), Sulfite reductase [NADPH] hemoprotein beta-component (A0KNL0), Phosphoadenosine phosphosulfate reductase (A0KNL1), Sulfate/thiosulfate import ATP-binding protein (A0KFW5), Sulfate adenylyltransferase subunit 1 (A0KP35), and Sulfate adenylyltransferase subunit 2 (A0KP36)], which are involved in sulfur metabolism among OXY susceptible and resistant strains were validated by Western blot analysis. The results have shown that the A0KP37, A0KNL1, A0KFW5, A0KP35 and A0KP36 proteins were substantially down-regulated in the 8 fold MIC OXY resistant A. hydrophila strain when compared with their respective controls. The sulfur metabolism is considered one of the important energy metabolism pathways in several bacteria and a drug target for potential bactericidal therapy. In our previous proteomic study, we have also observed that the increasing concentration of OXY substantially downregulated the energy metabolism pathway related genes such as R4VS58, R4VXT4 (TCA cycle), R4VW15 (glycolysis) and A0KFS2 (pyruvate metabolic pathway) in A. hydrophila . Therefore, the obtained data in this study further confirmed our previous hypothesis that decreasing energy generation might be a common phenomenon for A. hydrophila responses to OXY. Moreover, several recent studies have also reported that the down-regulation of the energy metabolism pathway may help bacterial survival in the antibiotics environment (Cheng et al., 2019;Palde et al., 2016), which indicates that the decreasing level of sulfur metabolism may protect A. hydrophila from a threatening environment. Two of the foremost weapons of the innate immune system to eradicate pathogenic organisms are the production of reactive nitrogen species (RNS) and reactive oxygen species (ROS), which are resultant from the nitric oxide and superoxide produced by the nitric oxide synthase and NADPH oxidase, respectively (Figueiredo et al., 2013). The hybrid cluster proteins (HCPs) are a family of bacterial proteins which protect microorganisms from RNS and ROS mediated toxicity. The research work of (Yurkiw et al., 2012) have revealed that HCPs are required to sustain the high amounts of nitrite reduction by the nitrite reductase in Desulfovibrio vulgaris; they also proposed the role of HCPs in protection from nitrite derived products. Flavohemoglobins (flavoHbs) are nitric oxide (NO) detoxifying proteins that metabolize the NO to nitrate and protect the microorganisms from growth inhibition and killing by NO producing immune cells. Furthermore, flavoHbs has a protective role against RNS and ROS and have the capability to protect microorganisms against nitrosative stress (Helmick et al., 2005;Koskenkorva-Frank and Kallio, 2009). The transcriptome sequencing and qRT-PCR analyses data have revealed that the hcp and hmp genes were up-regulated, suggesting that A. hydrophila protect themselves from various nitrosative stresses when treated with increasing concentrations of OXY.
ModA is a soluble periplasmic protein that is part of an ABC transporter system to uptake nutrients (Tam and Saier, 1993). Besides its functions in nutrients uptake, ModA was found to be down-regulated when responding to cell-wall-active antibiotics in Staphylococcus aureus (Utaida et al., 2003). Whole genome sequencing among vancomycin resistant and susceptible S. aureus strains also found loss of its functional mutant in the vancomycin resistant strain (Avison et al., 2002). In line with these reports, the decreased expression level of modA was observed in the OXY resistant A. hydrophila strain. This indicates that modA may play an important role in the antibiotic resistance evolution and perhaps are involved in antibiotic resistance by reducing Fig. 6. Validation of transcriptomic data by qRT-PCR. The differential expression of sulfur metabolism (Round), mannitol metabolism (Triangle) and efflux pump system (Square) related genes were assessed by qRT-PCR at the mRNA level. Of these, 12 genes (cysA, cysN, cysC, cysD, cysI, cysH, cysG, metQ, modA, AHA_2577, AHA_3371 and AHA_3565) related to sulfur metabolism were down-regulated, whereas 2 genes (hcp and hmp) related to sulfur metabolism, 3 genes (AHA_0549, AHA_0550 and AHA_0551) related to mannitol metabolism and 3 genes (AHA_0021, AHA_0022 and AHA_0023) related to the efflux pump system were up-regulated. Green indicates down-regulation, while red indicates up-regulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) membrane permeability.
Mannitol is one of the most abundant of six carbon polyols in microorganisms, which plays a vital role in stress tolerance, osmoregulation and antioxidant defence in them (Nguyen et al., 2019). Kumar et al. have elucidated the functional properties of the mannitol operon system and revealed their role on the virulence and pathogenicity of V. cholerae (Kumar et al., 2011). In this study, the obtained transcriptome data have shown the normal regulation or down-regulation of mannitol metabolism related genes at the 1 fold MIC concentration of OXY in A. hydrophila. In contrast, the increasing concentration of OXY showed the up-regulation of mannitol metabolism related genes for their protection against various oxidative stresses caused by high concentrations of antibiotic exposure. The obtained qRT-PCR data also validated the same, in which the three mannitol related genes such as AHA_0549, AHA_0550 and AHA_0551 were up-regulated in the 8-fold MIC OXY resistant A. hydrophila strain.
Acriflavin is a multidrug efflux pump system involved in the exclusion of antibiotics from bacterial cells. It also probably protects the bacterial cells from hydrophobic inhibitors. In addition, it is a homolog to the AheABC efflux pump, which belongs to the resistance nodulation cell division family. The AheABC efflux pump is associated with the intrinsic resistance to several antibiotics and compounds. (Cruz et al., 2015) and (Hernould et al., 2008) have identified that this AheABC efflux pump system extruded at least 9 antibiotics and 13 substrates in A. hydrophila ATCC 7966, thus indicating their important role in antibiotics resistance. Our current data further provides evidence that increasing AheB in both 4 fold and 8 fold MIC OXY resistant strains may promote the evolution of antibiotics resistance.
Furthermore, we have observed the up-regulation of AHA_0022 (RND transporter) in 4 fold and 8-fold MIC OXY resistant A. hydrophila strains. The RND transporter efflux pumps are common among bacterial pathogens. They afford a second-line barricade against the different classes of antibiotics and have a three-way assembly such as a periplasmic adaptor protein, inner membrane spanning transporter protein and outer membrane efflux protein. These three components help the bacteria force out the antibiotics promptly to the extracellular milieu through the water-filled channel produced by the three-way assembly of RND transporter efflux pumps (Sandhu and Akhter, 2018). Our data suggest that AHA_0022 may contribute to antibiotic resistance by pumping the OXY in A. hydrophila.
Besides these, AHA_0023 (coded as a TetR family transcriptional regulator) was previously reported to be involved in bacterial antibiotic resistance. It plays a major role in metabolism, quorum sensing and antibiotic resistance (Cuthbertson and Nodwell, 2013). In this study, we have observed the increased expression level of the AHA_0023 gene during the OXY resistance evolution. The up-regulation of the efflux pump system related gene cluster (from AHA_0021 to AHA_0023) indicates that they may force out the OXY or regulate the various metabolism and virulence factors.

Conclusion
Based on the Illumina HiSeq×10 RNA sequencing analysis, the resistance to OXY in A. hydrophila was identified to be more complex as it targets multiple metabolic pathways. Moreover, the sulfur metabolism related protein (sulfate adenylyltransferase subunit 1, sulfate adenylyltransferase subunit 2, hybrid cluster protein, flavohemoprotein, phosphoadenosine phosphosulfate reductase, sulfite reductase (NADPH) hemoprotein beta-component, molybdate ABC transporter periplasmic molybdate-binding protein), mannitol metabolism related proteins (mannitol operon repressor, mannitol-1-phosphate 5-dehydrogenase, PTS system mannitol-specific EIICBA component) and efflux pump system related proteins (acriflavin resistance plasma membrane protein, RND transporter protein and TetR family transcriptional regulator protein) may play significant roles in the OXY resistance of A. hydrophila. Furthermore, it is speculated that the mechanism of OXY resistance in A. hydrophila is closely related to a decrease of various nitrosative stresses, sulfur metabolism and increase of mannitol utilization system, various efflux pump systems such as acriflavin multidrug efflux system and RND transporter efflux system activation. Overall, these data provide a deeper insight to understand the intrinsic molecular mechanisms of OXY resistance evolution in A. hydrophila.

Author Contributions
All the authors contributed extensively to work presented in this manuscript. XL and HY designed the experiments; JY, SR, LC and FZ generated experimental data and wrote the manuscript. SR, XL, YZ, LL and SM conceived the work and critically reviewed the manuscript. The authors declare no competing financial interest.

Declaration of Competing Interest
All the authors declare no financial conflict of interest and have fulfilled the authorship criteria for the manuscript.