A. hydrophila isolates for disease challenge were recovered from diseased channel catfish from commercial aquaculture operations in western Alabama and the Delta region of west Mississippi. Briefly, catfish demonstrating symptoms typical of MAS were collected in a moribund state and submitted for diagnostic evaluation and necropsy at Auburn University in Auburn, AL or Mississippi State University's Aquatic Research and Diagnostic Laboratory at the Thad Cochran National Warmwater Aquaculture Center (NWAC) in Stoneville, MS. Liver and kidney tissues were sampled for aerobic bacterial cultures. Samples of tissue were homogenized in sterile phosphate buffered saline and portions of the homogenate streaked onto tryptic soy agar (TSA; Beckton Dickinson, Franklin Lakes, NJ) or brain heart infusion (BHI; Beckton Dickinson) for bacterial isolation. Pure cultures were identified as vAh strains by the vAh-specific qPCR method previously described and/or utilization of myo-inositol as a sole carbon source (Griffin et al., 2013; Hanson et al., 2014). From this sampling design, A. hydrophila isolates ML09-119, ML10-51K, S04-690, S14-296, and S14-452 were cryogenically preserved (mixed with 50% glycerol, stored at −80°C) and subsequently used in catfish immersion challenges.

Channel catfish were obtained as fry from the USDA-ARS Warmwater Aquaculture Research Unit housed at NWAC and reared to experimental size in 340 liter troughs supplied with 26 ± 2°C dechlorinated municipal water under pathogen-free conditions. All animal experiments were approved by, and conducted in compliance with, regulations of the Institutional Animal Care and Use Committee of the Aquatic Animal Health Research Unit (USDA-ARS) in Auburn, Alabama. Water temperature was maintained at 26 ± 2°C with a centralized heater. Prior to trials, heart, liver, head kidney, trunk kidney, spleen, brain, and skeletal muscle tissues were collected from 10 randomly sampled fish and sampled by culturing on BHI to verify fish were not presently infected with vAh.

Two hundred catfish fingerlings, with a mean weight of 111 ± 47 g and length of 19 ± 3 cm, were acclimated for 12 days in 56-L glass aquaria (10 fish per tank, 3 tanks per isolate, and 2 mock infected control tanks) containing about 50-L water prior to challenge. The immersion challenge was conducted using the recently described fin clip method (Zhang et al., 2016). At the time of infection, water volume was reduced to 15-L per tank. To sedate animals for handling, fish were netted from individual aquaria and placed into a container filled with 20-L of dechlorinated water containing 150 mg/L of buffered Tricaine-S (tricaine methanesulfonate; Western Chemical, Inc., Ferndale, WA). Once fish were anesthetized, the adipose fin was clipped at its base and fish were returned to respective aquaria for recovery from anesthesia.

For the bacterial challenge, 100 mL of tryptic soy broth (TSB) containing approximately 3.0 × 10⁹ CFU/mL of the respective A. hydrophila strains (ML-09-119, ML10-51K, S04-690, S14-296, S14-452, and ZC1) was added to each of three aquaria, resulting in an approximate challenge dose of ~2.0 × 10⁷ CFU/mL. Two tanks served as mock infected controls, receiving only 100 mL sterile TSB. After 1-h exposure, water flow to aquaria (0.5 L/min) was resumed. Fish mortality was monitored daily for 7 days. At least 50% of dead fish were sampled for confirmation for the presence of vAh in liver and kidney tissues using M9 minimal medium containing 0.3% (w/v) myo-inositol (M9I) agar (Hanson et al., 2014). Moribund fish were removed from aquaria daily and surviving fish were euthanized by at least 15 min exposure to 300 mg/L buffered Tricaine-S solution, then necropsied. Heart, liver, head kidney, trunk kidney, spleen, brain, and skeletal muscle tissues were harvested and fixed in 10% buffered formalin for histopathology. Tissues were also collected and used for bacterial identification and quantitation. Formalin fixed tissues were processed and embedded in paraffin. The tissues were cut in four micron sections, stained with hematoxylin and eosin and evaluated for microscopic lesions by light microscopy.

In total, 61 complete and draft Aeromonas spp. genomes were included in this study (Table 1) which included A. hydrophila representative genomes as well as Aeromonas spp. genomes that are now recognized as having had an erroneous affiliation with A. hydrophila (Figueras et al., 2014; Beaz-Hidalgo et al., 2015a). These genomes were retrieved from the US National Center for Biotechnology Information (NCBI) GenBank database and included two A. caviae isolates, one A. dhakensis isolate, one A. enteropelogenes isolate, one A. media isolate, one A. molluscorum isolate, one A. taiwanensis isolate, one Aeromonas sp. isolate, 26 A. hydrophila fish vAh disease isolates, and 28 non-vAh A. hydrophila isolates.

Strains representative of different vAh lineages were selected for Illumina sequencing based on results of vAh genotype-specific PCR (see below) from a total of 38 suspected vAh isolates recovered from diagnostic cases in MS between 2013 and 2015. Strains were selected so that isolates represented MAS outbreaks from multiple geographically discrete operations. Genome sequencing with 250 bp read-length using paired-end sequencing was performed on the Illumina MiSeq platform using the Nextera XT kit (Illumina, San Diego, CA) to prepare bar-coded fragment libraries according to the manufacturer's protocol. Sequence reads were trimmed and quality sequence reads were assembled de novo using the CLC Genomics Workbench (Qiagen, Redwood City, CA) using default settings. vAh strain draft genomes (n = 14) were generated for strains Ahy_Idx7_1, ALG15-098, IPRS-15-28, ML10-51K, S13-612, S13-700, S14-230, S14-296, S14-458, S14-606, S15-130, S15-242, S15-400, and S15-591 (Supplemental Table 1). In addition to a standard Illumina MiSeq run for vAh strain S14-452, a NxSeq 20 kb mate pair library was constructed and sequenced using an Illumina MiSeq at the Lucigen Corporation (Middleton, WI). The de novo assembly from the standard Illumina MiSeq sequences resulted in 13 contigs (80.15 average coverage) whereas the combination of these sequences together with the mate pair-derived sequences using de novo assembly with SPAdes (v3.5.0) resulted in a complete genome sequence.

Evaluation of the previously described vAh-specific qPCR primer set (listed as 2986L and 2986R in Table 2) (Griffin et al., 2013) was performed in silico using Geneious v. R9 with a maximum mismatch setting of two bases and a band prediction interval of between 100 and 1000 bases, which predicted that numerous vAh isolates (Ahy_Idx7_1, J-1, AL09-79, AL10-121, JBN2301, ML09-121, ML09-122, NJ-35, PB10-118, S13-612, and S15-591) would not produce an amplicon. To address this potential pitfall, we evaluated the primer set using touchdown PCR on all available and relevant vAh and non-vAh isolates (data not shown). Touchdown PCR was performed on an Eppendorf Mastercycler Gradient S with 50 ng of template gDNA extracted from each isolate (E.Z.N.A.® Bacterial DNA Kit; Omega Biotek, Georgia, USA), 13 μl of Econotaq Plus Green 2x MasterMix (Lucigen, Madison, Wisconsin, USA), and 20 picomoles of each primer in a 25-μl reaction. Cycling parameters comprised of an initial denaturation of 94°C for 3 min; 10 cycles of 94°C for 30 s, 68°C for 30 s (−1^oC per cycle), and 72°C for 1 min; followed by 25 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 1 min; and a final extension at 72°C for 5 min. Amplicons were resolved on a 1% agarose gel and stained with ethidium bromide.

In order to have a PCR assay that differentiates vAh from non-vAh strains while encompassing both US and Asian vAh isolates, new vAh genotype-specific primers (listed as vAh-SerF and vAh-SerR in Table 2) were designed using comparative analysis of annotations produced by the Rapid Annotations using Subsystems Technology (RAST) v2.0 server. These primers target a serine protease gene sequence and are predicted to produce a 502 bp amplicon (data not shown). Touchdown PCR conditions and validation methods were identical to those described previously.

In addition to primers that help to identify vAh strains, A. hydrophila lineage-specific primers were developed that were specific to vAh lineages represented by strains ML09-119 (NC_021290), S14-452 (SAMN05256776), and ZC1 (SAMN02404465) (Table 2). The ML09-119 lineage-specific primer set targets the tnsA endonuclease gene and produces a 246 bp amplicon, the S14-452 lineage-specific primer set targets the COG3339 genetic locus and produces a 350 bp amplicon, and the ZC1 lineage-specific primer set targets a hypothetical protein and produces a 400 bp amplicon.

To determine the vAh genotype affiliation of disease isolates in AL and MS, genomic DNA was isolated from each isolate (Gentra Puregene DNA isolation kit; Qiagen, Hilden, Germany) and used as a template in a 25-μl PCR that comprised of 13 μl of Econotaq Plus Green 2x MasterMix (Lucigen, Madison, Wisconsin, USA), 20 picomoles of each oligonucleotide primer and 50 ng of template gDNA. Samples were run on a C1000 Touch thermal cycler (BioRad, Hercules, California, USA) with an initial denaturation of 94°C for 3 min followed by 35 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 1 min, with a final extension at 72°C for 5 min and amplicon verification was performed as previously described.

A core genome was created using both coding and non-coding sequences of 61 genomes labeled as A. hydrophila within GenBank, some of which were revealed to be misclassified. Specifically, any contigs less than 10 Kbp in size were first filtered from draft genomes to increase computational efficiency and decrease false positives by limiting the mapping of smaller fragments to non-homologous and/or multiple regions. In general, smaller fragments contribute less to the production of core genomes as a function of being flanked by areas of heavily repetitive sequences. In the assembly stage, these regions stopped growing and remained short. While core genomes produced by this method are often smaller, the removal of smaller contigs did not compromise the ability to form a phylogenetically-informative core genome.

Filtered sequence data were then submitted as FASTA files to the multiple whole genome alignment tool Mugsy v1.2.3 (Angiuoli and Salzberg, 2011) under default parameters. The resulting alignment was subsequently processed with GBLOCK v0.91b (Castresana, 2000) in order to identify regions of high conservation across all isolates. Parameters for retention by GBLOCK are dictated by the input alignment and were: A minimum of 31 and 51 sequences for conserved and flanked positions, respectively, a maximum of 8 contiguous, but non-conserved positions, a minimal block length of 10, and one-half of the sequences were allowed to possess gapped positions within a block. From the final alignment, a maximum likelihood (ML) phylogeny for the 61 Aeromonas spp. isolates, including 54 isolates labeled in GenBank as A. hydrophila, was inferred using RAxML v8.2.8 (Stamatakis, 2014) under the General Time Reversible model of evolution with estimated proportions of invariable sites and rate variation among sites (i.e., GTR+I+G) and 1000 bootstrap replicates to determine branch supports. Trees were visualized using Archaeopteryx v.beta 0.9901.

Following generation of a consensus sequence, the National Microbial Pathogen Data Resource (NMPDR) Rapid Annotations using Subsystems Technology (RAST) v2.0 server was used in conjunction with the SEED v2.0 algorithm to annotate the core genome and generate metabolic models (Aziz et al., 2008; Overbeek et al., 2014). These predictive models were evaluated using both protein-protein Basic Local Alignment Search Tool (e.g., BLASTp and BLASTx) algorithms through GenBank as well as the Joint Genome Institute's Genomes OnLine Database (GOLD) v5.0.

To assess overall genetic similarity, the average nucleotide identity (ANI) comparison of 61 Aeromonas spp. genomes was evaluated using JSpecies (v1.2.1) and cross-validated with the Konstantinidis lab ANI calculator (Richter and Rosselló-Móra, 2009; Rodriguez and Konstantinidis, 2014). According to criteria used for the genus Aeromonas, ANI-values >96% indicated strains belong to the same species (Colston et al., 2014; Beaz-Hidalgo et al., 2015a).

To evaluate differences in gene content among the 61 genomes, data from the PAthogen Resource Integration Center (PATRIC) protein family sorter tool, RAST/SEED gene annotations, the Pathogen Host Interaction database (http://www.phi-base.org), and the Virulence Factors of Pathogenic Bacteria databases (http://www.mgc.ac.cn/VFs) were combined with copy number data for previously identified virulence genes (Rasmussen-Ivey et al., 2016). Results were evaluated by comparing predicted virulence-associated genes with closely and disparately related isolates to identify genes that were differentially present. Notably, genes shared between vAh and non-vAh isolates or those unique to an individual strain were removed from downstream analyses. These data were subsequently coupled with screening of virulence factors and vAh-associated genes using searches of the NCBI GenBank database with MegaBLAST and BLASTn algorithms (Wheeler et al., 2003). Thresholds for absence were specific to the respective gene and were restricted to mutations altering the predicted functional domains of proteins in which the protein sequence in question returned a functionally divergent protein. Once validated, these gene clusters and virulence factors were transformed for statistical analyses using Orange Data Mining software v.3.3.5 and/or R Studio v.0.99.896. After pre-processing in R Studio using packages MuMIn v.1.15.6, randomForest v.4.6-12, and k-means v.0.1.1, heat maps of resultant data were generated using Orange Data Mining software. To analyze subclade differences, vAh and non-vAh strains with known virulence properties (n = 25) were evaluated using k-means at 20 clusters (100% between sum of squares/total sum of squares). In addition to these analyses, the T346 Secretion System Hunter (version is not published) was used to identify T6SS-associated gene clusters using Glimmer v3.02, and HMMER3 v3.1b2 (Martínez-García et al., 2015).

The S04-690 genome was previously found to be the genetic intermediate (raw genetic distance) between US catfish vAh isolates (represented by strain ML09-119) and the Asian carp vAh isolate ZC1 (Hossain et al., 2014). A total of 38 vAh isolates recovered from MAS outbreaks in Mississippi from 2013 to 2015 and confirmed as vAh by phenotypic (myo-inositol usage) and/or genomic tests (qPCR) (Griffin et al., 2013; Hanson et al., 2014) were subjected to vAh lineage-specific PCR. Testing using lineage-specific primer sets based on representative isolates ML09-119, S14-452, and ZC1 showed that of the 38 MS isolates recovered from MAS outbreaks between 2013 and 2015, 20 isolates from a single farm were positive for the ML09-119 lineage-specific amplicon, while 18 isolates across five geographically discrete farms were positive for the S14-452 lineage-specific amplicon. No isolates produced amplicons with the ZC1 lineage-specific primer set (data not shown). Isolates that produced an amplicon with either the ML09-119 or S14-452 lineage-specific primer sets did not produce amplicons when assayed with other lineage-specific primer sets (data not shown).

Results from touchdown PCR assays demonstrate that the original 2968 primer set (Griffin et al., 2013) reliably distinguishes US vAh from non-vAh isolates (data not shown). However, the Asian carp isolates J-1 and NJ-35 were unable to be experimentally evaluated and in silico analysis predicts that Asian isolates J-1 and NJ-35 would not produce a PCR amplicon using the 2968 primer set. To address this possible limitation, a new vAh genotype-specific primer set (listed as vAh-SerF and vAh-SerR in Table 2) was designed to be inclusive of all Asian and US isolates, based on in silico analyses of vAh and non-vAh isolates. The vAh-SerF/R primer set was empirically confirmed to produce a 502 bp amplicon in all available US vAh isolates while returning a negative result for non-vAh isolates (data not shown).

Alignment of the complete and filtered draft genomes of the 61 Aeromonas spp. genomes via Mugsy produced a matrix of 19,817,762 positions. Following processing with GBLOCK, the core genome of these 61 strains contained 32,401 blocks and a consensus of 3,776,490 bp. This included 120,049 variable (>2 different nucleotides) and 79,507 phylogenetically informative sites (nucleotides that contributed to sorting phylogenetic groups), with percent G+C composition of 62.6%. Notably, these conserved regions collectively have a higher percent G+C content than the 61-isolate average of 61.1% (with a range of 60% for the genome of the type strain of A. molluscorum 848^T to 63.2% for the one of A. taiwanensis LMG24683^T; p-value = 0.003). Across the core analysis, within complete sites, the average transition to transversion ratio was 1.437 for all sequence pairs, with a minimum of 0 transitions and 1 transversion (A. hydrophila S04-690 and S13-700) and a maximum of 13 transitions and 2 transversions (A. hydrophila 173 and 14).

Conserved sequences from the core genome analysis were used to infer phylogenetic relationships among the 61 Aeromonas spp. isolates (Figure 1A). The core genome phylogeny indicates with 100% bootstrap support that vAh strains form a monophyletic group that is fundamentally distinct from other A. hydrophila (Figure 1A). In addition to supporting the monophyly of these ST251 isolates, this phylogenetic analysis also revealed support for distinct sub-clades among vAh isolates, with 2660 single nucleotide polymorphisms identified between the vAh sub-clades. For example, A. hydrophila ZC1, isolated from a grass carp in China (Figure 1B), is closely affiliated with vAh strains isolated from catfish in MS (S14-452, S14-458, S15-130, S15-400, and S15-591). A key genotype that differentiates the vAh sub-clades is the presence and genetic organization of T6SS components (Figure 2). Each of the ZC1-affiliated strains, as well as the other strains isolated from carp in China (i.e., J-1 and NJ-35), were found to contain at least 80% of the core proteins necessary for a complete T6SS, with ZC1 and S14-452 having two separate T6SS-associated gene clusters, whereas J-1 and NJ-35 each have a single complete T6SS cluster (Figure 2). In contrast, vAh isolates from catfish in AL as well as other MS isolates (i.e., S13-612, S13-700, S14-606, and S14-296) formed a distinct subclade that lack the majority of the core T6SS genes (Figures 1B, 2). Notably, while AL and MS vAh strains lacked the majority of T6SS components, they consistently retain other genes that are involved in T6SS in other bacteria such as a valine-glycine repeat protein G (VgrG), a T4 bacteriophage tail-like hole forming protein (Leiman et al., 2009); hemolysin coregulated protein (Hcp), a repetitive tubular protein that is similar to the phage major tail protein GpV (Pell et al., 2009); the chaperone protein ClpB, a chaperone and ATPase that interacts with Hcp to translocate effectors (Shrivastava and Mande, 2008); and VasH, a putative transcriptional regulator (Kitaoka et al., 2011) (Figure 2).

The average nucleotide identity values of the 61 Aeromonas spp. genomes was determined (Figure 3). High ANI-values (>99%) were found for all vAh-vAh pairwise comparisons, supporting the core genome phylogeny (Figures 1B, 3). In contrast, all vAh comparisons with non-vAh isolates possessed ANI-values less than 97%. In accordance with previous results, 14 strains appear to have a discrepancy between the species classification listed in GenBank and the species affiliation indicated by ANI-values (Table 1).

Challenge of channel catfish with vAh isolates from AL and MS resulted in ≥ 60% mortality (Figure 4). Within these US vAh isolates tested, there were no observed differences in virulence among ML09-119, ML10-51K, S04-690, S14-296, and S14-452, based on Duncan's multiple range test (p-value > 0.05). The carp isolate ZC1 was less virulent than AL and MS isolates with only 26.7% mortality observed (Figure 4). Most mortality (~96%) occurred within 48 h post challenge for all isolates, including ZC1. All dead fish (100%) sampled for confirmation were positive for the presence of vAh in liver tissue. Control fish yielded no mortality from the mock challenge.

Cutaneous lesions observed in fish infected with vAh by immersion challenge included extensive hyperemia over the pale ventrum of the fish, in tissues surrounding and within the mouth and at fin bases (Figure 5A). There was also extensive hyperemia around the eyes and exophthalmos in some fish. Cutting into the muscle of the lateral body wall revealed multifocal to coalescing foci of congestion/hemorrhage. Gill lesions were variable; gills were pale in some fish and reddened in others. Internally, there was widespread hyperemia of abdominal organs as well as petechial and ecchymotic hemorrhages scattered over mesenteric tissues. The spleen was moderately to severely swollen and dark red (Figure 5B). Head and trunk kidneys were moderately edematous and red and friable when harvested. The intestinal tract was mildly to moderately dilated and red (Figure 5B). The liver was mildly to moderately swollen with slightly rounded edges (Figure 5B). Glisson's capsule was peppered with variable numbers of petechial and ecchymotic hemorrhages. The atrial chamber of the heart was dilated and filled with blood.

Histopathologic lesions observed were strongly suggestive of a septic disease and included edema and necrosis in many internal organs. Splenic lesions included ellipsoidal necrosis and congestion/hemorrhage in the splenic red pulp (Figure 5C). In scattered necrotic ellipsoids variable numbers of short, rod-shaped bacteria were observed. In the liver there was often necrosis of the acinar pancreatic tissue surrounding hepatic vessels (Figure 5D). Scattered acinar cells were rounded and necrotic and extracellular zymogen granules could be observed in the necrotic exudate. Small numbers of inflammatory cells including macrophages, lymphocytes, and neutrophils were observable in scattered areas of pancreatic acinar necrosis. Multifocally within the hepatic parenchyma were variable sized foci of hepatocellular necrosis characterized by the presence of small aggregates of degenerating and necrotic hepatocytes. Hematopoietic cells in the renal interstitum and renal epithelial cells lining scattered renal tubules were undergoing degeneration and necrosis. In sections of head kidney, the tissue was edematous and congested/hemorrhagic. There was scattered degeneration and necrosis of red and white cell precursors. Intestinal lesions were minimal in fish infected in this fin clip model and were limited to mild congestion and hemorrhage of the vessels in the lamina propria and vessels of the muscularis and serosa. In some sections of heart there was mild necrosis of myofibers in the myocardium. Vessel of the brain were often moderately dilated and congested but the neuropil of the cerebrum, cerebellum, and brain stem was normal. The epithelium of the gills was normal but branchial capillaries were sometimes dilated and congested with erythrocytes.

In order to robustly define the vAh pathotype-specific loci in this study, previous results on established A. hydrophila virulence factors as well as genes unique to vAh strains were used in combination with a clustering approach and a random forest decision tree to identify gene products that may contribute to functional differences in virulence. This approach confirmed previously described gene clusters for L-fucose and O-antigen biosynthesis as well as myo-inositol catabolism. Additionally, we identified predicted virulence factors conserved within all vAh strains (Table 3) which includes virulence factors well-known to be important in A. hydrophila pathogenesis (Rasmussen-Ivey et al., 2016). Among predicted virulence factors conserved among vAh strains, there were many genes uniquely associated with vAh strains that were not present in other sequenced A. hydrophila strains, including a L-serine dehydratase, a N-acetylmannosamine kinase, a N-acetylneuraminate lyase, and a gene product predicted to be required for queuosine biosynthesis (Table 4). In a more comprehensive approach, all RAST/SEED predicted genes from the 41 confirmed A. hydrophila genomes (26 vAh and 15 non-vAh) included in this study were evaluated on the basis of linkage with the vAh pathotype using exhaustive iterations of random forest modeling, which resulted in 26 genes uniquely associated with vAh by either presence/absence or by differential copy number when compared to non-vAh (Figure 6).

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.