Genetic diversification of penaeid shrimp infectious myonecrosis virus between Indonesia and Brazil.

Infectious myonecrosis virus (IMNV) is a pathogen of penaeid shrimp, most notably the whiteleg shrimp Litopenaeus vannamei. First discovered in L. vannamei from Brazilian aquaculture farms in 2003, IMNV was additionally confirmed in L. vannamei from Indonesian farms in 2006 and has since been found in numerous provinces there. Only two complete sequences of IMNV strains have been reported to date, one strain from the Brazilian state of Piauí collected in 2003 and another from the Indonesian province of East Java collected in 2006. In this study, we determined the complete sequences of two additional Indonesian strains, one from Lampung province collected in 2011 and another from East Java province collected in 2012. We also determined partial sequences for six other strains to enhance phylogenetic comparisons, which have heretofore been limited by the small number of reported sequences, including only one for an Indonesian strain. The new results demonstrate clear genetic diversification of IMNV between Indonesia and Brazil, as well as within Indonesia. Analyses of conserved sequence motifs suggest a revised RNA pseudoknot prediction for ribosomal frameshifting.

Introduction

The whiteleg shrimp Litopenaeus vannamei is a major component of modern aquaculture. Though native to the east Pacific between Mexico and Peru, it is now farmed across much of the tropical and subtropical world. For example, L. vannamei was officially introduced on Brazilian farms in 1983 and on Indonesian farms in 2001 (Briggs et al., 2004). Since then in Indonesia, it has been farmed in at least 17 provinces. In 2012, Indonesia produced >400,000 metric tons of shrimp, with a value of >USD 2.5 billion, ∼65% of which was cultured L. vannamei (Ministry of Marine Affairs and Fisheries, Indonesia, 2013).

Like other species of cultured penaeid shrimp, L. vannamei is susceptible to numerous infectious diseases, which can severely impact aquaculture production. Viral diseases of importance in cultured shrimp include Taura syndrome, yellowhead disease, infectious hypodermal and haematopoietic necrosis, and white spot disease, respectively caused by viruses from families Dicistroviridae (plus-strand RNA), Roniviridae (minus-strand RNA), Parvoviridae (single-strand DNA), and Nimaviridae (double-strand DNA) (Flegel, 2012, Lightner et al., 2012, Seibert and Pinto, 2012, Walker and Winton, 2010). Added to this list in 2002–2003 was a newly recognized disease, infectious myonecrosis (IMN) (Lightner et al., 2004, Nunes et al., 2004), caused by a double-strand (ds) RNA virus, infectious myonecrosis virus (IMNV), tentatively assigned to family Totiviridae (Poulos et al., 2006).

IMN usually appears acutely on a farm, with elevated mortalities and gross signs among younger shrimp, but then progresses more chronically, with persistent mortalities accumulating to 40–70% (Lightner et al., 2004, Lightner, 2012). The gross signs include whitened areas in the abdominal segments and tail fan, both regions rich in skeletal muscle. Histopathologically, these lesions are characterized by coagulative muscle necrosis with edema during the acute phase, progressing to liquefactive necrosis with hemocytic infiltration and fibrosis during the recovery or chronic phase. In some cases, perinuclear basophilic inclusion bodies are seen in the muscle cells and others (Lightner et al., 2004), consistent with the viral etiology of this disease. Horizontal transmission of IMNV has been shown to occur by cannibalism and is thought likely to occur through water as well (Lightner, 2011, Lightner, 2012). Vertical transmission from broodstock to progeny is also suspected.

The IMNV genome is nonsegmented; spans 7561 base pairs (bp) of dsRNA; and comprises two long, partially overlapping open reading frames (ORFs) (Nibert, 2007, Poulos et al., 2006, Senapin et al., 2007). ORF1 encodes a polyprotein of 1605 amino acids (aa) thought to be processed co- or post-translationally into four separate polypeptides, designated 1, 2, 3, and 4 (Nibert, 2007, Tang et al., 2008). Polypeptide 1 (93 aa) is predicted to contain a classical dsRNA-binding domain (dsRBD) (Chang and Ramos, 2005, Nibert, 2007, Poulos et al., 2006), which seems likely to be involved in suppression of host defense responses; polypeptides 2 and 3 have unknown functions (but see below); and polypeptide 4 (901 aa) is the major capsid protein (MCP) (Poulos et al., 2006). The two ORF1 processing cleavages between polypeptides 1, 2, and 3 appear to be mediated by 2A-like translational stop/continue motifs (Firth and Brierley, 2012, Luke et al., 2008, Nibert, 2007) and are thus likely to be ribosome, rather than protease, mediated. The nature of the ORF1 processing cleavage between polypeptides 3 and 4 remains unknown. ORF2 is in the −1 frame relative to ORF1 and has the potential to encode a 912-aa product, designated polypeptide 5 (Nibert, 2007), which includes regions sharing strong sequence similarity with the RNA-dependent RNA polymerases (RdRps) of other dsRNA viruses (Poulos et al., 2006). Its translation is thought to require −1 ribosomal frameshifting in the region of ORF1/ORF2 overlap, yielding a polypeptide 4/5 fusion product (MCP/RdRp; 1734 aa) (Nibert, 2007). Indeed, a predicted RNA “slippery” sequence and pseudoknot characteristic of −1 ribosomal frameshifting (Giedroc and Cornish, 2009) are found in the region of ORF1/ORF2 overlap.

The three-dimensional structure of the IMNV virion has been determined by electron cryomicroscopy and icosahedral image reconstruction (Tang et al., 2008). Its capsid is ∼450 Å in diameter but notably contains lollipop-like fiber complexes extending by >80 Å beyond the capsid at the icosahedral fivefold axes. These fibers are probably functional in binding to cell-surface receptors and perhaps also in membrane penetration during cell entry (Nibert and Takagi, 2013, Tang et al., 2008). It is not yet known whether polypeptide 2 (284 aa), polypeptide 3 (327 aa), or both are responsible for forming these fibers (Tang et al., 2008). The IMNV capsid otherwise comprises 120 copies of the MCP arranged with T  = 1 symmetry (2 copies of MCP per asymmetric unit) (Tang et al., 2008), and presumably 1 or 2 of these 120 copies is a component of the MCP/RdRp fusion product. Each virion would thereby contain 1 or 2 RdRp domains projecting into its interior cavity, where the RdRp(s) could function to package, replicate, and transcribe the viral genome, as known for other dsRNA viruses (Castón et al., 1997, Zhang et al., 2003).

IMN was first recognized as a distinctive disease in Piauí, Brazil, in 2002 (Nunes et al., 2004), and has since spread to other states in the Brazilian Northeast Region (da Silva et al., 2010, Lightner, 2012, Mello et al., 2011). Notably, in 2006, IMN was also recognized in Indonesia (Senapin et al., 2007). The Indonesian outbreak was first noted on the island of Java, in East Java province, Situbondo regency, and by 2007 had spread east to the neighboring island provinces of Bali and West Nusa Tenggara. In the second half of 2008, IMN appeared to the west in the island province of Bangka Belitung and in Lampung province on the island of Sumatra. By 2009, the disease was also reported in Central Java, West Java, and North Sumatra provinces and across the Java Sea in West Kalimantan and South Sulawesi provinces. Although IMN has been rumored to occur elsewhere in Asia, this has been disputed (Senapin et al., 2011).

Only two complete sequences of IMNV strains have been reported to date, one strain from Piauí, Brazil, collected in 2003 (Poulos et al., 2006) and another from East Java, Indonesia, collected in 2006 (Senapin et al., 2007). The sequences reveal that these two strains share 99.6% nucleotide (nt) sequence identity, leading Senapin et al. (2007) and others (Lightner, 2012, Walker and Winton, 2010) to suggest that the IMNV appearing in Indonesia might have originated in Brazil and possibly to have spread by unregulated transfer of shrimp broodstock or postlarvae for aquaculture production, as thought to have occurred for some other shrimp-disease agents (Flegel, 2006). Partial sequences of three other strains, all from Brazil in 2004–2007, are also available in GenBank (see Section 2 for all accession numbers), but this small number of reported sequences has limited phylogenetic comparisons. To better elucidate the possible source and ongoing evolution of IMNV in Indonesia, we newly determined sequences of several IMNV strains sampled in 2011–2012. The results demonstrate clear genetic diversification of IMNV between Indonesia and Brazil, and conserved sequence motifs relating to polyprotein processing, ribosomal frameshifting, and protein functions.

Materials and methods

IMNV-infected L. vannamei samples

Symptomatic shrimp were collected in early 2012 from multiple aquaculture farms in East Java, Lampung, and Bangka Belitung provinces, Indonesia, and fixed in alcohol for transport and storage. A frozen 2011 sample from Lampung province was also collected at this time. Aside from province and year, the samples were blinded with regard to the specific farms from which they came, and no records identifying the farms were maintained. An additional sample from early 2011 was obtained from the state of Ceará in northeast Brazil.

Complete and partial sequencing of IMNV strains

Total RNA was extracted from IMNV-infected L. vannamei samples and purified by RNeasy mini kit (QIAGEN). RT-PCR (QIAGEN OneStep kit) was then performed using primer pairs described below. The resulting, overlapping RT-PCR products were purified by QIAquick kit (QIAGEN). RNA extraction, RT-PCR, and PCR purification were performed according to the manufacturer's protocols.

For complete (whole-genome) sequencing of strains ID-LP-11 and ID-EJ-12-1, primer pairs 1–15 as specified by Senapin et al. (2007) were used. The purified, overlapping RT-PCR products were directly submitted for Sanger DNA sequencing of both DNA strands on Applied Biosystems 3730xl DNA analyzers at the Dana-Farber/Harvard Cancer Center DNA Resource Core, using the same primers as for amplification. The residue at each nt position in the IMNV genome was read de novo from both DNA strands, except for short regions (46–68 nt) near the genome termini, including the regions fixed in sequence by the two terminal primers, which were read only in the outward direction.

For partial sequencing of strains ID-LP-12-1, ID-LP-12-2, ID-BB-12, ID-EJ-12-2, and ID-EJ-12-3, we used primer pairs 5–10 (covering nt positions 1921–4920) specified by Senapin et al. (2007). Because pair 6 failed to yield the expected product for some samples, the following pair was added: forward 5′-CTATTAGGGCCGACGCTGT and reverse 5′-ATTTTGACGTGCTGCCATC. Pair 10 also failed to yield the expected product for ID-EJ-12-3 but was not replaced. The purified, overlapping RT-PCR products were then separately ligated into the pGEM-T Easy Vector plasmid (Promega), using T4 DNA ligase and an insert:vector molar ratio of 3:1, followed by transformation into Escherichia coli DH5α. Insert-containing colonies were picked from each plate and screened for correctly sized inserts by PCR with M13 primers and JumpStart RedTaq polymerase (Sigma–Aldrich). After purification again by QIAquick, the PCR products were submitted for Sanger DNA sequencing of both DNA strands on Applied Biosystems 3730 DNA analyzers at the University of Arizona Genetics Core, using M13 primers. Each residue between nt positions 1921–4920 of the IMNV genome (1921–4420 for ID-EJ-12-3) was read de novo from both DNA strands of 1–3 plasmid clones. For partial sequencing of BZ-11-UAZ219, the following primer was used for both RT-PCR and sequencing: forward 5′-AACACAAAATCTGCCAGCAA and reverse 5′-CCCAACCACCCAAATTCATA. Sequencing was performed directly from the purified RT-PCR products at the University of Arizona Genetics Core.

GenBank accession numbers

Accession numbers and abbreviations for IMNV strains: BZ-03 (AY570982.2), ID-EJ-06 (EF061744), BZ-04-ZS2011001 (KC200075), BZ-07-1 (HM030799), BZ-07-2 (HM357803), plus new entries ID-LP-11 (KJ636782), ID-EJ-12-1 (KJ636783), ID-EJ-12-2 (KJ636784), ID-LP-12-1 (KJ636785), ID-EJ-12-3 (KJ636786), ID-LP-12-2 (KJ636787), ID-BB-12 (KJ636788), and BZ-11-UAZ219 (KJ636789).

Sequence and phylogenetic analyses

Pairwise sequence alignments were performed using EMBOSS 6.3.1:needleall as implemented at http://mobyle.pasteur.fr/. Multiple sequence alignment was performed using MAFFT 7.146b as implemented with defaults at http://mafft.cbrc.jp/alignment/server/. For analyses of the 656-nt and 2520-nt regions as described in Section 3 and Table S1, the multiple sequence alignment was manually trimmed to those limits. For the 656-nt region, neighbor-joining phylogenetic comparisons were performed as implemented in MEGA6 (Tamura et al., 2013). Distances were computed using the maximum composite likelihood method, and branch support was determined by bootstrapping (1000 replicates).

Supplementary table related to this article can be found, in the online version, at http://dx.doi.org/10.1016/j.virusres.2014.05.013.

Table S1

Thirteen available IMNV sequences compared in this study.

Prior to maximum-likelihood and Bayesian-inference analyses, testing of nt-substitution models was performed using jModelTest 2.1 (Darriba et al., 2012). The best models were identified according to the Bayesian information criterion (BIC): 656-nt region, model TPM1uf; 2520-nt region, model HKY+I; and 7561-nt region, model HKY+I. For the 656-nt region, model HKY scored fourth (after TPM1uf, TPM1uf+I, and TPM1uf+G) and was used in place of model TPM1uf for consistency with the other two regions.

Maximum-likelihood phylogenetic analyses were performed using PhyML 3.0 (Guindon et al., 2010) as implemented at http://www.hiv.lanl.gov/content/sequence/PHYML/interface.html. The nt-substitution model and associated parameter values as needed (I, proportion of invariable sites; G, gamma shape parameter (4 substitution rate categories); and F, empirical nt residue frequencies) were used as indicated by jModelTest 2.1 as described above. In each case, the starting tree was optimized according to both topology and branch length, and tree improvement was performed according to the best of NNI (nearest neighbor interchange) and SPR (subtree pruning and regrafting). Branch support was determined by the aLRT (approximate likelihood ratio test) with SH-like supports.

Bayesian-inference phylogenetic analyses were performed using MrBayes 3.2 (Ronquist and Huelsenbeck, 2003) as implemented at http://www.phylo.org/sub_sections/portal/. For the 656-nt region (model HKY), the following customized commands were used: lset nucmodel = 4by4 nst = 2; prset statefreqpr = fixed (empirical); mcmc ngen = 2,000,000 printfreq = 20,000 samplefreq = 20,000 nruns = 2 nchains = 8; sump burnin = 1,000,000 nruns = 2; sumt burnin = 1,000,000 nruns = 2. For the 2520-nt and 7561-nt regions (model HKY+I), the same customized commands were used except that the following were added: lset rates = propinvar; prset pinvarpr = fixed (0.93). In each case, running the analysis for ngen = 2,000,000 ensured that the average standard deviation of split frequencies fell below 0.03 and persisted there. Branch support was determined by posterior probability testing. Final trees obtained with each of the methods were refined for publication in FigTree 1.4.0.

Frameshifting/pseudoknot predictions were performed using KnotInFrame and pKiss as both implemented at http://bibiserv2.cebitec.uni-bielefeld.de/rna. Default settings were applied in each case, except for using the thermodynamic model parameters of Andronescu et al. (2007) and a temperature of 29 °C. Homology detection and structure prediction by HMM–HMM comparison for IMNV polypeptide 1/dsRBD was performed using HHpred as implemented with defaults at http://toolkit.tuebingen.mpg.de/hhpred.

Results

Complete sequencing of Indonesian IMNV strains

L. vannamei samples were collected from aquaculture farms in East Java, Bangka Belitung, and Lampung provinces in early 2012, including one sample that had been stored frozen since early 2011. This Lampung 2011 sample (ID-LP-11) and one East Java 2012 sample (ID-EJ-12-1) were chosen for complete (whole-genome) sequencing. After RNA extraction from the shrimp tissue, and RT-PCR to amplify overlapping regions of the IMNV genome, sequencing was performed directly from both strands of the purified RT-PCR products.

The genomes of both ID-LP-11 and ID-EJ-12-1 were found to span 7561 bp and to encompass two long, partially overlapping ORFs on one of the RNA strands, i.e., the plus strand. The upstream ORF, ORF1, covers 4818 nt (encoding 1605 aa) from the first in-frame AUG codon at nt positions 136–138 to the first in-frame stop codon at nt positions 4951–4953, leaving a 5′ nontranslated region of 135 nt. The downstream ORF, ORF2, covers 2739 nt (encoding 912 aa) from the first in-frame codon at nt positions 4752–4754 to the first in-frame stop codon at nt positions 7488–7490, leaving a 3′ nontranslated region of 71 nt not including the ORF2 stop. ORF1 and ORF2 thus overlap by 199 nt (positions 4752–4950) not including the ORF1 stop, and ORF2 is in the −1 frame relative to ORF1. These results agree with those reported by Senapin et al. (2007) for the IMNV strain from East Java collected in 2006 (herein referred to as ID-EJ-06).

The genome sequence of the first IMNV isolate, from Brazil in 2003 (Poulos et al. (2006) (herein referred to as BZ-03), was originally reported to be 7560-bp long, differing from the genome sequences of Indonesian strains ID-EJ-06, ID-LP-11, and ID-EJ-12-1 as noted above by a single nt deletion near the 3′ end of ORF2, slightly shortening ORF2 and the predicted length of its encoded polypeptide. Given this discrepancy, we asked the previous authors to recheck their original BZ-03 sequence, and they indeed found an error, indicating that the genome length of BZ-03 is also 7561 bp, with ORF2 encompassing the same nt positions (7488–7490) as in the Indonesian strains. This correction has been made to the GenBank entry for BZ-03.

Pairwise comparisons revealed that the complete sequences of the one Brazilian and three Indonesian IMNV strains aligned without gaps and shared nt-sequence identities of 98.7–99.6%, or only 30–102 nt substitutions per pair (Table 1 ). Pairwise scores for aa-sequence identities among these strains showed similar relative values (98.9–99.5%, or only 12–28 aa substitutions per pair). The two strains exhibiting the lowest pairwise-identity scores were ID-LP-11 and ID-EJ-12-1, suggesting they may represent diverging lineages of IMNV in Indonesia. Based on these findings and the phylogenetic results described below, a crude estimate of the nt-substitution rate for IMNV is ∼8 substitutions/genome/year between 2003 and 2012, or ∼1 × 10−3 substitutions/site/year, which approximates that for avian reoviruses suggested by Liu et al. (2003) and is at the higher end of that suggested for dsRNA viruses by Stenger et al. (2010).

Table 1

Pairwise comparisons of the four IMNV whole-genome sequences.

	BZ-03	ID-EJ-06	ID-LP-11	ID-EJ-12-1
BZ-03	100a	99.6 (30)	99.1 (66)	99.0 (76)
ID-EJ-06	99.5 (12)	100	99.3 (56)	99.3 (56)
ID-LP-11	99.2 (19)	99.4 (14)	100	98.7 (102)
ID-EJ-12-1	99.0 (25)	99.2 (19)	98.9 (28)	100

Scores (%): nt above diagonal, aa (ORF1/2 fusion) below diagonal. Parenthetical values are raw numbers of substitutions within each pair.

Partial sequencing of additional IMNV strains to enhance phylogenetic comparisons

To enhance phylogenetic and other comparisons, we performed partial sequence determinations for five more of the IMNV-infected samples collected in the East Java, Lampung, and Bangka Belitung provinces of Indonesia in 2012 (ID-EJ-12-2, ID-EJ-12-3, ID-LP-12-1, ID-LP-12-2, and ID-BB-12) and also for an IMNV-infected L. vannamei sample obtained from the northeast Brazilian state of Ceará in 2011 (BZ-11-UAZ219) (see Section 2). The three other GenBank-deposited IMNV sequences, partial sequences for Brazilian strains collected in 2004–2007 (BZ-04-ZS2011001, BZ-07-1, and BZ-07-2), were included in these analyses. Alignment of the complete and partial sequences of these 13 total IMNV strains showed them to be free of gaps, except for the regions missing from the partial sequences.

Due to the limited divergence among available IMNV sequences reflected in Table 1, we used nt, rather than aa, sequences for phylogenetic comparisons (Townsend et al., 2008). In addition, given the availability of both complete and partial IMNV sequences and therefore the possibility that the data missing from the partial sequences (Table S1) might reduce the reliability of the resulting trees (Shavit Grievink et al., 2013, Wiens and Morrill, 2011), we analyzed three different, progressively larger portions of the aligned sequences. First, a 656-nt region spanning positions 3124–3779 was fully represented in 12 sequences, but not at all in one sequence (BZ-07-02), and was analyzed as the largest available set of IMNV sequences with no missing data among those 12 strains. Second, a 2520-nt region spanning positions 1921–4440 was fully represented in nine sequences and partially represented in the four others (BZ-04-ZS2011001, BZ-07-1, BZ-07-2, and BZ-11-UAZ219), and was analyzed as a version of all 13 sequences with limited missing data among them. Third, the full lengths of the available sequences were analyzed as a version of all 13 sequences with a greater amount of missing data.

Trees for the 656-nt region reconstructed using three different phylogenetic approaches—neighbor joining (distance) (Fig. 1A), maximum likelihood (Fig. 1B), and Bayesian inference (Fig. 1C)—were notably similar. One main feature common to the three trees was that the Brazilian and Indonesian strains constituted separate subclades of IMNV, i.e., were not intermixed. In addition, the original Indonesian strain, ID-EJ-06, consistently branched near the base of the Indonesian subclade, consistent with its earliest sampling among the Indonesian isolates. Lastly, the two recent Indonesian strains for which complete sequences were determined, ID-LP-11 and ID-EJ-12-1, consistently localized on different branches within the Indonesian subclade, in line with their low pairwise-identity scores in Table 1. Differences among the three trees concerned specific branch positions for a few strains within the Indonesian subclade.

Fig. 1

Phylogenetic comparisons for 656-nt region of IMNV genome. Trees were reconstructed according to different phylogenetic approaches as described in Sections 2 and 3: neighbor joining (distance) (A), maximum likelihood (B), and Bayesian inference (C). All trees are unrooted and displayed as radial phylograms. Branch-support values (%) were obtained as described in Section 2, and branches with <70% support were collapsed to a polytomy at the preceding node. IMNV strains for which sequences have been previously reported are labeled in gray, those newly reported here are in black, and those for which whole-genome sequences are available are underlined.

For the 2520-nt and full-length (7561-nt) regions, which were characterized by different extents of missing data for the nine partial sequences (Table S1), neighbor-joining trees were not reconstructed due to the lack of overlap between the partial sequences for BZ-07-01 and BZ-07-02. Nevertheless, both maximum-likelihood trees (Fig. 2A and C) and Bayesian-inference trees (Fig. 2B and D) reconstructed for these larger regions of the IMNV genome were notably similar to the trees for the 656-nt region (Fig. 1), including the same three main features described in the last paragraph. Differences among the trees were again noted with regard to specific branch positions for a few strains within the Indonesian subclade. Thus, in this case, data missing from the partial sequences appeared to have little or no impact on the primary phylogenetic conclusions.

Fig. 2

Phylogenetic comparisons for 2520-nt and full-length (7561-nt) regions of IMNV genome. Maximum-likelihood trees (A, C) and Bayesian-inference trees (B, D) were reconstructed for the 2520-nt region (A, B) and the 7561-nt region (C, D) as described in Sections 3, 2. All trees are unrooted and displayed as radial phylograms. Branch-support values (%) were obtained as described in Section 2, and branches with <70% support were collapsed to a polytomy at the preceding node. IMNV strains are labeled as in Fig. 1.

Analysis of conserved sequence motifs among IMNV strains

Given the overall high degree of sequence conservation observed among the IMNV strains, it might have been expected that nt- and aa-sequence motifs previously identified in BZ-03 and ID-EJ-06 (Nibert, 2007, Poulos et al., 2006, Senapin et al., 2007) would be highly conserved. Regardless, availability of the new sequences determined for this report provided a useful opportunity to reexamine these motifs.

At the RNA level, previously identified nt-sequence motifs in IMNV pertaining to the −1 ribosomal frameshifting in the overlap region between ORF1 and ORF2 (Nibert, 2007) were encompassed by eight of the available sequences, including six first reported here (Fig. 3A). The proposed “slippery” sequence (GGGUUUU; nt positions 4803–4809), within which frameshifting is thought to occur, was found to be fully conserved in these eight sequences. The predicted pseudoknot sequence located immediately downstream (nt positions 4818–4862) was also conserved, except for a unique G → A substitution at position 4828 in all three strains from Lampung province. Notably, this nt substitution altered the pseudoknot prediction by introducing a disfavored A:C base pair into one of the previously predicted main stems (Nibert, 2007). Further considerations about this finding are discussed below.

Fig. 3

Sequence motifs. (A) Ribosomal −1 frameshifting signals. Sequences for nt positions 4803–4869 are shown. The proposed slippery sequence in the region of ORF1–ORF2 overlap is indicated by double underlines. The two main stems (S1A:S1B and S2A:S2B) of the predicted H-type pseudoknots are indicated by heavy underlines. Additional stems in loop 1 (s1a:s1b) or loop 2 (s2a:s2b) are indicated by regular underlines. MFE values calculated at 29 °C are shown at right. Gray shading highlights the unique G → A substitution in all three Lampung strains. At top is the original prediction for BZ-03 from Nibert (2007). In the middle are the predicted H-type pseudoknots from KnotInFrame for the eight IMNV strains whose determined sequences overlap this region. At bottom is a result from pKiss, shown only for BZ-03, to illustrate that the pseudoknot predicted for the Lampung strains can also be formed by the other strains, albeit with a lower predicted stability. (B) Polyprotein processing motifs as they occur between the intervening ORF1 polypeptides (pep1/dsRBD, pep2, pep3, and pep/MCP). Arrows indicate the “cleavage” sites between the polypeptides. The 2A-like elements (2ALE) at the pep1-2 junction (#1; aa positions 86–94) and the pep2-3 junction (#2; aa positions 370–378) are labeled; aa positions 695–714 are shown across the pep3-4 junction. Gray shading highlights the unique Val → Ile substitution in 2ALE#2 of ID-EJ-12-1. The partial sequences determined for ID-EJ-12-2, ID-EJ-12-3, ID-BB-12, ID-LP-12-1, and ID-LP-12-2 do not overlap the 2ALEs, as indicated by parentheses. (C) The aa sequence of polypeptide 1/dsRBD from IMNV strain ID-LP-11 is aligned with the 69-aa, second dsRBD from Xenopus RNA-binding protein A, for which a co-crystal structure with dsRNA (PDB ID, 1DI2) has been reported (Ryter and Schultz, 1998). The alignment and secondary-structure predictions (sspred; H/h, α helix; E/e, β strand; C/c, coil) between query (Q) and target (T) are those output by HHpred. The central line of text shows identities (|) and similarities (:,.) between the aligned aa sequences. Gray shading highlights the unique Lys → Glu substitution in ID-LP-11 relative to other IMNV strains. Heavy lines indicate the two regions of α helix and lighter lines indicate the three regions of β strand comprising the αβββα core seen in 1DI2. Residues that provide R-group contacts to dsRNA in the 1DI2 structure are shown by asterisks, and labeled with the dsRNA groove in which they bind.

At the protein level, one set of previously identified aa-sequence motifs in IMNV pertain to co- or post-translational processing of the 1605-aa, ORF1-encoded polyprotein (Nibert, 2007, Poulos et al., 2006). The two 2A-like, ribosomal stop/continue (Firth and Brierley, 2012) elements were encompassed by only the four complete sequences (Fig. 3B). The 2A-like element between polypeptides 1/dsRBD and 2 (aa positions 86–94 in the ORF1 polyprotein) was fully conserved in the four sequences (GDVESNPGP). On the other hand, the 2A-like element between polypeptides 2 and 3 (aa positions 370–378 in the ORF1 polyprotein) was conserved in three of the sequences (GDVEENPGP) but included a Val → Ile substitution in ID-EJ-12-1 (GDIEENPGP). Since Ile has been found in this position in 2A-like elements from several other RNA viruses (porcine teschoviruses (family Picornaviridae), turkey avisiviruses (family Picornaviridae), and Israeli acute paralysis viruses (family Dicistroviridae) (e.g., GenBank accession nos. AF296098, KC614703, and EU436455)), this change in ID-EJ-12-1 seems unlikely to be of major importance, though it might affect the relative efficiency of “cleavage” by this element. The ORF1 processing cleavage between polypeptides 3 and 4/MCP involves an unknown mechanism, but was fully conserved within the nine available IMNV sequences that encompass this region, including seven first reported here (Fig. 3B).

Polypeptide 1/dsRBD (aa 1–93) was encompassed by only the four complete IMNV sequences. It was found to be fully conserved in three of them, but in ID-LP-11 had a Lys → Glu substitution at aa position 12. Comparisons with other dsRBD proteins by using the programs HHpred and HHblits (Remmert et al., 2011) suggested this residue lies within the first β strand of the conserved αβββα structure and is not likely to be directly involved in dsRNA binding (Fig. 3C) (Chang and Ramos, 2005). This substitution may therefore have little effect on the function of this polypeptide, though that remains to be determined by direct testing. All aa substitutions found in ORF1 polypeptides 1/dsRBD, 2, 3, and 4/MCP in the four complete sequences are summarized in Table S2. Substitutions at two aa positions in polypeptide 4/MCP were the only ones that differentiated between multiple (three or four) Brazilian strains and multiple (eight) Indonesian strains (Table S2), and might therefore serve as markers for the two, geographically segregated IMNV subclades.

Supplementary table related to this article can be found, in the online version, at http://dx.doi.org/10.1016/j.virusres.2014.05.013.

Table S2

Substitutions at the aa level in the four IMNV whole-genome sequences, relative to original strain BZ-03.

The ORF2 product polypeptide 5/RdRp, which is expressed only in fusion with polypeptide 4/MCP after −1 ribosomal frameshifting, was also encompassed by only the four complete IMNV sequences. Poulos et al. (2006) noted that a central portion of this polypeptide contains the eight aa-sequence motifs characteristic of other dsRNA-virus RdRps (Routhier and Bruenn, 1998). All eight of these RdRp motifs were fully conserved in the four whole-genome sequences. Throughout polypeptide 5/RdRp (893 aa following the predicted ribosomal frameshift), a total of only nine aa substitutions were found among the four sequences, as also listed in Table S2. Four of these substitutions were found in the N-terminal region of polypeptide 5/RdRp preceding the characteristic RdRp motifs; two were found in the RdRp-motif region (between the eight motifs themselves); and three were found in the C-terminal region of polypeptide 5/RdRp following the RdRp motifs.

Discussion

Source of IMNV in Indonesia

Several authors have suggested that IMNV likely arrived in Indonesia from Brazil, through unregulated transfer of L. vannamei broodstocks or postlarvae for use in aquaculture production (Lightner, 2012, Senapin et al., 2007, Walker and Winton, 2010). The subclade relationships shown in Fig. 1, Fig. 2 are consistent with this suggestion, but are also consistent with the opposite direction of transfer, since the direction of time is undefined in unrooted trees such as these. The earlier appearance of IMN in Brazil (Lightner et al., 2004, Nunes et al., 2004, Senapin et al., 2007) is certainly consistent with the possibility that Indonesian IMNV may have originated in Brazil. Moreover, since L. vannamei has been cultured in Brazil since 1983 but in Indonesia since only 2001 (Briggs et al., 2004), transfer from Indonesia to Brazil seems less likely. Another possibility consistent with the current phylogenetic results, however, is that related ancestors of the Brazilian and Indonesian IMNV subclades may have been separately introduced to the two countries from a third, still unknown, geographical source.

A fundamental finding of the current study is the clear genetic diversification of IMNV between Indonesia and Brazil, reflected by their separation in distinct subclades. This finding does not rule out a transfer from Brazil as the original source of IMNV in Indonesia, but does appear to indicate that such transfers have not continued to occur, probably reflecting greater efforts to restrict international transfer of IMNV-infected aquaculture stocks once the new disease was recognized and more widely appreciated as an important problem. Despite the presence of these two subclades, it is important to recognize that they should continue to be considered members of the same viral species, since their extent of divergence remains so low to date (Table 1).

Presence of IMNV in a wild host has not been formally reported to date. Furthermore, since IMNV can experimentally infect at least two other species of cultured shrimp, L. stylirostris and Penaeus monodon (Tang et al., 2005), the possible role of such other species, including other crustaceans (e.g., crabs, lobsters, crayfish, or copepods), in the origin and spread of IMNV deserves consideration. The availability of IMNV-related sequences obtained from other host species, for example, would allow investigations relating to virus–host coevolution (Göker et al., 2011) that may shed light on the specific origin of IMNV in Brazilian and Indonesian L. vannamei.

Genetic diversification of IMNV within Indonesia

Given that IMNV-associated disease seems to have been first recognized in Indonesia in East Java province (Senapin et al., 2007), it is reasonable to conclude that the Indonesian outbreak initiated there. That conclusion remains less than certain, however, in that the outbreak might have initiated in another province, but not been recognized or reported there. The phylogenetic trees in Fig. 1, Fig. 2 are nevertheless consistent with an East Javan epicenter in that strain ID-EJ-06 consistently constitutes the earliest, or one of the earliest, branches in the Indonesian subclade.

The different trees in Fig. 1, Fig. 2 exhibit some discrepancies with regard to specific branch positions for a few strains within the Indonesian subclade. These discrepancies make it difficult to arrive at firm conclusions about the trajectory of IMNV diversification in Indonesia, but it seems likely that there are at least two distinguishable lineages within this subclade evolving in parallel, represented by the two Indonesian strains for which complete sequences were determined in this study, ID-LP-11 and ID-EJ-12-1. Even though one consistent feature of most of the trees is a branch comprising the three recent East Java strains, the sampling of strains remains limited to date, and any ongoing human exchange of L. vannamei broodstocks or postlarvae between different regions of Indonesia might be expected to obscure or complicate such regional patterns over longer periods or with more intensive sampling. Thus, although strains representing distinct lineages of Indonesian IMNV as suggested in Fig. 1, Fig. 2 might continue to co-circulate for years to come, descendants of one lineage might come to dominate for particular reasons. On the other hand, different environments or practices between aquaculture farms in different regions might provide selective pressures that would reinforce the regionally based genetic diversification of IMNV that we found in this study.

There is also the possibility of natural exchanges of IMNV between nearby regions via water or wildlife, which may contribute to regional mixing. For example, because of its nonenveloped virions (Poulos et al., 2006, Tang et al., 2008), IMNV, similarly to Taura syndrome virus, is likely to remain infectious in the gut and feces of seabirds that feed on dead or moribund shrimp at farms with ongoing IMN epizootics, and to be spread within and among farms by feces or regurgitated carcasses (Vanpatten et al., 2004). In contrast, exchange of IMNV between Indonesia and Brazil, or vice versa, almost certainly cannot occur by natural means, given the locations of these countries on opposite sides of the globe (though interestingly at similar latitudes just south of the Equator), emphasizing the likely role of humans in the dispersal of IMNV in the early 21st century.

Ribosomal frameshifting signals

The original pseudoknot prediction for BZ-03 (Nibert, 2007) was an H-type (two-stem) pseudoknot, but with an additional stem within loop 2 (labeled s2a:s2b in Fig. 3A), as also found for SARS coronavirus pp1a-1b (Giedroc and Cornish, 2009). The distance between the end of the frameshifting “slippery” sequence and the start of this predicted pseudoknot is somewhat long, 8 nt, and the newly calculated minimum free energy (MFE) of this pseudoknot at 29 °C (average temperature recorded during IMNV epizootics) is −12.0 kcal/mol, vs. −10.1 kcal/mol at 37 °C. According to the program KnotInFrame (Theis et al., 2008) and others, however, this particular H-type pseudoknot is not the preferred one; instead, in the newly predicted, preferred pseudoknot, the S1 stem is shifted such that the distance between the end of the slippery sequence and the start of the pseudoknot is only 5 nt, and the pseudoknot has a greater predicted stability, MFE −12.9 kcal/mol at 29 °C (Fig. 3A). Moreover, the G → A substitution in the three Lampung strains does not strictly eliminate this pseudoknot, but rather causes the S2 stem to be shifted and allows an additional stem to form within loop 1 (labeled s1a:s1b in Fig. 3A), yielding a related, preferred pseudoknot with only a somewhat lower predicted stability, MFE −9.6 kcal/mol at 29 °C. The reduced stability of this predicted pseudoknot in the Lampung strains might reduce their efficiency of frameshifting to allow translation of the downstream, RdRp-encoding ORF2; however, this possibility, and whether any of these predicted structures are in fact the relevant ones, requires direct testing.

Concluding summary

New results obtained in this study include partial to complete genome sequences for eight IMNV strains collected in 2011–2012 (seven from Indonesia, one from Brazil), as compared to those from five strains collected in 2003–2007 (one from Indonesia, four from Brazil) and previously reported to GenBank. Moreover, two of the sequences reported herein encompass the whole IMNV genome, doubling the previous number of complete sequences for this virus as well as corroborating the coding organization of IMNV and its known nt- and aa-sequence motifs of functional significance, including ones for −1 ribosomal frameshifting and ORF1 polyprotein processing. Phylogenetic comparisons involving the eight new and five previously reported IMNV sequences provide new insights into the possible source and ongoing diversification of IMNV between Indonesia and Brazil, as well as within Indonesia. Additional sequencing of IMNV strains in coming years should be expected to clarify the trajectory of IMNV evolution and the propensity of IMNV for spread to other parts of the world.