From mid-April 2022 to the beginning of May 2022, we investigated a case of severe mortality in juvenile size Labeo rohita (Length = 120.4 ± 6.8 mm, Weight = 12.51 ± 1.9 g) in the aquaculture farms of Moina, Purba Medinipur district (22^o18’33 N, 87^o50’20 E), West Bengal, India. In the aquaculture farms, mixed fish species, including L. rohita, Cirrhinus mrigala, L. catla, and Hypopthalmichthys molitrix were cultured. As the central nodal agency for fish disease surveillance, the team of researchers from ICAR-Central Inland Fisheries Research Institute, Kolkata, India visited the culture facility. During our visit, we estimated approximately 60% mortality in L. rohita through passive data collection. The fish were showing signs of disease, including extreme lethargy, hemorrhage, ulcer, discoloration, and redness in fins and all over the body surface. Suspecting a microbial infection, moribund juveniles of L. rohita showing clinical signs were transported to the laboratory, necropsied, and subjected to further microbiological tests.

Diseased rohu fish (Labeo rohita) showing clinical signs, including ulcerated skin, muscle lesions, and tail and fin rot, are collected and sacrificed for bacterial isolation. The fish were anaesthetized using clove oil (Merck, Germany) at a 50 μL/L dose, and the gut sample was collected aseptically from the fish and transferred to Nutrient Broth (NB) (Himedia, India) and incubated for 24 h at 37°C at 200 rpm. The growth of morphologically similar colonies on Nutrient Agar (NA) plates indicates the presence of a single bacterial strain. A single colony was streaked onto freshly prepared NA plates and incubated at 37°C for 24 h to obtain a pure strain. Subsequently, the colony was transferred to NB, and after 24 h of incubation at 37°C, the bacterial strain stock culture was prepared in 30% glycerol and stored at-20°C until further use.

Organization for Economic Cooperation and Development (OECD) guidelines were followed for the handling and care of experimental animals. The animal utilization protocol was approved by Institutional Animal Ethics Committee, ICAR-Central Inland Fisheries Research Institute, Kolkata, India, (IAEC/2021/04) for the experimental setup.

Healthy L. rohita (Length = 114.2 ± 2.5 mm, Weight = 12.1 ± 5.2 g) were purchased from a local fish hatchery. The fish appeared normal and healthy with no external clinical symptoms like ulcer, hemorrhage, descaling, discoloration, and redness in the body surface was used in the experiment. Additionally, the fish were randomly selected and screened for the possible presence of infectious microbes following the standard protocol (Johansen et al., 2006). The fish was acclimatized for 2 weeks in 200 L Fibre-reinforced plastic (FRP) tanks, supplied with proper commercial floating feed (CP: 30%, CL: 5%) @ 3–5% of body weight fed twice a day. The bacterial strain was cultured in 20 mL sterile Tryptone soya broth (TSB) in a 50 mL Erlenmeyer flask (Himedia, India) for 24 h at 28°C. The bacterial cells pellet was collected by centrifuging for 5 min at 5000 rpm and washed thrice with a sterile normal saline solution. Afterward, the pellets were resuspended in saline solution and the number of cells (CFU/ml) was estimated through the spread plate method. The experimental fish were intraperitoneally injected (20 numbers/ concentration) with 200 μL (1.2 × 10¹, 1.2 × 10², 1.2 × 10³, 1.2 × 10⁴, 1.2 × 10⁵, 1.2 × 10⁶ and 1.2 × 10⁷ CFU/mL) of the bacterial suspension. Control fish were injected with 200 μL of saline solution. Afterward, the fish were kept in an FRP tank and observed every 24 h until 120 h. To confirm Koch’s postulate, the bacteria were reisolated and identified from the liver, kidney, and blood of the moribund fish.

At first, the isolated strain was primarily distinguished as either Gram-negative or Gram-positive by the Gram-staining method. Later, the genomic bacterial DNA was isolated following the standard Sarkosyl method (Kumar et al., 2022). In brief, the quality of the DNA extracted was checked on agarose gel (1%) and quantified using Nano-drop (Eppendorf, Germany). The 16S rRNA gene was amplified by Gene Amp PCR system 9,700 thermal cycler (Applied Biosystems, Foster City, CA) using universal bacterial primers. 50 μL of the PCR reaction mixture consisted of 10× PCR buffer (5 μL), 50 mM MgCl₂ (1 μL), 10 mMdNTP (1 μL), 10 pmol of each primer (1 μL), isolated genomic DNA (100 ng), and Taq DNA polymerase (1 U). The PCR condition comprised of 2 min at 95°C initial denaturation, 35 cycles of 94°C for 30 s denaturation, 52°C for 60 s annealing and 72°C for 90 s extension and 7 min at 72°C of final extension. Subsequently, the PCR-amplified products were visualized on agarose gel (1.8%) (Behera et al., 2017). The amplified gene was sequenced in forward and reverse directions using ABI 373xl capillary sequencer (Applied Biosystem, Foster City, CA). Contig was prepared by aligning forward and reverse sequences using DNA baser 7.0.0. The amplified 16S rRNA gene from the strain was sequenced and submitted to GenBank, and accession number JF330411 was obtained. The BLAST analysis against the non-redundant database showed >99% identity with Aeromonas hydrophila (GenBank Accession Numbers: MT384379 and CP046954).

For RNA extraction, bacterial cells were harvested through centrifugation at 5000 rpm for 2 min, then washed with DEPC-treated water. The collected cells were frozen in liquid nitrogen, and we proceeded to isolate total RNA. Total RNA was isolated by a customised protocol combining the RaFlex total RNA isolation kit (Genei, India) and the Xcelgen RNA isolation kit (Xcelris Labs, India). The total RNA (1 μg) was estimated on a 1% denaturing agarose gel. The total RNA (4 ng/μl) was analysed on an Agilent Pico chip on an Agilent Bioanalyzer (Agilent, United States). The bioanalyzer profile of total RNA with a RIN value of 9 was obtained, which assured the presence of intact RNA in the sample. Total RNA of A. hydrophila isolated from bacterial cells was applied to enhance the RNA transcripts by collectively depleting ribosomal RNA molecules (rRNA) by applying MICROB Express bacterial mRNA purification kit (ThermoFisher Scientific, United States). Following rRNA-depleted mRNA isolation, the samples were arranged to form a transcriptome library employing the SOLiD Total RNA Seq kit (ThermoFisher Scientific, United States). The rRNA depleted the total RNA, i.e., the mRNA was split by applying the RNase III enzyme. Subsequently, the splitted mRNA was assorted with adapters following connection. The whole mRNA was converted to cDNA by reverse transcription. Size assortment was performed, and in-gel PCR was executed for library amplification. The extension was confirmed on 2% E-Gel size-select gel. Consequently, the transcriptome library was enumerated by a high-sensitivity DNA chip on a bioanalyzer. The sample was run in triplicate in support of the study of the transcriptome.

The execution of the sequencing run was performed by SOLID TOP SEQ-FRAG LIB KIT MM 50 (ThermoFisher Scientific, United States) for the construction of 40–50 million reads for each sample. Subsequent to the quality filtration (where the mean quality score is > = 20) and the adaptor trimming, the higher superiority A. hydrophila reads were amassed via the Velvet-Oases pipeline. A total of 9,273 transcript contigs were acquired. A total of 9,273 A. hydrophila transcript contigs were recorded on the sorted reads for the removal of misassembled transcript contigs. Following mapping, entire transcript contigs accepted the support conditions and were further employed for functional illustration.

Transcript contigs (9,273) of A. hydrophila were then used for functional annotation. The functional annotation was executed by lining up the transcript contigs with the bacterial non-redundant database of NCBI using the BLASTX program. A gene ontology study was performed to identify the significant genes in the A. hydrophila transcriptome analysis.

A transcriptome data study employed bioinformatics software to obtain complete information on the genes involved in A. hydrophila. BLAST (Altschul et al., 1990) was run on the databases of UniProt-TrEMBL as well as Swiss-Prot (Apweiler et al., 2004) in search of the broadly significant genes implicated. GO was initially characterised at different graph intensities on the basis of three important ontologies: biological processes, cellular components, and molecular functions. Gene Ontology interpretations were allocated employing the QuickGO¹ tool available at the European Bioinformatics Institute (Huntley et al., 2009).

Survival data of Labeo rohita were arcsine-transformed to satisfy normality and homoscedasticity requirements, as necessary. The data were then subjected to one-way analysis of variance followed by Duncan’s multiple range test using the Statistical Software Statistical Package for the Social Sciences version 24.0. p-values smaller at p < 0.05.

L. rohita exhibited hemorrhagic lesion and redness over the body surface which indicates the possible involvement of bacterial pathogen in mortality. Hence, the collected gut samples were processed for bacterial enumeration and characterization. The isolated bacterial strain displayed a Gram-negative character when subjected to Gram staining. Later, the bacterial strain was subjected to 16S rRNA gene and phylogenetic analysis. The amplified 16S rRNA gene from strain were sequenced and submitted to GenBank and the Accession number is OM900178. The BLAST analysis against non-redundant database showed 99.22% identity with Aeromonas hydrophila (GenBank Accession Number: MT384379 and CP046954). Afterwards, a phylogenetic tree was prepared with the sequence accessed from NCBI.

We next sought to investigate the virulence of isolated A. hydrophila strain possibly involved in the disease development and mortality of L. rohita. The control fish did not have mortality or clinical sign of disease during the experimental challenge. However, the fish challenged by intraperitoneal injection with bacterial suspension mostly developed subcutaneous haemorrhagic ulcers of about 0.7–1.8 cm diameter. Reddening at the injection sites and ulcer on the mouth region were also observed.

The A. hydrophila isolated from fish gut samples were used for RNA isolation, cDNA and library preparation. The schematic workflow of the whole transcriptome analysis of A. hydrophila was depicted (Figure 1). The total RNA integrity was analysed by applying the Bio-Analyzer, which supported the purity of the RNA for further processing (Figure 2). SOLiD sequencing was executed for the transcriptome library preparation of A. hydrophila.

To prepare the transcriptome of A. hydrophila, assembly was executed by applying the entire high-quality reads (after the elimination of replica reads in both samples). Following the quality filtration, i.e., considering the mean quality score > =20 as well as adaptor trimming, the high-quality reads of the A. hydrophila sample were accumulated through a de novo assembly pipeline employing Velvet v1.2.01 followed by Oases v0.2.04. BLASTX 2.2.26 and Blast2Go 2.5.0 software were used to further analyse the transcriptomic data (Conesa et al., 2005). A sum of 51,461,864 raw reads was produced with 2,573,093,200 numbers of bases with a mean read length of 50 bp, and 41,461,864 filtered reads were generated with 2,073,093,200 numbers of bases with a mean read length of 50 bp. A total of 9,273 transcript contigs were obtained. The assembly statistics and the allocation of transcript contigs conferring their length have been provided in Tables 1, 2. The average GC content of A. hydrophila transcript contigs was 52%. A. hydrophila has a high fraction of transcript contigs (33.58%) through GC content in an array of 50–55%. The allocation of transcript contigs depending on their length and the GC content distribution has been provided in Supplementary Figures S1, S2. BLASTX hits for every transcript contig reclined inside the cut-off value of 1e-06, and elevated significant resemblances were shown in the transcript contigs. Similarity distribution demonstrated that approximately 80% of transcript contigs had affirmative position lengths involving 50–95% of A. hydrophila transcript contigs. The E-value and similarity distribution of A. hydrophila have been provided in Supplementary Figures S3, S4.

The A. hydrophila 9,273 transcript contigs were recorded on the sorted reads to remove misassembled transcript contigs. A minimum of three plotted reads per transcript contig was required for corroboration of transcript contigs. After mapping, all the transcript contigs approved the confirmation standard and were further used for functional illustration. The A. hydrophila 9,273 transcript contigs were subsequently intended for functional annotation. The functional illustration was executed by lining up the transcript contigs with the NCBI non-redundant database of bacteria using the BLASTX program. BLASTX was affected in the illustration of 8,013 transcript contigs out of 9,273. As a result, 1,260 transcript contigs had notable negative BLAST hits. These interpreted transcript contigs were recorded on the GO database, owing to which 6,533 were GO annotated while 1,480 were not allocated any GO terms.

A species distribution study showed that the utmost proportion of A. hydrophila transcript contigs demonstrated momentous resemblance with Salmonella entrica followed by Klebsiella pneumoniae, Escherichia coli, Klebsiella sp., and Streptomyces sp. (Figure 3).

GO categorises different functions at different levels of graphs basically on the basis of three main processes, i.e., molecular function, cellular components, and biological processes. The significance of the enrichment is represented by an FDR adjusted value of p. The key illustrates the three areas that Gene Ontology covers, such as molecular activity, biological processes, and cellular components. Only annotations that have a significant false discovery rate (FDR) adjusted value of p of less than 0.05 are shown. In the case of A. hydrophila, the cellular component category shows maximum annotations at the 5th level, and the molecular function category shows maximum annotations next to the 7th GO level, while the biological process category shows maximum annotations next to the 6th GO level (Supplementary Figure S5). Gene ontology allocation might be used to better understand the distribution of annotated transcript contigs in precise ontology spheres such as molecular functions, cellular components, or biological processes. Since the total figure of annotated transcripts was 8,013 transcripts, they were mapped to GO. The maximum amount of information in transcript contigs was allocated towards molecular utilities. Individual transcript contigs might be characterised keenly over one GO domain. The transcript contigs allocated for molecular functions are 40, 38% for biological processes, and 22% for cellular components (Figure 4).

GO string allocations help stipulate every interpreted node encompassing GO efficient clusters. Transcripts coupled with comparable tasks are dispensed to identical GO functional clusters. Metabolic and cellular processes are the major biological processes in A. hydrophila, since approximately 76% of the transcript contigs in the sample are connected to these ontologies. The main metabolic process includes cellular metabolic processes, oxidation–reduction progression, primary metabolic processes, etc., whereas the main cellular progression includes cellular metabolic processes, cell division, and the cell cycle. Binding activity is the major molecular function within A. hydrophila since about 38% of the transcript contigs within the sample are interrelated to this ontology. The binding activity comprises protein binding, nucleotide binding, binding to ions, binding to nucleic acid, binding to cofactors, binding to lipids, binding to vitamins, and so forth. Approximately 50% of transcript contigs are accountable for the catalytic activities. This embraces transferase action, hydrolase action, oxidoreductase action, ligase action, lyase action, isomerase activity, and so forth. The cell is the most important cellular constituent in A. hydrophila and composes about 71% of the entire transcript contigs. This comprises the cell part. Around 16% of transcript contigs are related to macromolecular complexes, including protein complexes, protein-DNA complexes, and ribonucleoprotein complexes, while organelles constitute about 12% of transcript contigs. The GO level 2 sequence distribution of A. hydrophila is depicted in Figure 5.

GO conditions were collected into different stages for the three ontology spheres, i.e., the biological processes, the molecular functions, and the cellular components. Preliminary intensities denote the general function of the transcript contigs. When the intensity progresses, the transcript contigs’ role is explicitly added. All transcript contigs might be multi-functional and can recline within in excess of a single function. Therefore, the number of transcripts contigs intended for each stage varied significantly (Figure 6).

Several expressed genes have been identified following the transcript annotation and gene ontology analysis of A. hydrophila. Different virulent factors of A. hydrophila are encoded by several virulent genes. The genes thus identified are the potential virulent genes of Aeromonas sp., which contribute to their degree of virulence. The virulent genes include the outer membrane protein gene (ompW, ompA, ompTS, omp 38, omp 48), cytotoxin, amylase lipase, and haemolysin gene (Table 3), which have been reported earlier by Dash et al. (2014). A subset of these virulent genes might be a key player in the ability of the bacterium to cause disease. The products of such genes will facilitate the successful colonisation and survival of the bacterium in or cause damage to the host are considered as virulence or pathogenicity determinants. Hence, these genes can significantly interpret the pathogenic mechanisms exerted by A. hydrophila infection.