9
views
0
recommends
+1 Recommend
0 collections
    0
    shares
      • Record: found
      • Abstract: found
      • Article: found
      Is Open Access

      Investigations into the presence of nidoviruses in pythons

      research-article

      Read this article at

      Bookmark
          There is no author summary for this article yet. Authors can add summaries to their articles on ScienceOpen to make them more accessible to a non-specialist audience.

          Abstract

          Background

          Pneumonia and stomatitis represent severe and often fatal diseases in different captive snakes. Apart from bacterial infections, paramyxo-, adeno-, reo- and arenaviruses cause these diseases. In 2014, new viruses emerged as the cause of pneumonia in pythons. In a few publications, nidoviruses have been reported in association with pneumonia in ball pythons and a tiger python. The viruses were found using new sequencing methods from the organ tissue of dead animals.

          Methods

          Severe pneumonia and stomatitis resulted in a high mortality rate in a captive breeding collection of green tree pythons. Unbiased deep sequencing lead to the detection of nidoviral sequences. A developed RT-qPCR was used to confirm the metagenome results and to determine the importance of this virus. A total of 1554 different boid snakes, including animals suffering from respiratory diseases as well as healthy controls, were screened for nidoviruses. Furthermore, in addition to two full-length sequences, partial sequences were generated from different snake species.

          Results

          The assembled full-length snake nidovirus genomes share only an overall genome sequence identity of less than 66.9% to other published snake nidoviruses and new partial sequences vary between 99.89 and 79.4%. Highest viral loads were detected in lung samples. The snake nidovirus was not only present in diseased animals, but also in snakes showing no typical clinical signs.

          Conclusions

          Our findings further highlight the possible importance of snake nidoviruses in respiratory diseases and proof multiple circulating strains with varying disease potential. Nidovirus detection in clinical healthy individuals might represent testing during the incubation period or reconvalescence. Our investigations show new aspects of nidovirus infections in pythons. Nidoviruses should be included in routine diagnostic workup of diseased reptiles.

          Related collections

          Most cited references19

          • Record: found
          • Abstract: found
          • Article: found
          Is Open Access

          Viruses Infecting Reptiles

          A large number of viruses have been described in many different reptiles. These viruses include arboviruses that primarily infect mammals or birds as well as viruses that are specific for reptiles. Interest in arboviruses infecting reptiles has mainly focused on the role reptiles may play in the epidemiology of these viruses, especially over winter. Interest in reptile specific viruses has concentrated on both their importance for reptile medicine as well as virus taxonomy and evolution. The impact of many viral infections on reptile health is not known. Koch’s postulates have only been fulfilled for a limited number of reptilian viruses. As diagnostic testing becomes more sensitive, multiple infections with various viruses and other infectious agents are also being detected. In most cases the interactions between these different agents are not known. This review provides an update on viruses described in reptiles, the animal species in which they have been detected, and what is known about their taxonomic positions.
            Bookmark
            • Record: found
            • Abstract: found
            • Article: not found

            Isolation, identification, and characterization of novel arenaviruses, the etiological agents of boid inclusion body disease.

            Boid inclusion body disease (BIBD) is a progressive, usually fatal disease of constrictor snakes, characterized by cytoplasmic inclusion bodies (IB) in a wide range of cell types. To identify the causative agent of the disease, we established cell cultures from BIBD-positive and -negative boa constrictors. The IB phenotype was maintained in cultured cells of affected animals, and supernatants from these cultures caused the phenotype in cultures originating from BIBD-negative snakes. Viruses were purified from the supernatants by ultracentrifugation and subsequently identified as arenaviruses. Purified virus also induced the IB phenotype in naive cells, which fulfilled Koch's postulates in vitro. One isolate, tentatively designated University of Helsinki virus (UHV), was studied in depth. Sequencing confirmed that UHV is a novel arenavirus species that is distinct from other known arenaviruses including those recently identified in snakes with BIBD. The morphology of UHV was established by cryoelectron tomography and subtomographic averaging, revealing the trimeric arenavirus spike structure at 3.2-nm resolution. Immunofluorescence, immunohistochemistry, and immunoblotting with a polyclonal rabbit antiserum against UHV and reverse transcription-PCR (RT-PCR) revealed the presence of genetically diverse arenaviruses in a large cohort of snakes with BIBD, confirming the causative role of arenaviruses. Some snakes were also found to carry arenavirus antibodies. Furthermore, mammalian cells (Vero E6) were productively infected with UHV, demonstrating the potential of arenaviruses to cross species barriers. In conclusion, we propose the newly identified lineage of arenaviruses associated with BIBD as a novel taxonomic entity, boid inclusion body disease-associated arenaviruses (BIBDAV), in the family Arenaviridae.
              Bookmark
              • Record: found
              • Abstract: found
              • Article: found
              Is Open Access

              Ball Python Nidovirus: a Candidate Etiologic Agent for Severe Respiratory Disease in Python regius

              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. RESULTS 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. TABLE 1  Lesions identified in ball pythons with respiratory tract disease BPno. Collection Age/sex Sinusitis/rhinitis Stomatitis/pharyngitis Tracheitis Mucousor caseousmaterial in airpassageways Thickenedlung Pulmonaryhemorrhage Pneumonia 1 1 Juvenile UDS a NE b No No No No No IPP c ; moderate 2 1 Juvenile UDS NE Moderate-severe; subacute Moderate Moderateto severe Moderateto severe No IPP; severe 5 3 Adult female NE No Mild tomoderatehyperplasia No No Yes IPP; moderate-severe 6 3 Adult female Mild; subacute Mild; subacute Mild tomoderatehyperplasia No No Yes IPP; moderate-severe 7 3 Adult female Mild; subacute Mild; subacute Mild tomoderatehyperplasia No No Yes IPP; moderate-severe 8 4 Juvenile female Moderateto severe;subacute Moderate withepithelialhyperplasia Moderateto severe;subacute No No No IPP; mild 9 5 Juvenile male No No Moderateto severe withepithelial hyperplasia Severe Moderate No IPP; mild 10 6 Adult female No No Moderate Moderateto severe Moderate No IPP; severe 11 7 Adult male NE No Necrotizing;severe No Moderate No IPP; severe a UDS, undetermined sex. b NE, not examined. c IPP, interstitial proliferative pneumonia. 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). FIG 1  Representative macroscopic lesions. (A) Common wild-type ball python (no. 1), Python regius. The palatine mucosa is both thickened and necrotic, and there is an accumulation of caseous material in the internal choanae. (B) Mojave variant ball python (no. 11), Python regius. The lung is thickened and edematous with abundant mucoid to caseous material (arrow) in air passageways. 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. FIG 2  Representative photomicrographs of diseased lung section. (A) Ball python 1 (BP1). Photomicrograph of a cross section of the lung. The smooth muscle (SM) bundle at the luminal end of each septum is hyperplastic and surrounded by inflammatory cells. As the septa approach the base of each faveolus, the septal wall between adjacent faveoli is of normal thickness. No inflammatory cells are seen in air passageways. H&E stain; bar = 500 µm. (B) BP1. Photomicrograph of a cross section of lung. At a higher power, numerous round cells (RC) have infiltrated the submucosa and overlying epithelial cells are hyperplastic (arrows). H&E stain; bar = 100 µm. (C) BP11. Photomicrograph of a cross section of lung. Severe proliferative interstitial pneumonia. Inflammation and epithelial hyperplasia (arrows) markedly thicken the trabecular septa and faveolar walls. Abundant mucinous exudate fills the lumina. H&E stain; bar = 200 µm. (D) BP11. Photomicrograph of a cross section of the lung at a higher magnification. Faveolar septum, severe proliferative interstitial pneumonia. The epithelium (EP) is hypercellular, thick, and densely packed with pseudostratified columnar cells, the majority of which produce mucus. Macrophages, lymphocytes, and granulocytes infiltrate (IN) the epithelium and submucosa. Faveolar capillaries are below the epithelial layer. H&E stain; bar = 50 µm. 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). FIG 3  Transmission electron photomicrographs of pulmonary epithelial cells of ball pythons, Python regius, with proliferative pneumonia. (A) BP6. Multiple bacilliform (arrows) and spherical nucleocapsid cores of virions are seen within the cytoplasm. Bar = 200 nm. (B) BP6. At a higher magnification, bacillary nucleocapsids are surrounded by granular cytoplasmic material (white arrows), presumed to be a component of the envelope. Bar = 100 nm. (C) BP6. Progression of immature forms in the cytoplasmic matrix to more mature forms (arrow) within a cytoplasmic vesicle. Bar = 200 nm. (D) BP2. Mature enveloped spherical virions (arrows) within a cytoplasmic vesicle. Spikes can be seen on the surface of several virions (arrowhead). Bar = 200 nm. (E) BP6. Spherical and filamentous mature virions (arrows) can be seen within cytoplasmic vesicles that are subjacent to microvilli projecting from the cell membrane Bar = 200 nm. (F) BP6. Spherical and bacillary forms of mature virions (arrows) can be seen in between microvilli and cilia of a pulmonary epithelial cell. Bacillary bacteria lacking a trilaminar cell wall (arrowheads) are also present. Scale bar corresponds to 1 µm. 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). FIG 4  Ball python nidovirus genome and replicase polyprotein organization. (A) A cartoon showing the genome organization of ball python nidovirus. Major open reading frames are indicated. RFS, position of putative ribosomal frameshift “slippery” sequence. (B) Position of Sanger sequencing PCR and RACE amplicons used to validate the NGS-based assembly (C) Concordance plot of read pairs mapping to assembly. Read pairs from snake no. 1 data set were mapped to the genome assembly using the Bowtie2 aligner. The inferred size of each library molecule is plotted. (D) Read coverage. Reads were aligned as in panel C, and the average number of bases supporting each position in 50-nt sliding windows of the alignment is indicated. The dip at approximately nt 5000 corresponds to a repeat-containing region with decreased mappability. (E) Conserved domains present in the replicase polyproteins (pp1ab). Domains were identified by searching the PFam database using the HMMER3 tool as described in the text and Table 2. Scale bars are indicated. 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. TABLE 2  Identification of functional domains in replicase polyprotein Domain a Description Position in replicasepolyprotein (AA) b Basis for assignment c Alignmentsupport d ADRP ADP-ribose-1′′-phosphatase (macrodomain) Not detected HMMER/Pfam Macro PKinase Protein kinase 2340–2494 HMMER/PFam PKinase 2 × 10−6 TM2 Transmembrane domain 4628–4737 TMHMM Mpro 3C-like main protease 4864–4997 HMMER/PFamPeptidase_S39 6 × 10−5 TM-3 Transmembrane domain 5121–5296 TMHMM RdRp RNA-dependent RNA polymerase 6200–6513 HMMER/PFam RdRP_1 1 × 10−12 Zn-Hel Zn-finger/5′-to-3′ helicase 6899–7178 HMMER/PFam Viral_helicase1 8 × 10−14 ExoN 3′-to-5′ exoribonuclease 7262–7423 HMMER/PFam NSP11 e 1 × 10−7 N7 Guanine-N7-methyltransferase Not detected HMMER/PFam NSP11 e NendoU Nidoviral uridylate-specific endoribonuclease 7709–7822 HMMER/PFam EndoU_bacteria 8 × 10−3 O-MT Ribose-2′-O-methyltransferase 7841–8081 HMMER/PFam NSP13 3 × 10−28 a Domain nomenclature according to reference 1. b Coordinates in amino acids of region within the BPNV pp1a or pp1ab polyprotein as defined by alignment boundaries. Closely juxtaposed TMHelix domains were considered to be part of a single TM domain. c The particular Pfam domain used to identify (or to exclude the presence of) each domain is indicated. Note that Pfam domain nomenclature does not necessarily correspond to current preferred nidovirus nomenclature. d The alignment’s E value from an HMMER3 search of the Pfam database. e The Pfam “NSP11” domain corresponds to the coronavirus nsp14 protein, which contains the ExoN and guanine-N7-MT domains. 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 [53]). 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). FIG 5  Analysis of subgenomic RNA species. (A) A cartoon showing the location of primers used to amplify regions of sgRNAs. Twenty-three kilobases have been removed for clarity. (B) Agarose gel electrophoresis analysis of PCR products obtained from reactions using the indicated primer pairs. Positions of molecular size standards (in bp) are indicated. (C) Sequences of PCR products spanning the junctions between upstream ORF sequences and 5′ UTR (leader) sequence. The top sequence in each triplet is the sequence of the genomic RNA (gRNA) beginning at base 150; the middle sequence is from the amplicons in panel B, corresponding to the indicated subgenomic mRNAs; the bottom sequence is that of the genomic RNA beginning at the indicated base. Regions of homology and putative start codons are underlined. For the ORF5a/b triplet, the start codons of both ORFs are indicated, and 180 nt of intervening sequence has been removed. 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. FIG 6  Detection of ball python nidovirus RNA is associated with clinical manifestation. (A) Relative levels of viral RNA were measured by qRT-PCR in case and control lung samples. Primers targeted ORF1a or ORF2 (S) sequences. Each row represents values from a single case or control animal of the indicated species. Levels were normalized to levels of cellular GAPDH mRNA and are represented as a heat map. Corallus annulatus and boa constrictor samples have been described previously (59). A polymorphism in the S sequence primer-binding site of snake 8 decreased amplification efficiency. (B) Relative levels of viral RNA in different tissues from infected snakes. For each snake, values were calculated and quantified as in panel A. Lung from a healthy control ball python served as a negative control. 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 98% global pairwise identity (77). Host-derived sequences were then filtered, first by using the BLASTn alignment tool (version 2.2.25+ [30]) to query a database of snake ribosomal and mitochondrial sequences and then by using the Bowtie2 alignment tool (version 2.0.0-beta7 [26]) 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 [42]). 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. SUPPLEMENTAL MATERIAL Table S1 Oligonucleotides used in this study Table S1, PDF file, 0.1 MB. Table S2 High-throughput sequencing summary Table S2, PDF file, 0.1 MB. Table S3 NCBI accession numbers of sequences of representative nidoviruses used for phylogenetic analyses Table S3, PDF file, 0.1 MB. Figure S1 Representative microscopic lesions in upper respiratory tract. (A) Common wild-type ball python (no. 7), Python regius. Photomicrograph of the oral mucosa, revealing a mild to moderate stomatitis. There is a focal area of necrosis (arrows) and heterophilic inflammation in oral mucosa. H&E stain; bar = 100 µm. (B) Cinnamon variant ball python (no. 8). Photomicrograph of the oral mucosa. There is diffuse necrosis of the oropharyngeal mucosa with numerous lymphocytes and fewer plasma cells and macrophages infiltrating the submucosa. H&E stain; bar = 75 µm. (C) Common wild-type ball python (no. 7). Photomicrograph of a nasal cavity revealing mild heterophilic and lymphocytic inflammation in sinus mucosa. (D) Cinnamon variant ball python (no. 8). Photomicrograph of the tracheal mucosa, revealing areas of hyperplasia and necrosis (arrow). The submucosa (SM) is severely and diffusely thickened, primarily with round cells. Cartilage (CA) surrounds the submucosa. H&E stain; bar = 75 µm. Download Figure S1, PDF file, 3.6 MB Figure S2 The putative AAAAAAC “slippery sequence” in the region of ORF1a/b overlap is well supported by Illumina and Sanger sequencing. The 55-nt region of overlap between ORF1a and ORF1b is shown, as are Sanger and Illumina sequences mapping to this region. Bases are numbered according to the overall genomic position. Sanger sequences are from directly sequenced PCR products. Twenty randomly selected Illumina reads (from 85× average coverage) are depicted, and disagreements to the consensus sequence are highlighted. Download Figure S2, PDF file, 1.1 MB Figure S3 Dot plot showing regions of similarity between BPNV and FHMNV replicase polyproteins. Lines in the dot plot represent 100-amino-acid windows with pairwise global alignment scores of ≥40 (BLOSUM50 matrix). A cartoon of BPNV pp1ab with predicted functional domains highlighted (see the text), with highlighting extended into the dot plot, is shown. Blue and orange colors correspond to the ORF1a and ORF1b encoded domains. Download Figure S3, PDF file, 1.5 MB Figure S4 Transmembrane domains and glycosylation sites in proteins predicted to be encoded by ORF2 to ORF6 of snake nidovirus. Domains and glycosylation sites were predicted as described in the text. Predicted protein lengths are indicated. Download Figure S4, PDF file, 0.1 MB Figure S5 Multiple alignment of conserved domain of replicase polyproteins of indicated representative nidoviruses: RdRp domain (A), helicase domain (B), and Mpro domain (C). Columns with >75% identity are highlighted. Sequence accession numbers are indicated. Abbreviations are defined in Table S3. Download Figure S5, PDF file, 0.3 MB Figure S6 Phylogenies of conserved replicase subunits of representative nidoviruses. Unrooted Bayesian phylogenies based on multiple sequence alignments of indicated conserved regions of replicase polyproteins. At select nodes, Bayesian posterior probabilities are indicated in percentages. The five major nidovirus families and subfamilies are indicated, as are select genera and the currently uncategorized BPNV and possum nidovirus taxa. Known host phyla are indicated in parentheses. The branch leading to BPNV is highlighted red for emphasis. Download Figure S6, PDF file, 0.4 MB
                Bookmark

                Author and article information

                Contributors
                kore.schlottau@fli.de
                Journal
                Virol J
                Virol. J
                Virology Journal
                BioMed Central (London )
                1743-422X
                17 January 2020
                17 January 2020
                2020
                : 17
                : 6
                Affiliations
                [1 ]Chemical and Veterinary Investigation Office, Westerfeldstraße 1, D-32758 Detmold, Germany
                [2 ]GRID grid.417834.d, Institute of Diagnostic Virology, Friedrich-Loeffler-Institut, ; Südufer 10, D-17493 Greifswald-Insel Riems, Germany
                [3 ]GRID grid.1016.6, Commonwealth Scientific and Industrial Research Organisation, ; Canberra, ACT 2601 Australia
                Author information
                http://orcid.org/0000-0002-3999-0393
                Article
                1279
                10.1186/s12985-020-1279-5
                6969405
                31952524
                8a61eb03-9095-4ccf-83a8-350d414ea562
                © The Author(s). 2020

                Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License ( http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

                History
                : 10 October 2019
                : 9 January 2020
                Categories
                Research
                Custom metadata
                © The Author(s) 2020

                Microbiology & Virology
                nidovirus,python,snake,respiratory disease,pneumonia,rt-qpcr,reptiles,high-throughput sequencing

                Comments

                Comment on this article