INTRODUCTION

The order Nidovirales encompasses a diverse group of viruses that includes significant veterinary and human pathogens (1–6). These viruses cause a variety of diseases that range from mild enteric infection to severe respiratory disease or hemorrhagic fever (7, 8). Examples of disease-causing nidoviruses include the severe acute respiratory syndrome (SARS) coronavirus, a number of other coronaviruses that cause typically mild respiratory disease in humans, and agriculturally important animal pathogens, such as equine arteritis virus, porcine reproductive and respiratory syndrome virus, and yellow head virus. Nidoviruses are characterized by their overall genome architecture, distinct pattern of gene expression, and presence of a conserved set of functional domains in their nonstructural polyproteins. The nidoviruses cluster into five major groups, which have been taxonomically categorized into four families: Arteriviridae, Roniviridae, Mesoniviridae, and Coronaviridae. Viruses in the Coronaviridae family (subfamilies Torovirinae and Coronavirinae) have the largest known RNA genomes, an attribute thought possible because of a virally encoded proofreading exonuclease (ExoN) that increases replication fidelity (9–12). Although nidoviruses are known to infect mammals, birds, fish, and crustaceans, no nonavian reptile nidovirus has been previously described.

Ball pythons (Python regius) have become one of the most popular types of reptiles sold and kept as pets (13). Native to West Africa, these snakes make popular pets because of their relatively modest size (≤1.5 m), docile behavior, and ease of care. Selective captive breeding has resulted in a tremendous variety of colors and patterns (morphs), many of which command high prices. Since the late 1990s, veterinarians have been aware of respiratory tract disease as a common syndrome affecting ball pythons. This syndrome is often characterized by pharyngitis, sinusitis, stomatitis, tracheitis, and a proliferative interstitial pneumonia. The clinical and epidemiological characteristics suggested an infectious etiology.

In this study, we investigated the pathology and etiology of this disease. We obtained case samples from 7 collections around the United States, performed necropsies, and collected multiple tissues for light microscopy and samples of lung for transmission electron microscopy (TEM). Although TEM of the lung suggested a viral etiology, traditional molecular diagnostic methods did not identify an agent. Metagenomic sequencing was used to identify and assemble the genome of a novel virus in the order Nidovirales. Here we describe clinical and pathological manifestations of this disease, ultrastructural findings, tissue tropism, disease association, and subgenomic RNA expression and analyze the genome of this virus in the context of related viruses.

Pathological findings.

Fixed and frozen samples from affected snakes from collections in Wisconsin, Texas, Florida, Oklahoma, and Pennsylvania were collected between 2006 and 2013 (see Materials and Methods and Table 1). Necropsies were performed on 9 snakes with clinical signs of respiratory disease, and multiple tissues were collected for histopathology; samples of lung from two snakes were collected for TEM.

Lesions were evident in the upper and lower respiratory tracts of all diseased animals (Table 1 and Fig. 1). Four snakes had stomatitis/pharyngitis, with one having caseous material in the choanae (sinusitis) (Fig. 1A). Three snakes had pulmonary hemorrhage, three had either mucoid or caseous material in air passageways, and four snakes had moderately to severely thickened lungs (Fig. 1B). At a microscopic level, four snakes had mild to severe stomatitis (see Fig. S1A in the supplemental material). Eight of nine snakes had tracheitis, with multifocal mild to severe changes, including infiltrates of heterophils, macrophages, and lymphocytes in the lamina propria, and in four of these cases, there was hyperplasia in the tracheal epithelium (see Fig. S1B).

In comparison to the histology of normal snake lung, all nine snakes had a proliferative interstitial pneumonia of various severities (Table 1). In some cases, inflammatory cells were primarily around the terminal smooth muscle bundle of each septum that projected into the central lumen of the lung (Fig. 2A and B). Hyperplastic alveolar cells completely proliferated over capillary beds that are normally distributed at the surface of the primitive faveoli (Fig. 2C and D). Epithelial cells were tall and columnar with abundant apical mucus production (mucous metaplasia). Periodic acid-Schiff (PAS) and alcian blue staining were positive, consistent with the presence of mucopolysaccharides. Heterophils, lymphocytes, and in some cases macrophages infiltrated the interstitium of the pulmonary septa (Fig. 2D) and often were seen in air passageways along with sloughed epithelial cells. In deeper areas of the faveolar spaces, the lung was histologically normal. In more severe cases, the proliferative changes occurred diffusely throughout the lung and were seen from the proximal to distal surfaces lining faveolar structures.

Three snakes had mild to marked bronchial epithelial hyperplasia, with mild to moderate smooth muscle hypertrophy of the terminal muscle bundles. In these cases, the bronchus-associated lymphoid tissue was mildly to moderately hyperplastic. Sections of lungs of two cases were stained with Warthin-Starry stain and Brown-Hopps Gram stain; there was no positive staining for bacteria between cilia and microvilli of pneumocytes overlying muscle bundles. Three of five nasal cavities that were examined revealed mild subacute rhinitis/sinusitis. The vomeronasal organ in one snake was almost obliterated by foamy macrophages. Ziehl-Neelsen acid-fast staining was negative for acid-fast bacteria, and methenamine silver staining was negative for yeast and hyphae. Three cases had mild to moderate stomatitis/pharyngitis. The changes seen in the respiratory tract were consistent with a viral infection.

Of the 9 necropsied snakes, lesions were also observed in other areas of the body. These included mild nonsuppurative encephalitis characterized by periventricular perivascular lymphocytic cuffing (in 1 of 9 snakes examined), mild to moderate acute nephritis (1/9), mild nephrosis (1/9), mild to moderate multifocal subacute dermatitis (1/9), mild acute salpingitis (1/9), and hepatic lipidosis (5/9). In the eye of one snake, there was moderate necrotizing conjunctivitis, with diffuse corneal ulceration, mild superficial keratitis, and mild inflammation and necrosis of the inner layer of the spectacle. One snake had mild bacterial colitis. It is unclear whether these other lesions were related to the observed lung pathology.

Ultrastructure of virus morphogenesis.

Two cases with moderate to severe diffuse proliferative pneumonia were selected for ultrastructural examination by TEM (snakes 2 and 6) (Table 1) (14). Virus-like particles were observed within pneumocytes lining the faveoli of both snakes. These particles were observed within ciliated and mucous epithelial cells overlying smooth muscle bundles at the tips of faveolar septa and alveolar type II cells lining faveolar spaces of both snakes (Fig. 3A). The particles corresponded to stages of a virus with both circular and bacillary (elongated rod-shaped) nucleocapsids and were seen within electron-lucent areas of the cytoplasm. At a higher magnification, bacillary nucleocapsids contained a lucent core that was surrounded by fine granular cytoplasmic material (Fig. 3B). The uncoated intracellular capsids measured 10 to 12 nm. Nucleocapsids lined up along membranes of cytoplasmic vesicles of uncertain origin, and as they progressed into a vesicle, a membrane coat was acquired (Fig. 3C). Cross sections of mature virions having a lipid envelope and surface spikes were also identified in cytoplasmic vesicles (Fig. 3D). Bacillary virions were observed within vesicles subadjacent to epithelial cell membranes on the free cell surface (Fig. 3E) and extracellularly between cilia and microvilli projecting into the air passageway (Fig. 3F). The mature particles measured approximately 50 by 180 nm. In some sections, filamentous bacteria (200 nm by 200 µm) were occasionally observed between cilia and closely associated with the cytoplasmic membrane (Fig. 3F).

Molecular diagnostics.

Pathology and EM findings were consistent with a viral etiology, so reverse transcription-PCR (RT-PCR) was used to attempt to identify an etiologic agent. Consensus primers targeted viruses in the families Paramyxoviridae (genus Ferlavirus), Rhabdoviridae, Reoviridae (genera Orthoreovirus and Aquareovirus), and Filoviridae (15–19). RNA was extracted from frozen lung tissue from animals 1 and 6 and used as a template for reactions, which were uniformly negative, though positive-control primers targeting conserved host sequences were positive.

Metagenomic sequencing and genome assembly.

Because consensus PCR approaches did not identify any candidate etiologic agents, we pursued an unbiased metagenomic shotgun sequencing strategy. Suitable frozen tissue was available for 8 of 9 snakes, and total RNA was extracted and converted into sequencing libraries (see Materials and Methods). The libraries were sequenced on an Illumina HiSeq 2500 instrument to an average depth of 2.9 million pairs of 100-nucleotide (nt) sequences per sample (see Table S2 in the supplemental material). We then used a stepwise analysis pipeline to efficiently identify possible pathogen-derived sequences. First, we removed low-quality and low-complexity sequences and trimmed adapter sequences from reads. For each library, identical reads were collapsed to yield sets of unique sequences. Next, we filtered out snake ribosomal sequences and sequences aligning to the Burmese python and boa constrictor genomes (20, 21). After these operations, approximately 2% of each data set remained.

We then searched for possible pathogen-derived sequences in what remained of the data sets. We performed BLASTx alignments using the remaining reads to search a database of virus protein sequences. In case samples but not in control samples, we observed sequences with similarity to viruses in the Nidovirales order. The data set from snake 2 had the most nidovirus-like sequences: 92 read pairs with BLASTx alignments to nidovirus protein sequences with an expect (E) value of ≤10−3 (see Table S2 in the supplemental material). These sequences were used to seed targeted de novo genome assembly using the PRICE metagenomic assembler (22), which produced a draft assembly of what appeared to be the complete viral genome of 34.5 kb (Fig. 4A). Sanger sequencing of the entire genome and rapid amplification of cDNA ends (RACE) were used to validate the assembly and to confirm that it contained the proper end sequences (Fig. 4B) (23–25). We used the Bowtie2 software alignment tool to retrospectively map paired-end reads to the draft genome assembly and found it to be well supported by deep sequencing data (26, 27). Coverage levels were even across the genome, and read pairs generally mapped concordantly (Fig. 4C and D). We found that the deep-sequencing-based assembly represented nearly the entire genome, except for 165 bases of 3′ untranslated region (UTR) sequence, which were determined by 3′ RACE. This genome sequence has been deposited in GenBank (see below).

Comparative analysis.

We next performed comparative and phylogenetic analyses to better understand this virus sequence in the context of related species. The genome of this virus is presumed to be positive-sense single-stranded RNA (ssRNA), and is 33,452 nt long. This is 1,766 nt longer than the genome of the beluga whale coronavirus SW1 (accession no. NC_010646.1), which was previously the longest described RNA genome, at 31,686 nt (28). The genome contains 8 positive-sense open reading frames (ORFs) longer than 400 nt (Fig. 4A). There are also 6 ORFs longer than 400 nt internal to larger ORFs, which are of unknown significance. The overall pattern of ORFs matches the pattern observed for related viruses (1–3, 9, 10, 29). There are 2 large (17,412 and 6,966 nt) 5′ ORFs that are predicted to encode nonstructural proteins (see below), designated ORF1a and ORF1b. The next ORF, ORF2, is predicted to encode the “S” or spike glycoprotein. Then, there are 5 “3′ ORFs,” designated ORF3 to ORF6, that are predicted to encode additional structural proteins (see below). The 5′ and 3′ UTRs are 1,017 and 916 nt, respectively, and RACE analysis corroborated the prediction that the 3′ end of the genome is polyadenylated.

We used several tools to study the predicted proteome of this ball python virus. We used the BLASTp alignment tool to search for similar sequences in the NCBI nonredundant protein database (nr) (30). We also used the more sensitive HMMER3 and HHPRED hidden Markov model-based alignment and structure prediction software tools to detect more distant homologies and identify functional domains within the replicase polyprotein (31, 32) (http://hmmer.org). We used TMHMM to predict transmembrane domains and the NetNGlyc and NetOGlyc tools to predict N- and O-linked glycosylation sites in protein sequences (33, 34) (http://www.cbs.dtu.dk/services/NetNGlyc/).

In other nidoviruses, the replicase gene, comprised of ORF1a and ORF1b, encodes nonstructural proteins involved in viral genome replication and modulation of host cell activities (1–3, 9, 10, 29). ORF1a is translated to produce the pp1a polyprotein. A fraction of translating ribosomes undergo a −1 frameshift near the end of ORF1a, resulting in the continuation of translation though ORF1b and the production of the pp1ab polyprotein. This virus is likely to utilize a similar mechanism of translation. ORF1a and ORF1b are offset by a −1 frame difference and overlap by 52 nt. This overlapping region contains a putative “slippery” sequence (AAAAAAC) that may be the site of frameshifting (Fig. 4A; see also Fig. S2 in the supplemental material) (35, 36). In this virus, pp1a is predicted to be a 5,803-amino-acid (aa) protein with a molecular mass of 663 kDa, and pp1ab is predicted to be 8,108 residues long with a mass of 925 kDa (Fig. 4E). As is the case for other nidoviruses, it is likely that these polyproteins are proteolytically cleaved into multiple functional domains (37).

Nidoviruses are characterized in part by the presence and organization of a set of functional subunits within their pp1ab replicase polyproteins (1, 5, 9, 10, 38, 39). We first queried the snake virus pp1ab sequence against the NCBI nr database using the BLASTp tool. The best alignments produced were to pp1ab sequences from viruses in the Torovirinae subfamily, with bafinivirus sequences producing slightly higher scoring alignments than did torovirus sequences. The best alignment was to the pp1ab sequence of fathead minnow nidovirus (FHMNV) (see Fig. S3 in the supplemental material) (40). This alignment covered nearly the entire second half of pp1ab in 2 contiguous pieces from residues 3973 to 5060 and 5429 to 8029 of snake virus pp1ab. In these regions, the polyproteins shared 22% and 29% pairwise identities, respectively. Thus, the overall organization and domain content of the second half of the snake virus’s pp1ab protein resembled most closely those of bafinivirus replicase polyproteins (see Fig. S3).

We next used the HMMER3 sequence analysis tool (http://hmmer.org) to search the PFam database for particular pp1ab domains, and for the most part, these domains were present and organized as expected given the overall similarity to bafinivirus replicase polyproteins (Fig. 4E and Table 2; see also Fig. S3 in the supplemental material) (1, 9, 32, 41, 42). The expected domains present in snake virus pp1ab included a 3C-like protease located between two predicted transmembrane regions near the end of pp1a (TM2-Mpro-TM-3), an RNA-dependent RNA polymerase (RdRp), a helicase (Zn-Hel), a 3′ to 5′ exoribonuclease (ExoN), a nidoviral uridylate-specific endoribunuclease (NendoU), and a ribose-2′-O-methyltransferase found at the carboxy terminus of pp1ab (O-MT) (Fig. 4E and Table 2) (12, 37, 38, 43, 44). The replicase polyproteins of some taxa of nidoviruses encode a guanine-N7-methyltransferase (1, 9, 43, 45). As is the case for viruses in the Torovirinae subfamily, snake virus pp1ab appears to lack this enzyme. Functional experiments will be necessary to validate the cleavage patterns, biochemical activities, and roles in the virus life cycle of these predicted replicase proteins.

One distinguishing feature of the snake virus replicase polyprotein is the apparent lack of an ADP-ribose binding “macro” domain that is present in nidoviruses in the family Coronaviridae (Fig. 4E and Table 2) (46–48). The activity of this domain is not essential for replication in vitro but may help counteract the innate immune response in vivo (49, 50). At approximately the same position of snake virus pp1a is a predicted protein kinase (PFam PKinase domain) that is not present in any other nidovirus as determined by HMMER analysis (HMMER3 E value, 1.5 × 10−6) (Fig. 4B and Table 2). This prediction is also supported by BLASTp analysis, with the best alignments from a search of the NCBI nr database with pp1ab residues 2300 to 2500 being to cellular serine/threonine protein kinases (lowest E value, 3 × 10−5).

In nidoviruses, the ORFs 3′ to the replicase genes encode structural and accessory proteins. The number and organization of these genes vary widely across nidovirus taxa, and in many cases there is little or no detectable similarity between sequences of more distantly related species (1, 9).

Nidovirus S proteins are involved in receptor binding and membrane fusion via their N- and C-terminal S1 and S2 subunits (1, 2, 51, 52). The best alignment from a BLASTp search of the snake virus protein against the nr database was to thrush coronavirus HKU12-600 spike glycoprotein (YP_002308497). This alignment had 18% pairwise identity over 333 aa in the C-terminal third of the protein, i.e., in the S2 membrane fusion domain (52). The putative S1 receptor-binding domain of this protein possesses no clear similarity to other proteins by BLAST or HMMER analysis (residues 25 to 625 queried). Consistent with its putative identity as a spike glycoprotein, this protein is predicted to contain several transmembrane (TM) domains, including one characteristically near the C terminus, and is predicted to be glycosylated (see Fig. S4 in the supplemental material).

The lipid bilayers of nidoviruses contain one or more proteins in addition to S, including the membrane (M) protein present in some nidoviruses (1, 2). The protein encoded by ORF4 is likely to be a homolog of the M protein of other large nidoviruses. This prediction is based on several attributes of the deduced protein sequence. First, the size of the snake virus protein (215 aa) is just under the range for other Coronaviridae M proteins (216 to 268). Second, the protein contains 3 predicted TM domains in the N-terminal half, a characteristic feature (see Fig. S4 in the supplemental material). Third, the protein possesses sequence similarity to the putative M protein of the white bream virus nidovirus that is detectable by BLASTp analysis of the nr database. The sequences are 21% identical over 116 aa, and the E value of this alignment is slightly lower than those of alignments to various bacterial protein sequences. Whether this putative M protein performs the same functional role as its distant relatives remains to be validated experimentally.

The 152-aa protein predicted to be encoded by ORF5a is similar to nucleocapsids (N) from other nidoviruses. This similarity is detectable by BLASTp, with the best such alignments from a search of the NCBI nr database being to porcine reproductive and respiratory syndrome virus (PRRSV) N proteins (22% pairwise identity over 125 aa) and by HMMER3 analysis, which identifies the Arteri_nucleo PFam domain in this sequence (E value, 5.6 × 10−9). Although the best local alignments from a BLASTp search of nr are to arterivirus N sequences, by pairwise global alignments the ORF5a-encoded N sequence is more similar to bafinivirus and torovirus N sequences (e.g., 24% pairwise identity to white bream virus N) than to arterivirus N sequences (e.g., 18% pairwise identity to PRRSV N). This protein is predicted to be basic, with an isoelectric point of 10.6, and to contain no transmembrane domains, consistent with a putative function as an RNA-binding nucleocapsid protein.

ORF3, ORF5b, and ORF6 encode proteins with no clear homology to other nidovirus proteins or to other viral or nonviral proteins. These proteins are all predicted to contain transmembrane domains and to contain multiple glycosylation sites (see Fig. S4 in the supplemental material). Along with S and M, these proteins are likely additional structural components of the virus particle. In some phyla of large nidoviruses, one of the 3′ ORFs encodes a hemagglutinin-esterase (HE) protein. The snake virus does not appear to encode an HE homolog. Similarly, no putative homolog of the small envelope (E) protein encoded by some nidoviruses was identified for this virus. E homologs have been identified in coronaviruses and arteriviruses but not in other nidovirus lineages (1).

Molecular characterization of subgenomic RNAs.

Nidoviruses generate subgenomic mRNA species (sgRNAs) that are competent for the translation of the downstream 3′ ORFs. In most of these viruses, these subgenomic mRNAs share 5′ and 3′ terminal sequences with the genomic RNA and are transcribed from minus-strand templates that are produced by a process termed discontinuous RNA synthesis (1, 2, 54). In fact, the prefix “nido” refers to the nested nature of these RNA species. In our deep-sequencing data sets, we detected pairs of reads where the first read mapped to the 5′ UTR and the second read mapped near the beginning of the downstream ORFs, suggesting that this virus might also generate sgRNA species by a similar process.

In order to characterize these RNAs, we performed PCR using primer pairs separated by more than 23 kb of genomic sequence but that would amplify putative sgRNAs (Fig. 5A) (55). PCR and sequencing of amplicons revealed sgRNA species for each of the predicted 3′ ORFs, except for ORF5a and -5b, which appear to share a single mRNA. ORF5a and -5b overlap by 251 nt, and their start codons are separated by 208 nt (Fig. 5B and C). The sgRNAs share a 5′-terminal “leader” sequence (approximately nt 1 to 170 of genome) (Fig. 5C). Short regions of homology are evident at the junction points of the sgRNAs (Fig. 5C). These data are consistent with this virus using a strategy of discontinuous RNA synthesis and translation similar to that employed by other nidoviruses (54).

Disease association and tissue tropism.

Using this complete genome sequence as a reference, we assessed whether other case and control samples contained sequences from similar viruses as determined by BLASTn (30). Nearly identical sequences were present in the fully filtered data sets of all 8 sequenced diseased snakes (see Table S2 in the supplemental material). The read coverage in the other data sets was not sufficient to generate complete genome assemblies, but the average pairwise nucleotide identity of reads aligning was ≥95% in all cases, suggesting that other diseased snakes were infected by closely related strains of the virus (see Table S2). Conversely, we never observed similar sequences in 57 data sets from tissues from unaffected ball pythons and from other snakes that were collected for other studies. These control data sets corresponded to various tissues (mainly not lung) from a number of snake species, including ball pythons (n = 2), other python species (n = 3), boa constrictors (n = 39), other boa species (n = 10), black-tailed rattlesnakes (n = 1), and California king snakes (n = 2).

To independently confirm the metagenomic sequencing results, we used quantitative RT-PCR (qRT-PCR) to measure viral RNA levels in these samples. Viral RNA was detected in all affected animals but not in control animals (Fig. 6A). Thus, detection of viral RNA by both metagenomic sequencing and qRT-PCR perfectly correlated with disease status, though larger studies will be necessary to confirm this association.

Data sets did not contain consistent evidence of other possible bacterial or viral pathogens. Data sets contained retrovirus-like sequences most similar to that of Python molurus endogenous retrovirus (ERV) sequences (56). These were present in nearly every case and control Python regius data set and most likely derived from ball python ERVs. Also, as with all metagenomic data sets from nonsterile sites, bacterial sequences were present. However, no particular bacterial taxon was associated with cases, in contrast to the nidovirus sequences, which were perfectly associated. Nevertheless, from these sequencing data sets alone, it is not possible to formally rule out a role in disease causality by other viral or bacterial pathogens.

We next sought to determine the distribution of viral load within the tissues of individual infected snakes. Frozen samples from multiple tissues were available for 4 snakes. We extracted RNA from these tissues and used qRT-PCR to quantify viral RNA relative to expression of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) mRNA (Fig. 6B). Consistent with the pathology observed, viral load was consistently highest in respiratory tract tissue. Viral RNA levels were at approximately the same level in intestine-derived RNA from the single snake for which this tissue was available (no. 11), a finding consistent with respiratory/enteric tropism of many viruses in the Coronaviridae family (7, 8). Viral RNA was also detectable in other tissues (liver, kidney, heart, spleen, and brain) at levels that were 3 to 5 orders of magnitude lower than in the lung. While these results support the notion that the respiratory tract is the primary location of viral replication, consistent with the disease pathology, additional systematic collection and analysis of samples from diseased animals will be necessary to definitively characterize tissue tropism for this virus.

Phylogenetic analysis.

We performed phylogenetic analyses to understand the evolutionary relationship between this and other nidoviruses. We obtained protein sequences from 33 representative species, and created multiple sequence alignments of the relatively conserved protease, RdRp, and helicase domains of pp1ab (see Table S3 and Fig. S5 in the supplemental material) (42). Bayesian and maximum-likelihood (ML) phylogenies based on individual domain and concatenated alignments had nearly identical topologies (Fig. 7; see also Fig. S6). In these phylogenies, the 5 major groups of nidoviruses (Torovirinae, Coronavirinae, Roniviridae, Mesoniviridae, and Arteriviridae) formed well-separated branches. Ball python nidovirus (BPNV) formed a monophyletic clade with the viruses in the Torovirinae subfamily of the Coronaviridae family. The snake virus is approximately equally distant from the two genera in the subfamily, the toroviruses, which infect mammals, and the bafiniviruses, which infect ray-finned fish (Fig. 7; see also Fig. S6). We noted in the protease alignment that two putative active site motifs previously identified by Ulferts et al. (57), one containing a histidine and a GX[S/C] G motif, aligned in all sequences. A third putative Asp/Glu-containing motif shifted drastically in alignments, even within Coronavirinae. Nevertheless, the topology of the tree generated using the protease alignment was in close agreement with other phylogenies. This analysis did not find a monophyletic grouping of the Torovirinae and Coronavirinae subfamilies of the family Coronaviridae. This paraphyly was found to be significant (P < 0.01) by log likelihood (Shimodaira-Hasegawa) testing (58).

Attempted isolation of virus in culture.

We attempted to isolate the virus in several snake cell lines, including the boa constrictor JK cell line and the viper VSW and VH-2 lines (59, 60). We prepared an inoculum from frozen infected tissue and applied it to these lines. Viral RNA was detectable in the inoculum but disappeared from the culture supernatants as media were replaced over the course of 14 days, and no cytopathic effects were observed. It is possible that different cells, alternate culture conditions, or inoculum from freshly harvested infected tissue would permit replication.

DISCUSSION

In this report, we describe a severe respiratory disease affecting ball pythons, and we have identified a nidovirus as a candidate etiologic agent by means of association. We collected cases from around the United States and performed a variety of pathological and diagnostic analyses. Respiratory tract pathology was consistent with a viral etiology, as was the presence of virus-like particles in lung epithelium of affected snakes. Using unbiased deep sequencing, we identified and assembled the genome of a previously uncharacterized nidovirus. This virus shares several attributes with other nidoviruses but has several distinguishing characteristics, including its record-setting genome size and reptile host.

The findings presented here raise a number of questions related to the biology of this virus and its relationship to disease. Foremost among these is whether this virus causes the observed respiratory pathology. Several observations suggest a causal relationship. First, detection of viral RNA correlates with clinical signs. Second, the virus load is most profound in the respiratory tract, the site of disease. Third, other animal nidoviruses cause severe respiratory disease. Nevertheless, formal evidence of disease causality remains to be demonstrated. Although virus isolation attempts have yet to succeed, experimental infections using material prepared directly from infected tissues could fulfill Koch’s postulates.

The evolutionary relationship between ball python nidovirus and other nidoviruses is partially clear. In analyses based on the highly conserved replicase subunits, this virus clearly clusters with viruses in the subfamily Torovirinae, the toroviruses found in mammals and the bafiniviruses found in ray-finned fish, with 100% Bayesian posterior probability and ML bootstrap support (Fig. 7; see also Fig. S6 in the supplemental material). This suggests that at least for the core replicase gene, the snake virus shares a common evolutionary origin with these mammal and fish viruses. Although it has not been proven that the virus sequence corresponds to the particles observed in electron micrographs, the ultrastructural appearance and size of the observed virus are notably similar to those of virions of white bream virus in the genus Bafinivirus (61, 62). We propose that ball python nidovirus establish a new genus in the Torovirinae named Barnivirus, for bacilliform reptile nidovirus, a naming convention borrowed from the genus Bafinivirus.

While some studies have found support for the monophyly of Torovirinae and Coronavirinae, in accordance with the current International Committee on Taxonomy of Viruses (ICTV) status of Coronaviridae, our analysis did not find phylogenetic support for this grouping, in agreement with other recent analyses (1, 40, 42, 63). The paraphyly of Coronaviridae was supported by log likelihood testing. The availability of a more complete representation of existing species for comparison results in greater phylogenetic resolution (64). Significant errors can occur in phylogenetic analyses due to incomplete taxum sampling, even if very large sequence length is assessed (65). This further emphasizes the need for exploration of viral diversity and increased sampling of the Nidovirales.

Current standards in herpetoculture encourage disease transmission and may foster evolution of increased pathogen virulence. Captive snake breeding operations typically operate at high stocking densities, and breeders commonly attend trade shows, where animals from different sources are juxtaposed. In addition, animals from geographically and ecologically diverse areas are commonly imported and mixed with minimal quarantine. These practices increase pathogen exposure and lower barriers to transmission. The outcome is evident in the case of farmed turtles, which typically have very high Salmonella carriage rates, in contrast to the low prevalence in wild turtle populations (66–68). It is also possible that these practices select for increased virulence, as has been the case for feline calicivirus (FCV), which has independently evolved multiple times from a pathogen that causes a relatively benign upper respiratory disease to one that causes hemorrhagic disease with high mortality (69). This increase in FCV virulence has been documented only in high-density environments, such as shelters. Surveillance of wild ball pythons in Africa will shed light on whether a similar phenomenon has occurred here. Further investigation of the pathogens of reptiles and improved biosecurity practices and disease surveillance of the herpetoculture industry are indicated.

In addition to disease causality, a number of questions related to the biology of ball python nidovirus and to the epidemiology and natural history of this disease remain open. The precise host range of this virus remains undefined. It is not clear that ball pythons are the primary natural host for this virus, nor whether it can replicate or cause disease in other snakes, other reptiles, or other animals. The routes of transmission for this virus remain unknown. Additionally, the prevalence of this disease is yet to be determined, although it is apparently widespread, since cases were collected from geographically disparate locations in the United States from 2006 to 2013. Additional studies will help to answer these central questions and uncover additional details about the biology of this putative pathogen. The current study is another example where the application of genomic techniques to the study of infectious disease in reptiles and other less well studied organisms has identified new and unexpected relatives of human pathogens (59, 70).

Collections, animals, and pathological evaluations.

Carcasses and tissues of 9 ball pythons (BP) with clinical signs of respiratory tract disease (stomatitis, pharyngitis, rhinitis, glossitis, tracheitis, bronchitis, and pneumonia) from seven collections in Florida, Oklahoma, Pennsylvania, Texas, and Wisconsin were submitted for necropsy and light microscopic evaluation. Transmission electron microscopy (TEM) and/or molecular diagnostics were performed on selected cases.

(i) Collection 1.

Two snakes from a private collection in Texas (BP1 and -2) were submitted in 2007 for pathological evaluations. The owners reported that affected snakes exhibited acute signs of respiratory disease characterized by audible wheezing along with a clear mucoid nasal discharge, notable “cough-like” forced expiration, and excessive amounts of foamy mucus in the oral cavity. BP1 was euthanatized and immediately necropsied. BP2 was found dead and was refrigerated for 2 days prior to necropsy. Lung was collected from both snakes, and portions were frozen or placed in neutral buffered 10% formalin (NBF) and processed for light microscopy. Lung was sectioned at 6 µm and stained with hematoxylin and eosin (H&E). A portion of lung of BP1 that was fixed in NBF was cut into 1-mm cubes and transferred to 4% paraformaldehyde–2.5% glutaraldehyde in 0.1 M sodium cacodylate buffer, pH 7.24, and processed for TEM (see below). Molecular diagnostics were performed on frozen tissues from two snakes (snakes 3 and 4) from a second collection in Texas (collection 2), but they were negative and are not discussed further in this report.

(ii) Collection 3.

Three ball pythons (5, 6, and 7) from a private collection in Wisconsin were submitted to a veterinarian due to signs of respiratory disease. Ball python 5 was found dead and submitted for necropsy. The following tissues were collected, fixed in NBF, and processed for light microscopy: lungs, kidney, liver, stomach, spleen, and ovary with oviduct. Ball pythons 6 and 7 had acute onset of respiratory signs over a 2- to 3-day period characterized by an increase in respiratory rate and effort. There also was elevation of the cranial third of the snake’s body (from the ground), ultimately resulting in opisthotonos (stargazing). A decision was made to euthanatize and necropsy BP6 and -7. The following tissues were collected from each snake, fixed in NBF, and subsequently embedded in paraffin: trachea, lungs, heart, esophagus, stomach, intestines, liver, pancreas, kidneys, oviduct, brain, spinal cord, nasal sinuses, vomeronasal organ, and skin. The paraffin-embedded tissues were sectioned at 6 µm, processed routinely for light microscopy, and stained with hematoxylin and eosin. Sections of the vomeronasal organ of one case were stained with Ziehl-Neelsen acid-fast stain, which was negative for acid-fast bacteria; methenamine silver staining was negative for yeast and hyphae. A portion of lung from BP6 was fixed in 4% paraformaldehyde–2.5% glutaraldehyde in 0.1 M sodium cacodylate buffer, pH 7.24, and processed for TEM (see below).

(iii) Collection 4.

A full necropsy and pathological evaluation were performed on a 1-year-old female ball python (BP8) from a private collection in Texas that died after a history of respiratory disease for 1 week and respiratory distress for 2 days. The snake was from a group of four snakes. There was a history of respiratory disease in the group, and previously one other snake had died. The snake exhibited marked dyspnea with crackles and wheezes audible without a stethoscope. Other clinical findings included hypersalivation and a swollen, erythematous glottis. Culture of a tracheal swab isolated Shewanella putrefaciens. No response to treatment with the antimicrobial enrofloxacin was noted. The following tissues were analyzed: lungs, kidney, liver, stomach, small intestine, large intestine, heart, brain, oral cavity, pharynx, nasal cavity, trachea, eye, spleen, and ovary with oviduct were placed in NBF and processed routinely for light microscopy. Six-micrometer sections of paraffin-embedded tissues were stained with hematoxylin and eosin.

(iv) Collections 5 and 6.

Two BPs, one from private collection 5 in Florida (BP9) and the other from private collection 6 in Oklahoma (BP10), were submitted to an exotic animal practice for evaluation due to signs of respiratory disease and anorexia. BP9 was part of a small breeding collection where very few sick animals were previously observed. BP9 was euthanized, and the following tissues were collected during necropsy, fixed in NBF, and processed for light microscopy: trachea, lung, kidney, salivary gland, liver, brain, and small intestine. Portions of trachea, liver, and lung were collected, frozen at −20°C, and subsequently submitted for molecular diagnostics. Six-micrometer sections of paraffin-embedded tissues were stained with hematoxylin and eosin.

The second case (BP10) was a 5-year-old female common BP from a large collection (500-plus) of ball pythons in Oklahoma. Prior to 2009, no significant health issues were recognized in this collection. Starting in 2009 several ball pythons showed signs of respiratory disease, and from 2010 to early 2011, three ball pythons died with clinical signs of respiratory disease. Ball python 10 had signs of respiratory disease and was euthanized for pathological evaluation. The following tissues were collected, fixed in NBF, and processed for light microscopy: brain, kidney, large intestine, liver, lung, small intestine, and trachea. Paraffin-embedded tissues were sectioned at 6 µm and were stained with hematoxylin and eosin. Portions of lung were collected, frozen at −20°C, and subsequently submitted for molecular diagnostics.

(v) Collection 7.

Ball python 11 was a 6-year-old adult from a private collection in Pennsylvania. The owner had about 74 snakes in total, all of which were ball pythons. Seven snakes in close proximity to one another from his collection became sequentially ill with clinical signs of respiratory disease, including dyspnea and open-mouth breathing, All ill snakes were immediately quarantined in a separate room from the general population. Of the 7 snakes that became ill, 6 died within a few days of first manifesting signs of respiratory disease. One of the dead snakes was submitted to a private practitioner for necropsy. Portions of trachea, lung, kidney, liver, integument, and small intestine were fixed in NBF and submitted to a pathology service processed for light microscopy. Paraffin-embedded tissues were sectioned at 6 µm and were stained with hematoxylin and eosin, alcian blue, and PAS. Portions of lung, liver, and kidney were collected, frozen at −80°C, and subsequently submitted for molecular diagnostics.

Ethics statement.

No live animals were used in this study. Tissues and carcasses were from animals from private collections. Snakes either died or were euthanized by the attending veterinarian due to the severity of clinical signs and poor prognosis as described above. Tissues or carcasses were collected during necropsy by attending veterinarians and submitted for the purposes of this study, with the owners’ consent. The acquisition of fresh, frozen, and formalin-fixed tissue samples was authorized under University of Florida Institutional Animal Care and Use Committee Protocol A116.

Transmission electron microscopy.

Formalin-fixed lung samples from BP1 and BP6 were cut into 1-mm cubes and transferred to 4% paraformaldehyde–2.5% glutaraldehyde in 0.1 M sodium cacodylate buffer, pH 7.24. Fresh lung from BP6 was also cut into small cubes and fixed directly in 4% paraformaldehyde–2.5% glutaraldehyde in 0.1 M sodium cacodylate buffer, pH 7.24. Fixed tissues were processed with the aid of a Pelco BioWave laboratory microwave (Ted Pella, Redding, CA, USA). The samples were washed in 0.1 M sodium cacodylate buffer, pH 7.24, postfixed with buffered 2% OsO4, water washed, and dehydrated in a graded ethanol series, 25%, 50%, 75%, 95%, and 100%, followed by 100% acetone. Dehydrated samples were infiltrated in graded acetone-Spurrs epoxy resin (Electron Microscopy Sciences, Hatfield, PA) 30%, 50%, 70%, and 100% and cured at 60°C. Cured resin blocks were trimmed, thin sectioned (70 nm), and collected on Formvar-coated copper grids, poststained with 2% aqueous uranyl acetate and Reynolds lead citrate. Sections were examined with a Hitachi H-7000 TEM (Hitachi High Technologies America, Inc., Schaumburg, IL), and digital images were acquired using a Veleta 2k by 2k camera and iTEM software (Olympus Soft-Imaging Solutions Corp., Lakewood, CO).

Clinical molecular diagnostics.

Portions of fresh frozen lung specimens were selected from several of the above cases for identification of nucleic acid sequences for several viruses implicated as causative agents of respiratory tract disease of snakes. RNA was extracted from lung of BP1 and tested by using previously described consensus PCR methods for the genus Ferlavirus in the Paramyxoviridae, the Rhabdoviridae, and the genus Orthoreovirus in the Reoviridae (71–73). Additionally, RNA was extracted from portions of lungs of BP6 and -7 and tested for gene sequences for members of the family Rhabdoviridae and Filoviridae. For molecular diagnostic detection of filoviruses (Marburg viruses and Ebola viruses), RNA was sent to the Viral Special Pathogens Branch, Centers for Disease Control and Prevention. Two independent panfilovirus RT-PCR assays were utilized as previously described (74, 75) using SuperScript III reverse transcriptase (Life Technologies) according to the manufacturer’s instructions. Briefly, RNA was reverse transcribed for 30 min at either 45°C or 50°C, respectively, and then denatured for 2 min at 94 to 95°C. For RNAs amplified according to the protocol described elsewhere (74), reaction mixtures were thermocycled 40 times at 94°C for 15 s, followed by successive incubations at 45°C and 68°C for 30 s each. For RNAs amplified according to the previously described protocol (75), reaction mixtures were treated identically except that the annealing and extension temperatures were 53°C and 72°C degrees, respectively. Positive- and negative-control reactions were carried out in parallel to validate the performance of the testing.

Library preparation and sequencing.

Total RNA was extracted from frozen tissues for sequencing library generation and molecular analysis. RNA was extracted as previously described (59). Two hundred fifty nanograms of RNA was added to 20 µl RT reaction mixtures containing 200 pmol oligonucleotide MDS-286 (random hexamer), 1× reaction buffer, 5 mM dithiothreitol, 1.25 mM (each) deoxynucleoside triphosphates (dNTPs), and 100 U SuperScript III (Life Technologies). The sequences of all oligonucleotides are listed in Table S1 in the supplemental material. Reaction mixtures were incubated for 60 min at 42°C and then for 15 min at 70°C. Then, 1 U RNase H (NEB) diluted in 5 µl 1× RT buffer with 100 pmol MDS-286 was added to reaction mixtures, which were incubated for 30 min at 37°C and then 94°C for 2 min. Then, primer extension was used to convert single-stranded cDNA to double-stranded DNA (dsDNA): 2.5 U Klenow DNA polymerase (3′ to 5′ exo-; NEB) diluted in 5 µl 1× RT buffer with 2 mM (each) dNTP was added to reaction mixtures, which were incubated at 37°C for 15 min. DNA was purified using DNA clean and concentrator columns (CC-5; Zymo Research) according to the manufacturer’s protocol and using a 3:1 ratio of binding buffer/sample. The DNA concentration was determined fluorometrically, and 10 ng was used as a template in 10 µl Nextera transposon reactions, which contained 1× buffer and 0.5 µl transposase mix (Illumina), which were incubated at 55°C for 10 min and then placed on ice. Tagmented DNA was cleaned with CC-5 columns (4:1 buffer ratio), eluted in 20 µl H2O, and used as a template in PCRs that added full-length bar-coded adapters to library molecules. PCR mixtures contained 1× Kapa real-time library amplification PCR mix (Kapa Biosystems), 0.33 µM (each) primers MDS-143 and -445 and 0.017 µM (each) adapter1 and adapter2 bar-coded primers, and 5-µl library template. Thermocycling conditions were 72°C for 3 min, followed by 95°C for 30 s, and then 5 cycles of 98°C for 10 s, 63°C for 30 s, and 72°C for 3 min. Relative concentrations of libraries were then quantified in quantitative PCR (qPCR) reaction mixtures containing 10 mM Tris, pH 8.6, 50 mM KCl, 1.5 mM MgCl2, 5% glycerol, 0.08% IGEPAL CA-630, 0.05% Tween 20, 0.2 mM (each) dNTP, 1× Sybr green (Life Technologies), and 0.5 U Taq polymerase, and 0.1 µM (each) primers MDS-143 and -445. Equivalent amounts of each sample were then pooled, and the pool was cleaned with a DNA CC-5 column according to the manufacturer’s protocol. The pooled libraries were then size selected (350 to 450 nt) using a LabChip XT device (Caliper), and purified again using a CC-5 column. Size-selected pooled libraries were subjected to a final round of amplification in PCR mixtures containing 1× Kapa real-time library amplification mix, 200 pmol each MDS-143 and -445, and 5 µl library template. Reaction conditions were 98°C for 1 min and 16 cycles of 98°C for 15 s, 63°C for 30 s, and 72°C for 3 min. The pooled libraries were finally purified using a DNA CC-5 column and quantified spectroscopically and by using the Illumina library quantification kit (Kapa Biosystems). Sequencing was performed on an Illumina HiSeq 2500 instrument.

Sequence analysis.

The sequence analysis pipeline was similar to that employed previously (59, 76). Low-quality sequences and low-complexity sequences (defined as having a ratio of the length of the Lempel-Ziv-Welch compressed sequence to the uncompressed length of less than 0.46) were removed from further analysis. The first 5 bases of each sequence were trimmed, as was the last base. The CD-HIT sequence clustering tool was then used to collapse reads with >98% global pairwise identity (77). Host-derived sequences were then filtered, first by using the BLASTn alignment tool (version 2.2.25+) to query a database of snake ribosomal and mitochondrial sequences and then by using the Bowtie2 alignment tool (version 2.0.0-beta7) to query databases composed of draft assemblies of the Burmese python and boa constrictor genomes (20, 21). Sequences aligning with an expect value of less than 10−12 (BLASTn) or with –local mode alignment scores greater than 86 (Bowtie) were filtered. Similarly, sequences that aligned to the Illumina adapter sequences (see Table S1 in the supplemental material) or to PhiX-174 control sequence were removed. This filtering removed on average 98% of sequences (see Table S2). The remaining sequences were searched against custom databases of viral protein sequences using the BLASTx alignment tool. Sequences aligning to any viral protein sequence with an expect value of less than 0.25 were further examined. False positives were reduced by using BLAST to align putative viral sequences to the NCBI nonredundant nucleotide (nt) and protein (nr) databases. Only sequences whose best hit and whose pair’s best hit were still to viral sequences were retained. The PRICE de novo targeted genome assembler was used to generate initial contiguous virus sequences (22). The reference assembly (NCBI accession no. KJ541759) was assembled using reads from a single snake (no. 2); other animals lacked sufficient coverage to assemble the complete genome. To validate this assembly and generate coverage and pairwise identity data, reads were aligned to Sanger validated assemblies using the BLASTn algorithm. Sequencing data have been deposited in the UCSF Integrated Data Repository.

Sanger sequencing and RACE.

The sequencing-based virus genome assembly was validated using RACE and Sanger sequencing. PCR mixtures contained 1× iProof High-Fidelity DNA polymerase mix (Bio-Rad), 0.5 µM (each) primer, and 2 µl diluted cDNA, prepared as described above. Thermocycling conditions were 98°C for 30 s and then 40 cycles of 98°C for 10 s, 60°C for 20 s, and 72°C for 1 min, with a final 5-min 72°C extension. Primers were designed to tile across the assembly in overlapping ~1,200-bp amplicons (Fig. 4B) (78). Amplicons were verified on agarose gels and sequenced directly. In some cases, amplicons were gel purified with the Zymoclean kit (Zymo Research), cloned into the pCRBlunt vector (Invitrogen) according to the manufacturer’s protocols, and sequenced. 5′ and 3′ RACE was used to validate end sequences, essentially as described previously (24, 25), with primers listed in Table S1 in the supplemental material. Multiple RACE amplicons were cloned and sequenced as described above. The complete and validated genome sequence has been deposited with the NCBI (see below). Subgenomic RNA junctions were amplified and sequenced similarly.

Quantitative PCR.

To measure relative levels of viral RNA in samples, RNA was extracted and reverse transcribed as described above for library preparation. cDNA was diluted 1:10 in 10 mM Tris–0.1 mM EDTA, pH 8. qPCR reaction mixtures contained 10 mM Tris, pH 8.6, 50 mM KCl, 1.5 mM MgCl2, 5% glycerol, 0.08% IGEPAL CA-630, 0.05% Tween 20, 0.2 mM (each) dNTP, 1× Sybr green (Life Technologies), 0.5 U Taq polymerase, and 0.25 µM (each) primer (see Table S1 in the supplemental material). Viral RNA levels were normalized to levels of cellular GAPDH (see Table S1). Reaction efficiencies were determined using serial dilutions of templates (79).

Tissue culture and virus isolation.

Cells were cultured as described previously (59). To prepare inocula, frozen tissues were thawed on ice, and 1-g portions were minced with scalpels and disrupted by Dounce homogenization in ice-cold minimum essential medium (MEM) (Gibco) supplemented with 25 mM HEPES buffer, pH 7.4. Homogenate was clarified by centrifugation at 10,000 × g for 1 min and filtered through a 0.45-µm filter. Filtrate was diluted 1:10 in MEM plus HEPES and added to cultures of 80% confluent cells. Culture medium was harvested and replaced after 20 h and every 2 to 3 days thereafter for an additional 14 days and stored at −80°C until further processing. RNA was isolated from clarified culture supernatant using the ZR viral RNA kit (Zymo Research), according to the manufacturer’s protocol. RNA was reverse transcribed and analyzed by qRT-PCR as described above.

Comparative and phylogenetic analyses.

ORFs in the viral genome were identified using the Geneious software program (version 6.1; BioMatters). Homologs of predicted protein sequences were detected using the BLASTp tool (version 2.2.25+) to search the NCBI nonredundant protein database (nr) (30), and by using the HMMER3 alignment tool (version 3.1b1) to search the PFam database (http://hmmer.org) (32). For HMMER analysis of replicase polyprotein, sequence was split into 400-aa sliding windows offset by 60 aa. For sequences with no detectable similarity by those measures, the HHpred homology detection and structure prediction tool (version 2.0) was also used (31). TMHMM (version 2.0) was used to predict transmembrane domains. The NetNGlyc (1.0) and NetOGlyc (4.0) tools were used to predict N- and O-linked glycosylation sites in protein sequences (33, 34) (http://www.cbs.dtu.dk/services/NetNGlyc/). In cartoons depicting sequence features, all features are to scale and accurately positioned.

Sequences from a set of 33 representative nidoviruses based on a previously published analysis were downloaded (see Table S3 in the supplemental material). Conserved protease, RdRp, and helicase domains in replicase polyproteins were identified using HMMER3 and BLAST alignments. Multiple sequence alignments of these domains were created using MAFFT version 7 with default parameters (see Fig. S5) (80). Amino acid substitution models were evaluated for each region using corrected Aikake information criteria in the software program ProtTest 3.2.2 (81). LG was identified as the best-fit model (82, 83). Bayesian analyses of each alignment were performed using the software program MrBayes 3.2.2 (84) on the CIPRES server (85), with gamma-distributed rate variation and a proportion of invariant sites, and amino acid models were also assessed by model jumping. The first 10% of 1,000,000 iterations were discarded as a burn-in, based on examination of trends of the log probability versus generation. Two independent Bayesian analyses were run to avoid entrapment on local optima.

Maximum-likelihood (ML) analysis was performed using the RAxML software tool on the CIPRES server (86), with gamma-distributed rate variation and a proportion of invariant sites. The amino acid substitution model with the best fit using ProtTest was selected. Gill-associated virus was used as an outgroup. Bootstrap analysis was used to test the strength of the tree topology (1,000 resamplings) (87). Numbers of bootstrap replicates were determined using the stopping criteria of Pattengale et al. (88).

To test for paraphyly of the Coronaviridae, likelihood ratio testing was conducted (58). Trees were constrained to Coronaviridae monophyly (including ball python nidovirus with the Torovirinae) and compared by the Shimodaira-Hasegawa test to the best unconstrained tree identified in RAxML.

Nucleotide sequence accession number.

The complete and validated sequence assembled using reads from snake no. 2 has been deposited with the NCBI under accession no. KJ541759.