Transcriptome differences between enrofloxacin-resistant and enrofloxacin-susceptible strains of Aeromonas hydrophila

Enrofloxacin is the most commonly used antibiotic to control diseases in aquatic animals caused by A. hydrophila. This study conducted de novo transcriptome sequencing and compared the global transcriptomes of enrofloxacin-resistant and enrofloxacin-susceptible strains. We got a total of 4,714 unigenes were assembled. Of these, 4,122 were annotated. A total of 3,280 unigenes were assigned to GO, 3,388 unigenes were classified into Cluster of Orthologous Groups of proteins (COG) using BLAST and BLAST2GO software, and 2,568 were mapped onto pathways using the Kyoto Encyclopedia of Gene and Genomes Pathway database. Furthermore, 218 unigenes were deemed to be DEGs. After enrofloxacin treatment, 135 genes were upregulated and 83 genes were downregulated. The GO terms biological process (126 genes) and metabolic process (136 genes) were the most enriched, and the terms for protein folding, response to stress, and SOS response were also significantly enriched. This study identified enrofloxacin treatment affects multiple biological functions of A. hydrophila. Enrofloxacin resistance in A. hydrophila is closely related to the reduction of intracellular drug accumulation caused by ABC transporters and increased expression of topoisomerase IV.


Introduction
Following the decline in the capture fishing industry and diminishing wild fish stocks, the aquaculture industry has become an important source of food fish [1]. However, bacterial diseases hinder desirable production outputs. The gram-negative bacterium Aeromonas hydrophila is one of the major causative agents of disease and can cause serious damage in many animals [2], especially fish [3][4] as well as humans [5]. A. hydrophila, which is a representative of the Aeromonadaceae family, is an emerging aquatic pathogen that is distributed in a wide variety of aquatic systems [6][7]. It primarily inhabits freshwater and the intestines of freshwater animals. Farmers use a wide range of antibiotics or chemicals to control A. hydrophila infection [8]. Enrofloxacin is a third-generation fluoroquinolone with a broad antibacterial spectrum and high potency that is commonly used to treat bacterial infections afflicting aquaculture [9]. Enrofloxacin has been used primarily to control A. hydrophila infections in aquaculture, and recently, A. hydrophila has developed strong resistance to enrofloxacin among other drugs [10]. This resistance has rendered it increasingly difficult to treat diseases caused by A. hydrophila in aquaculture animals. Moreover, heavy antibiotic use is associated with negative effects such as antibiotic resistance in the environment and fish [11]. Thus, eventually, antibiotic use may be detrimental to the environment and human health.
Based on previous reports, quinolone-resistant bacteria adopt the following three main strategies of antibiotic resistance: First, chromosome-mediated changes in the topoisomerase target sites (changes in the amino acids in the quinolone resistance-determining region) [12]; second, reduction in intracellular drug accumulation caused by efflux pump [13]; and third, bacterial protection conferred by plasmid-encoded qnr protein [14]. However, the mechanisms by which A. hydrophila is resistant to enrofloxacin remains unclear, and little is known about its molecular mechanisms of resistance. Here, we conducted de novo transcriptome sequencing for the comprehensive analysis of the global transcriptomes of enrofloxacin-susceptible and enrofloxacin-resistant strains. Transcriptomic profiling is used to analyze gene expression and signaling pathways in specific tissues or cells. The recent rapid development of next-generation sequencing technologies such as the Solexa/Illumina technology offers great advantages in analyzing the functional complexity of the entire transcriptome [15]. Next-generation sequencing techniques have been used for transcriptome analyses to simultaneously provide data on sequence polymorphisms and the levels of gene expression involved in cellular development, cancer, and immune responses [16][17].
In the present study, we examined the genetic diversity of A. hydrophila using de novo transcriptome sequencing and investigated the molecular mechanisms of enrofloxacin resistance in A. hydrophila.

Culture of susceptible strain and induction of an enrofloxacin-resistant strain
A. hydrophila ATCC 7966, which was used in the present study, had been maintained at the National Pathogen Collection Center for Aquatic Animals, China. Because the genome background of the ATCC 7966 standard strain is known, it is known to be sensitive to quinolones; therefore, it was selected as the sensitive strain. This strain was inoculated on Luria Broth (LB) agar and incubated at 28˚C for 24 h.
The type strain ATCC 7966 was used to develop an enrofloxacin-resistant strain by culturing it in the presence of gradually increasing concentrations of enrofloxacin in vitro, and enrofloxacin-resistant strain 7966QR was judged by the critical concentration suggested by the Clinical and Laboratory Standards Institute (CLSI, 2011). The results of the drug sensitivity test were determined using the disc diffusion method. Resistance and sensitivity to enrofloxacin were determined on the basis of the drug sensitivity evaluation criteria issued by the American Association of Clinical Laboratory Standards in 2011 [18]. Briefly, ATCC 7966 cells were inoculated on LB agar containing 1/2 the minimum inhibitory concentration (MIC) of enrofloxacin (Shanghai Guoyao Chemical Reagent Co. Shanghai, China). Every 2 days, the cultured strains were inoculated into fresh LB agar containing 2× the previous concentration of enrofloxacin. The strains were inoculated to LB agar that containing no enrofloxacin until the MIC was conformed to the drug resistance determination. The resultant resistant strains were screened on LB for 12 generations until the resultant strain (7966QR) could be deemed resistant.
Suspensions containing the susceptible (ATCC 7966) and resistant (7966QR) strains were collected and centrifuged, and the final products were store at 4˚C until analysis.
RNA isolation, RNA-Seq library construction, and sequencing Frozen samples were ground in a mortar with liquid nitrogen, and total RNA was extracted from approximately 1 mL of bacterial suspension with TRIzol reagent (Invitrogen, USA), according to the manufacturer's instructions. DNA contaminants were removed by treatment with RNase-free DNase I (Takara Biotechnology, Dalian, China). The resultant total RNA was dissolved in 200 μL RNase-free water. The concentration of total RNA was determined using a Nano-Drop2000 spectrophotometer (Thermo Scientific, USA), and its integrity was checked using an RNA 6000 Pico LabChip with the Agilent 2100 bioanalyzer (Agilent, USA) at 37˚C for 1 h, and the sample volume was diluted to 250 μL using nuclease-free water. Messenger RNA (mRNA) was further purified with a Micropoly (A) Purist kit (Ambion, USA) according to the manufacturer's protocol. mRNA was dissolved in 100 μL RNA Storage Solution (Ambion), purified using oligo (dT) magnetic beads, fragmented by treatment with divalent cations and heat, and reverse transcribed into cDNA using reverse transcriptase and random hexamer primers. This was followed by second-strand cDNA synthesis using DNA polymerase I and RNaseH. The resultant double-stranded cDNA was end-repaired using T4 DNA polymerase, Klenow fragments, and T4 polynucleotide kinase followed by a single (A) base addition using Klenow 3 0 to 5 0 exopolymerase. This was then ligated with an adapter or index adapter using T4 Quick DNA ligase. The size range of the adapter-modified fragments was selected by gel purification, and these sizes were used as templates in PCR amplification. The cDNA library was validated with an Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-time PCR system and sequenced on a flow cell using an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA).

Sequencing, data processing, and quality control
We filtered low-quality reads and removed 3 0 -adapter sequences using Trim Galore. The obtained reads were cleaned using FastQC software (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/), and the content and quality of the nucleotide bases in the sequencing data were evaluated. Next, we conducted a comparative analysis with the reference genome (Aeromonas hydrophila subsp. Hydrophila ATCC 7966). For each sample, sequence alignment with the reference genome sequences was carried out using Tophat [19].

Assembly and functional annotation
High-quality reads were obtained after removing the adapter sequence, low-quality reads (reads with ambiguous bases N), and duplicate sequences using Trim Galore and FastQC, and then FastQC software was used to clean reads and evaluate the performance of different kmers. Next, the clean reads were combined using de Bruijn graphs and SOAPdenovo software based on sequence overlap to form longer fragments (without ambiguous 'N' reads), to create contigs [20]. Furthermore, the contigs were connected into transcript sequences and joined into scaffolds using paired-end reads. The paired-end reads were also used to fill the gaps in scaffolds, where the unigenes have the least Ns and cannot be extended on both ends. Based on the results of the assembly evaluation, the best results were selected and used for clustering analysis using TGI Clustering tools to achieve a unigene database [21]. The obtained unigenes were compared with the National Center for Biotechnology Information (NCBI), non-redundant protein (Nr), and UniProt databases using BLASTx (Basic Local Alignment Search Tool) search with an E value < 0.00001. Based on the results of the Nr annotation, we used Blast2GO software (https://www.blast2go.com/) to analyze functional annotations by gene ontology terms (GO; http://www.geneontology.org) [22]. The unigenes were also aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Group (COG), and Swiss-Prot databases to predict and classify gene functions to perform pathway annotation searching for unigenes with similarity >30% and an E value <0.00001, and all the information was merged.

Analysis of differentially expressed unigenes
To estimate the expression level (relative abundance) of a specific transcript expressed as fragments per kilobase per million fragments mapped (FPKM), we used RSEM software with the default parameter settings [23]. The expression level of each transcript was transformed using base 2 log 2 (FPKM+1). The fold changes in the expression of a transcript and differentially expressed gene (DEG) were estimated using DESeq software [24]. Two-fold changes in expression level and differences with a p value of <0.05 were considered significant.

GO functional and pathway enrichment analysis of DEGs
We annotated DEGs to analyze the transcriptome differences between enrofloxacin-resistant and enrofloxacin-susceptible strains of A. hydrophila. To this end, we used GO terms in accordance with previously published procedures [25]. This p = analysis first mapped all DEGs to GO terms in the database by calculating gene numbers for every term followed by an ultrageometric test to identify significantly enriched GO terms in DEGs compared to the transcriptome background. The following formula was used: where N represents the number of all genes with GO annotation, n represents the number of DEGs in N, M represents the number of all genes annotated to specific GO terms, and m represents the number of DEGs in M. The calculated p value was subjected to Bonferroni correction. A corrected p value <0.05 was defined as the "threshold." GO terms were considered significantly enriched in the DEGs.

Pathway analysis of DEGs
Pathways of DEGs were annotated against the KEGG database using the BLASTall program (http://nebc.nox.ac.uk/bioinformatics/docs/blastall.html). Enriched DEG pathways were identified according to the same formula as that used in the GO analysis. In this case, N represented the number of all genes with KEGG annotations, n represented the number of DEGs in N, M was the number of all genes annotated to specific pathways, and m was the number of DEGs in M [25].

Verification of DEGs using qRT-PCR
Quantitative RT-PCR (qRT-PCR) was used to verify the expression levels of DEGs that were identified by RNA-Seq analysis. Primers were designed using Primer 5 software and SpTub-b was used as the reference gene [25,26]. Reactions were performed in a 25-μL reaction volume composed of 2 μL cDNA, 0.5 μL each of forward and reverse primers (10 μM), 12.5 μL SYBR Premix Ex Taq (2×), and 9.5 μL RNase-free H 2 O. The thermal cycle protocol was as follows: 95˚C for 30 s followed by 40 cycles of 95˚C for 5 s, 60˚C for 30 s, and 72˚C for 30 s. Melting curve analysis was performed at the end of qRT-PCR to confirm PCR specificity.

Results
Illumina sequencing and quality assessment  Table 2). We then compared the unigenes of the sample species with the common data genes, and functional annotation was performed based on the similarity of the genes. The protein sequences were compared with the KOG, GO, and KEGG databases. The annotation of unigenes in the Swiss-Prot and TrEMBL databases accounted for 79.77% and 99.93% of the total unigenes, respectively (Table 3 and Fig 1).
The unigene annotations in the COG, GO, and KEGG databases were about 82.19%, 79.57%, and 62.3%, respectively (Table 3). Transcripts were analyzed by COG classification. There were 3,388 unigenes clustered into 25 functional categories (S1 Fig). The "amino acid transport and metabolism" and "signal transduction mechanisms" clusters represented the majority of transcripts (276 transcripts, 8.15%; S1 Fig). GO and KEGG database analysis of unigenes revealed that most unigenes were enriched in cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems (S2 and S3 Figs).

Analysis of DEGs
The Cuffdiff program was used to generate A. hydrophila gene expression profiles to identify genes that are differentially expressed between the resistant and susceptible strains of A. hydrophila (Figs 2 and 3). The program identified that among the DEGs, 135 genes were markedly upregulated and 83 were markedly downregulated, indicating that the gene expression had changed in the drug-resistant strains.

GO annotation of DEGs
We used GO to determine the biological functions in which the DEGs are involved. GO functional enrichment analysis also involved cluster analysis of expression patterns. Thus, the expression patterns of DEGs annotated with a given GO term were easily obtained. All the annotated genes were classified into three GO domains: biological process, cellular component, and molecular function. Dissimilar expression profiles of the DEGs in the treated and control groups were used to determine the molecular mechanisms of enrofloxacin resistance in A. hydrophila. The expression profiles of the three GO domains were as follows (Fig 4 and S1 Table). Biological process: formate metabolic process (5 genes), histidine catabolic process to glutamate and formate (4 genes), amine metabolic process (7 genes), histidine catabolic process to glutamate and formamide (4 genes), formamide metabolic process (4 genes), histidine catabolic process (4 genes). Cellular component: HslUV protease complex (2 genes), proteasome complex (2 genes), bacterial-type flagellum basal body, distal rod (2 genes), cytosolic proteasome complex (2 genes). Molecular function: oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (3 genes); anion transmembrane-transporting ATPase activity (3 genes); oxidoreductase activity, acting on single donors with incorporation of molecular oxygen (3 genes); oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation (3 genes); carbon-nitrogen lyase activity (3 genes); ATPase activity, coupled to transmembrane movement of ions (5 genes); and dioxygenase activity (3 genes).

KEGG pathway analysis of DEGs
To further explore the biological functions of the DEGs, DEGs were mapped to the KEGG database and enriched to important pathways based on the whole transcriptome. A total of 218 genes were mapped to 67 pathways. Many of the genes were found in multiple pathways, whereas many others were restricted to a single pathway. These pathways included metabolism, genetic Enrofloxacin-resistant and enrofloxacin-susceptible strains of Aeromonas hydrophila information processing, cellular processes, organismal systems, and environmental information processing. The metabolism-related pathways were the most significantly (Fig 6 and  S2 Table).
The expression of certain genes such as AHA_2490 (PI3K-Akt signaling pathway); AHA_2490 (NOD-like receptor signaling pathway); AHA_0608, AHA_2812, AHA_0913, AHA_3728, AHA_1964, AHA_1687, AHA_4285, AHA_1595, and AHA_2813 (ABC transporters) were upregulated, while the AHA_1419 (environmental information processing), Enrofloxacin-resistant and enrofloxacin-susceptible strains of Aeromonas hydrophila AHA_1331 (drug metabolism-cytochrome P450), AHA_2360, AHA_1331 (glycolysis/gluconeogenesis), and AHA_1331 (metabolism of xenobiotics by cytochrome P450) genes were downregulated. The ABC transporter genes were expressed at higher levels than drug metabolism-cytochrome P450, indicating that A. hydrophila resistance to enrofloxacin may be mediated by a mechanism involving ABC transporters. These finding are consistent with the results of the GO enrichment analysis. Enrofloxacin-resistant and enrofloxacin-susceptible strains of Aeromonas hydrophila Overall, the results of the DEG pathway analysis support the viewpoint that enrofloxacin inhibits A. hydrophila growth by affecting multiple biological functions, such as energy biogenesis, protein synthesis, and metabolism.

Verification of the differential expression of DEGs
Based on the results of the GO and KEGG analyses, the primers of eight genes significantly differed between the resistant and susceptible A. hydrophila strains and were therefore considered to be related to drug metabolism. With clear functional implication, these primers were designed to verify the expression of these DEGs identified in the RNA-Seq analysis. All primer sequences are listed in Table 4.

Discussion
Transcriptome sequencing is a powerful technique for studying the mechanisms of changes in biological characteristics of an organism and has been used successfully in some species [27][28][29]. In our study, we examined the transcriptome of A. hydrophila using the Illumina sequencing platform and explored the molecular mechanism of enrofloxacin resistance in A. hydrophilia. Compared with the reference genome, the total mapped rates of reads were 94.19% and 93.29% in the A. hydrophilia transcriptome of the enrofloxacin-susceptible (ATCC 7966) and enrofloxacin-resistant (ATCC 7966QR) strains, respectively, indicating that the quality of sequencing data met the demand for follow-up studies.
We obtained 218 DEGs and classified them into 1,052 GO terms consisting of three domains: biological process, cellular component, and molecular function. Of these, 176 GO terms were found to have dramatic changes in expression. We mapped the DEGs to 68 pathways, of which 10 were significantly enriched. We divided the genes into five branches based on the KEGG metabolic pathway involved: cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems. The metabolism-related biological pathways and biosynthesis of amino acids were the most significantly enriched pathways in this analysis; these pathways are responsible for the main biological functions of A. hydrophila. The results of the KEGG pathway analysis revealed that a considerable percentage of genes was enriched in ABC transporters, metabolism of xenobiotics by cytochrome P450, and drug metabolism-cytochrome P450. All genes enriched in ABC transporters were upregulated, whereas all the genes enriched in metabolism of xenobiotics by cytochrome P450 and drug metabolism-cytochrome P450 were downregulated. We speculate that Enrofloxacin-resistant and enrofloxacin-susceptible strains of Aeromonas hydrophila enrofloxacin promotes the protein expression of ABC transporters and inhibits the protein expression of cytochrome P450. This may be related to enrofloxacin resistance of A. hydrophila. Our study supplements the previous studies by Seshadri et al. who investigated the virulence genes in A. hydrophila by sequencing its genome [30]. The results of the GO annotation enrichment revealed that most DEGs were related to transmembrane transport and biosynthesis and degradation of amino acids. Furthermore, some of the genes were mapped to cellular response to DNA damage stimulus, and among the DEGs, gyrA was upregulated, which was verified by qRT-PCR. The qRT-PCR results were generally consistent with the results of the transcriptome analysis. Our results are in agreement with the results of Shakir et al., who reported that changes in the topoisomerase target sites of chromosomes (amino acid changes in quinolone resistance-determining regions, QRDRs) could produce drug resistance and, A. hydrophila strains exhibiting high levels of resistance to antibiotics are common in the strains with gyrA and parD double mutations in QRDRs [5]. Together, these results indicate that A. hydrophilia resistance to enrofloxacin occurs primarily due to alterations in multiple biological functions, energy biogenesis, protein synthesis, and metabolism. Our findings support that the mechanism of enrofloxacin resistance in A. hydrophila is closely related to reduction of intracellular drug accumulation caused by ABC transporters and topoisomerase IV. The data revealed that the upregulation or downregulation of these six genes was consistent with the RNA-Seq results. Together, these results indicate that the qRT-PCR and RNA-Seq results were reliable overall; however, further studies to determine the molecular mechanisms of resistance to enrofloxacin are required (Fig 7). https://doi.org/10.1371/journal.pone.0179549.t004 Most of the genes in A. hydrophilia encode putative proteins whose functions are not clear. In order to facilitate further research, we selected 8 of the DGEs whose functions were clearly known and subjected them to qPCR analysis to verify the RNA-Seq results of. The qPCR results were consistent with the results of the transcriptome analysis, except for two genes. Therefore, we believe that the qualities of the A. hydrophila transcriptomes are adequate for further studies on its functional genes. These findings greatly extend the existing sequence resources relating to A. hydrophilia and provide abundant genetic information that can be applied to further understand the molecular mechanisms of enrofloxacin resistance in A. hydrophilia in aquaculture.
Supporting information S1 Table.