Bioinformatics and Molecular Analysis of the Evolutionary Relationship between Bovine Rhinitis A Viruses and Foot-And-Mouth Disease Virus.

Bovine rhinitis viruses (BRVs) cause mild respiratory disease of cattle. In this study, a near full-length genome sequence of a virus named RS3X (formerly classified as bovine rhinovirus type 1), isolated from infected cattle from the UK in the 1960s, was obtained and analyzed. Compared to other closely related Aphthoviruses, major differences were detected in the leader protease (L(pro)), P1, 2B, and 3A proteins. Phylogenetic analysis revealed that RS3X was a member of the species bovine rhinitis A virus (BRAV). Using different codon-based and branch-site selection models for Aphthoviruses, including BRAV RS3X and foot-and-mouth disease virus, we observed no clear evidence for genomic regions undergoing positive selection. However, within each of the BRV species, multiple sites under positive selection were detected. The results also suggest that the probability (determined by Recombination Detection Program) for recombination events between BRVs and other Aphthoviruses, including foot-and-mouth disease virus was not significant. In contrast, within BRVs, the probability of recombination increases. The data reported here provide genetic information to assist in the identification of diagnostic signatures and research tools for BRAV.

Introduction

Bovine respiratory disease (BRD) is a leading cause of death of feedlot cattle in the USA. Annual costs associated with BRD total more than one billion dollars. Environmental stress, compromised host immunity, and virus infection predispose the animal to bacterial lung infection (broncho-pneumonia).1 Primary virus infections damage the respiratory tract and subdue the host immune system. The respiratory tract commensal bacteria cause opportunistic secondary infections and lung pneumonia under immunecompromised conditions (caused by viral infection and other factors). In fact, previous studies have reported bovine viral diarrhea virus, bovine respiratory syncytial virus, bovine herpesvirus 1, parainfluenza 3 virus, adenovirus, and bovine coronavirus from cattle with BRD.2–6 Although more than 40 years ago, studies concluded that bovine rhinitis viruses (BRVs) (then known as bovine rhinoviruses) cause mild respiratory disease in cattle,7–10 it is a surprise that BRVs are sparsely studied. This is partly due to the fact that the virus grows poorly in cell culture and is known to cause only self-limiting mild respiratory illness.

Based on the genomic organization and serologic and molecular characteristics, bovine rhinitis A virus (BRAV) and bovine rhinitis B virus (BRBV) along with foot-and-mouth disease virus (FMDV) and equine rhinitis A virus (ERAV) are classified as four species in the Aphthovirus genus of the family Picornaviridae.11,12 The Picornaviridae family consists of some of the most devastating human and animal pathogens, all of which are single-stranded positive-sense RNA viruses (+ssRNA). Presently, the following sequences of BRV isolates are available in GenBank: BRAV (140032-1, KP236129; SD-1, KP236128; H-1, JN936206; and BSRI4, KP264974) and BRBV (EC 11, EU236594; 140032-2, KP236130; BSRI1, KP264980; BSRI2, KP264976-KP264979; and BSRI3, KP264975).7,8,13

In the new era of deep sequencing, metagenomic profiling of biological samples has yielded a greater diversity of pathogens than was previously detectable. For example, recently, Ng et al.7, using metagenomic approach, detected bovine influenza virus, bovine adeno-associated virus, bovine parvovirus, picobirnavirus, and multiple strains of BRAV and BRBV along with earlier reported viruses from BRD cattle samples. In another independent study, Hause et al.8 showed that BRAV and BRBV commonly circulate in BRD cattle. These findings raise a concern that similar to other BRD-associated viruses, primary infection with BRVs could facilitate bacterial infection. Alternatively, BRVs may themselves cause BRD in some cases.

Previously, we characterized the near full-length genome length sequence of BRBV strain EC11 (EU236594, previously known as bovine rhinovirus type 2) and compared virus-encoded RNA elements and proteins with the related Aphthoviruses.13,14 BRVs can be safely handled in a biosafety level 2 laboratory and have thus received more attention after Osiceanu et al.15 proposed BRBV as a surrogate model for FMDV anti-viral development, the latter being a virus that requires strict biosafety level 3 conditions for any manipulation. The reason for the strict laboratory guidelines for handling FMDV stems from the fact that it is the etiological agent of an economically devastating disease and is one of the most infectious animal viruses known. Like BRV and other Aphthoviruses, of which FMDV is the prototypic member, FMDV is a +ssRNA virus with a large open reading frame (ORF) flanked by two highly structured 5′ and 3′ untranslated regions (UTRs). As with other Aphthoviruses, the FMDV ORF encodes a polyprotein that is subsequently proteolytically cleaved into a series of intermediate and mature proteins, and often this ORF is subdivided into three regions namely P1, P2, and P3. Just upstream of the P1 region exists the coding sequence for the Lpro, which is preceded by two different AUG start codons that results in two different isoforms of the enzyme (Lab and Lb). The P1 region encodes four structural proteins: VP4 (1A), VP2 (1B), VP3 (1C), and VP1 (1D), while P2 and P3 encode nine nonstructural proteins: 2Apro, 2B, 2C, 3A, 3B1, 3B2, 3B3, 3C protease (3Cpro), and 3D RNA-dependent RNA polymerase (3Dpol).

The Lpro of FMDV cleaves a variety of host factors in the cell cytoplasm, which prime the host cell environment for virus replication and help to subvert the host innate immune response. As such, Lpro is considered as a significant virulence factor. In an effort to produce a useful vaccine platform, “leader less” FMDV constructs were characterized, showing considerable attenuation relative to parental virus. In a separate study, the Lpro of FMDV was functionally exchanged with BRBV Lpro to develop an FMDV–BRBV chimera vaccine.16 This vaccine showed promise in controlling FMDV infection under controlled experimental conditions. However, we speculated that the observed similarity between FMDV and BRV has another face to reckon with. Organisms containing similar (homologous) sequences may undergo recombination. In fact, multiple in-depth examinations of publicly available Picornavirus sequences have revealed evidence of significant recombination events among the +ssRNA genomes of Picornaviridae members.17–21

A large percentage of the published in silico experiments were conducted using viruses belonging to Enterovirus, Aphthovirus, and Teschovirus genera, which are the most highly represented in the public sequence databases. Recombination events and potential “recombination breakpoint hotspots”17 in the FMDV genome are reportedly confined to regions encoding the non-structural proteins and sequences flanking the genes for the capsid proteins.17,18 The delineated patterns are reportedly mirrored among other Picornavirus genera.18,20,21 Nucleotide segments encoding structural proteins appear to be largely free of recombination events, and such “recombination immunity” has been hypothesized to be the result of selection against potentially deleterious effects to genetic fitness. Indeed, the 1A–1D structural genes of Aphthoviruses are flanked by “recombination breakpoint hotspots” at the leader–1A and 1D–2A junctions.

A potential recombination between BRV and FMDV has never been observed in the past, thereby making BRVs an attractive target to explore its possible exploitation in designing chimeric vaccine candidates. Therefore, we sought to address the probability of recombination between the FMDV and BRV genomes and, in particular, BRAV RS3X. To this end, we performed several in silico comparisons of the nucleotide sequences of each virus. First, we examined the full genomic sequences of BRAV RS3X (also known as BRV type 1), FMDV A24 Cruzeiro, BRBV EC 11, and ERAV; and then different BRVs in each of the two species. Subsequently, we performed recombination analyses with software packages recombination detection program (RDP) version 4.65.22–26 Finally, positive selection pressure was first determined in the Aphthovirus genus involving BRAV RS3X and then within BRV species to infer the evolutionary ancestry of BRAV RS3X in relation to Aphthoviruses and BRVs.

This study details new and important genomic information regarding Aphthoviruses in general and specifically, we report a new member of the bovine rhinoviruses, which is BRAV. Several bioinformatics tools employed in this study help to elucidate the molecular diversity of Aphthoviruses, and help to distinguish BRAV from other members of the virus lineage. The knowledge gleaned herein and the applications of these bioinformatics tools will assist other researchers to investigate novel viral pathogens.

Materials and Methods

RNA isolation, cDNA synthesis, and nucleic acid sequencing

Viral RNA was extracted from a field sample of BRAV RS3X (Ide and Darbyshire, 1969) using an RNeasy mini kit (Qiagen NV). First strand cDNA synthesis was performed using SuperScript® III First-Strand Synthesis System for reverse transcription-polymerase chain reaction (RT-PCR; Thermo Fisher Scientific) with an oligo-dT primer and approximately 1 μg of total RNA as the template. Different fragments of RS3X cDNA from the 5′ UTR end to the 3′ UTR end were amplified in a model PTC-200 thermal cycler (MJ Research Inc.) and using the Phusion High Fidelity PCR kit (Thermo Fisher Scientific) with specific primer pairs complementary to BRAV. The PCR reaction conditions used were according to the manufacturer's recommendations for amplification products between 5 kilobase (kb) pairs and 9 kb pairs. The 5′ end of the genome between the putative poly(C) and VP3 coding regions was isolated by using a commercially available rapid amplification of cDNA ends (RACE) kit (SMART™ RACE; Clontech Laboratories Inc.). A 3′ anchored cDNA was synthesized using the complimentary gene specific primer (VP3-R1: 5′ TTGGGCCTCACTCAGAGTGGTGGGGGGAT-3′) and by following the manufacturer's protocol. By using the Advantage® 2-PCR enzyme system, 5′ RACE reactions were carried out using anchor (universal primer mix) and gene-specific primers (VP3-R2: 5′-TGGGTCCGCGGTGATGGGACTAGTGGTGC-3′) according to the manufacturer’s recommendations for amplification of products ranging from 1 kb to 5 kb. PCR products were purified using a PCR-purification kit (Qiagen NV), and the integrity of amplicons was confirmed on a 1% agarose gel. Sequences corresponding to the ends of the purified amplicons were obtained by direct sequencing with specific primers designed from previously determined partial BRAV (SD-1) sequence. Subsequent sequence data were determined by using a “primer walking” strategy in which primers for sequencing were designed based on ongoing sequence determination. All sequencing reactions were carried out using the Big Dye Terminator cycle sequencing kit (Thermo Fisher Scientific) and analyzed on a PRISM 3730xl automated DNA sequencer (Thermo Fisher Scientific).

Sequence assembly and genome annotation

Nucleotide sequences were assembled and analyzed with Sequencher (Gene Codes Corporation). The sequence reported in this work has been deposited in the GenBank database under accession number KT948520. All other GenBank accession numbers are indicated either in figures or figure legends.

Analysis of nucleic acid and amino acid sequences in relation to related Aphthoviruses

Molecular phylogenetic analysis by neighbor-joining method

The evolutionary history of BRAV RS3X was inferred by building P1 and 3Dpol coding region phylogeny in Picornaviridae using the neighbor-joining method based on the Jones–Taylor–Thornton (JTT) matrix-based model for amino acid substitution27 included in MEGA 6 software package.28 The bootstrap consensus tree inferred from 1,000 replicates is taken to represent the evolutionary history of the taxa analyzed.29 Positions containing gaps and missing data were eliminated (pairwise deletions only). There were a total of 1247 sites for the P1 region and 552 sites for 3Dpol in the final dataset.

Selection pressure analysis

Selection pressure was evaluated by determining the natural selection mechanisms acting on the codons of the ORFs of Aphthoviruses and BRVs. These mechanisms were determined using hypothesis testing using phylogenies package under the Datamonkey web-server (www.datamonkey.org).30 The dN/dS ratios (v) were calculated using three different codon-based maximum likelihood approaches: the single-likelihood ancestor (SLAC), fixed-effects likelihood (FEL), and the internal branch fixed-effects likelihood (IFEL).31,32 The mixed effects model of evolution (MEME) method, a branch-site model, was also employed for studying the selection pressure. This method is a generalization of FEL, which models variable (v) across lineages at an individual site that could to detect smaller proportions of evolving branches subject to positive selection that would otherwise be detected as “negatively” selected by FEL.33 For all the methods employed for the ORF datasets, the HKY 95 model was used as a nucleotide substitution bias model. Trees were inferred by the neighbor-joining method and significance levels were set at P < 0.05 or Bayes factor 0.50.

Recombination analysis

Possible recombination events among different Aphthoviruses and BRVs were assessed separately using RDP v.4.65.34 In default mode, RDP, GENECONV, CHIMAERA, MAXCHI, BOOTSCAN, PHYLPRO, LARD, SISCAN, and 3SEQ algorithms were utilized to detect potential recombination events between the input sequences.34

Structural modeling

Homology models of BRAV RS3X 2B and FMDV A24 Cruzeiro 2B were built using the hepatitis C virus structure (protein data bank (PDB): 2MTS). VP1 of FMDV A24 Cruzeiro and BRAV RS3X were modeled using the FMDV type O capsid (PDB: 1FOD Chain 1) as a template. All of the homology models were prepared on SWISSMODEL workspace.35–39 The stereochemical quality of the models was further validated with PROCHECK.40 Structures were rendered using USCF-Chimera 1.10.41

Results

Genome and genome encoded proteins

We sequenced the near full-length genome of BRAV RS3X (7,267 nucleotides) from a putative poly(C) tract at the 5′ UTR end to the poly(A) tail at the 3′ end.

5′ UTR

Starting with the 5′ UTR poly(C) of BRAV, RS3X shares 64.12%, 54.00%, and 41.03% nucleotide sequence identity to FMDV (A24 Cruzeiro), BRBV, and ERAV, respectively.

Lpro

ORF scan analysis suggests that Lpro of BRAV RS3X is smaller in comparison with the Lab form observed in other Aphthoviruses including FMDV, BRBV, and BRAV isolates BSRI 4 and 140032-1. As shown in Figure 1, the Lpro active site (highlighted in yellow color)42 and translation initiation factor-binding sites (highlighted in cyan color) are conserved among BRAVs. Importantly, in BRAV RS3X, these sites are identical to other members of the Aphthovirus genus (Figs. 1A and B).43–45

Figure 1

(A) Alignment of polyprotein of Aphthoviruses: FMDV A24 Cruzeiro strain was used as a template for alignment and feature annotation and marked as FMDV A24 Cruz. BRAV RS3X (GeneBank: KT948520), BRBV EC-11 (GeneBank: EU236594) and ERAV (GeneBank: DQ272578) were aligned to FMDV A24 Cruz (GeneBank: AY 593768). Marks the start of a protein sequence (N-terminus of), and the text following it identifies the protein. Yellow highlights suggest functionally critical (active site residue). (·) denotes a conserved residue in comparison to FMDV A24 Cruz. (−) denotes a gap introduced in the alignment. Cyan highlights indicate functionally important residues of a protein other than active site residues. Green highlights denote the RGD cell surface receptor. (B) Alignment of polyprotein of BRVs: Color code of highlighting and feature marking is the same as (A). Text at the beginning of each line is the sequence identifier for the given virus strain. For example, BRAV RS3X indicates BRAV RS3X sequence. The GeneBank accession numbers of the sequences used in the alignment are KT948520, KP236128, KP236129, JN936206, KP264974, EU236594, KP236130, KP264975, and KP264980 for BRAV RS 3X, BRAV Sd-1, BRAV 140032-1, BRAV H-1, BRAV BSRI4, BRBV EC-11, BRBV 140032-2, BRAV BSRI3, and BRAV BSRI1, respectively.

P1 region

The structural protein VP4, the first encoded protein from the P1 region, shows the highest conservation among Aphthoviruses as well as BRVs. However, BRVs encode a longer VP4 protein than other Aphthoviruses. The N-terminus of VP2 is remarkably conserved among all Aphthoviruses. The rest of the VP2 protein and the entire VP3 and VP1 of BRVs and ERAV are distinctively different from that of FMDV. Additionally, the characteristic RGD tripeptide, which serves as the cell surface-binding site in the VP1 G–H loop of FMDV,46–48 is lacking in all other Aphthoviruses, suggesting an alternative mechanism of their interaction/attachment to the cell surface (Figs. 1A and B). Considering the absence of the RGD tripeptide in BRAV RS3X VP1, structural models of FMDV A24 Cruzeiro and BRAV RS3X were prepared and compared. Structural models of the two proteins when superimposed on each other show that the G–H loop region is completely disordered in the BRAV RS3X VP1 molecule, which reinforces the speculation that a different virus–receptor interaction occurs on the cell surface (Fig. 2Bi).

Figure 2

(A) Prediction of the transmembrane domain in: (i) FMDV A24 Cruz 2B and (ii) BRAV RS3X 2B. Color codes for different features are explained adjacent to each image. (B) Modeled structures of VP1 structural proteins of (i) FMDV A24 Cruz and BRAV RS3X are colored cyan, the G–H loop is marked with a yellow oval, and the RGD tripeptide is marked with cyan text. N- and C-terminus of the proteins in the structure are marked. (ii) The prediction of the 2B viroporin of: FMDV A24 Cruz (purple) and BRAV RS3X (cyan) were superimposed on HCV p7. N- and C-terminus of the proteins in the structure are marked. Structure-based sequence alignments are positioned to the left of each image. All structures are rendered in cartoon representation using UCSF-Chimera ver 1.10

Non-structural proteins

2A containing a C-terminal ribosome-skipping motif NPG↓P encoded by all Aphthoviruses (and many other Picornaviruses) is similar in size. The N-termini of 2A show greater variation, suggesting a lesser functional significance for this region.49,50

The most obvious differences between BRVs (including BRAV RS3X reported here) and FMDV (Figs. 1A and B) are the sizes of the 2B and 3A proteins as well as the number of 3B (VPg) peptides encoded. The N-terminus of the 2B “viroporin”51–53 in all Aphthoviruses, but FMDV bears a significant deletion toward the N-terminus. TMMH (membrane protein topology prediction method) web server analysis revealed that 2B of BRAV RS3X lacks a trans-membrane region that is observed in FMDV A24 Cruzeiro (Figs. 2Ai and ii). However, when we built a homology model of BRBV RS3X 2B using hepatitis C virus (HCV) p7 viroporin (PDB: 2MTS), the viroporin domain appeared to be conserved in BRAV RS3X and was very similar to HCV p7. In fact, the pore lining histidine residue is only present in BRAV RS3X, and FMDV 2B lacks the pore lining α-helix altogether, suggesting a different mechanism for FMDV 2B interaction with intracellular membranes (Fig. 2Bi).

The non-structural protein 3A contains two hydrophobic residues (L38 and 41 in FMDV) that are supposedly important for its attachment to intracellular membranes. Although the hydrophobicity of the amino acid residues corresponding to residues 38 and 41 of FMDV 3A are conserved among all the Aphthoviruses, the smaller 3A of BRAV RS3X and ERAV followed by BRBV indicates a divergence in the membrane association of 3A of these viruses from FMDV.54 The functional significance of smaller 2B and 3A proteins remains undefined and invites further investigation.

Furthermore, the 2C-ATPase is conserved in BRAV RS3X as well as other Aphthoviruses when compared to FMDV. Importantly, the amino acid residues corresponding to residues 116, 160, and 207 of FMDV 2C that are critical to 2C activity are conserved in BRV RS3X and other Aphthoviruses compared here (Figs. 1A and B).

FMDV is unique among Picornaviruses with respect to encoding three non-identical copies of the 3B (VPg) protein in its P3 region, although possible examples of two copies of VPg are speculated for members of the Aquamavirus, Mosavirus, and possibly Passerivirus genera (N.J.K., unpublished observations). Unlike FMDV, BRAV RS3X and all other Aphthoviruses encode only one copy of the nonstructural protein 3B that closely resembles the third non-identical copy of FMDV 3B (3B3). Importantly, the tyrosine residue at position 3, a target for uridylation and a substrate for virus genome replication, is conserved in this single copy gene.55

Finally, as expected, both the 3Cpro and 3Dpol (highlighted in yellow color) are most conserved and catalytically important residues of the two proteins (highlighted in yellow color) are absolutely conserved in all the viruses compared herein.56,57

Phylogeny

We constructed neighbor-joining phylogenetic trees of the structural protein region P1 as well as non-structural protein 3Dpol, the former being most diverse and latter being most conserved. The BRAV RS3X sequence reported in this study and all of the classified Picornaviruses (54 species in 31 genera, as shown in Fig. 3) were utilized to illustrate the evolutionary space of BRAV RS3X in Picornaviridae.11 The amino acid sequences of the P1 and 3Dpol regions were first aligned with the Muscle algorithm implemented in MEGA 6, and later were manually curated to ensure the accuracy of the alignment. Finally, a phylogeny tree was constructed as specified in the Materials and methods section. The resulting P1 phylogenetic trees (Fig. 3A) show clear divergence between the BRAVs separated by distinct nodes. The RS3X strain and isolate BSRI4 cluster together, while isolate SD1 forms another distinct node, and isolate H1 and strain 140032-1 form a third distinct cluster. Other Aphthoviruses including FMDV strain A24 Cruzeiro, BRBV, and ERAV each form separate nodes characteristic of their genera. However, as would be expected, all of the Aphthoviruses group together, thus proving their close evolutionary relationship within the genus Aphthovirus.

Figure 3

Neighbor-joining phylogeny trees of the P1 structural protein region (A) and 3Dpol non-structural protein (B) of family Picornaviridae including BRAV RS3X. Neighbor-joining phylogeny trees involving BRAV RS3X and all known members of Picornaviridae were constructed with the JTT model for amino acid substitution (pairwise deletions only) for (A) P1 structural proteins and (B) 3Dpol non-structural protein. As shown, 1,000 bootstraps with only those of 70% and above are displayed. Scale of resolution is indicated at the bottom of each figure.

Note: “ ” indicate proposed new species and genera in Picornaviridae.

The NJ–JTT phylogeny tree of the 3Dpol region showed similar clustering of Aphthoviruses (Fig. 3B). BRAV RS3X and the isolates BSRI4, SD1, H1, and 140032-1 form a cluster distinct from FMDV A24 Cruzeiro, BRBV, and ERAV. Finally, the order of relatedness exactly matched with that of the P1 region phylogeny.

Selection pressure analysis

Positive selection is a major mechanism of RNA virus evolution. Therefore, we first analyzed the selection pressure between Aphthoviruses that are the closest relatives to BRAV by including the BRAV RS3X strain, FMDV A24 Cruzeiro, BRBV, and ERAV by various methods included in the Datamonkey web-server. The results shown in Table 1 show that the SLAC method did not detect any codon site subject to positive selection between these viruses. However, the FEL method detected codon sites 489 and 832 in the structural protein P1 region and codon site 1958 in the 3Dpol region under positive selection pressure. The IFEL method that computes site-wise selection only on the internal branches of the phylogeny tree detected 14 sites subject to positive evolutionary selection, with ten sites in the 3Dpol region and four in the structural protein P1 region. Together, these data suggest that the 3Dpol region has the higher prevalence of evolutionary selection. SLAC and FEL methods detected 363 and 425 sites under negative evolutionary pressure, whereas IFEL unexpectedly found only 3 sites under negative selection. It appears that the selection pressure detected by IFEL, by virtue of working on the internal branches of the phylogeny tree, detects mostly the selection events within species and hence BRAV RS3X and BRBV, being more closely related, would yield such unexpected results.

Table 1

Distribution of positively or negatively selected sites among Aphthoviruses.

ANALYSIS	NEGATIVE SITES				POSITIVE SITES
	NO.	CODON SITE	AVG. dN/dS	P-VALUE	NO.	CODON SITE	AVG. dN/dS	P-VALUE
SLAC	363	N.A.	7.65015e-05	0.004–0.05	0	N.A.	N.A.	$0.05
FEL	425	N.A.	0.001946	0.002–0.05	3	489, 832, 1958	Infinite	0.038–0.045
IFEL	3	1119, 2001, 2439	0.0003	0.002–0.05	14	393, 431, 465, 993, 1123, 1471, 1742, 1880, 1901, 1914, 1920, 2055, 2059, 2331	Infinite	0.008–0.05
MEME	N.A.	N.A.	N.A.	N.A.	6	49, 735, 884, 2054, 2331, 2394	N.A.	0.015–0.048

Analysis of selection pressures acting among BRVs by the SLAC method did not detect any positive selection. The FEL method detected one codon site, 680 in the structural protein region under positive selection. However, the IFEL method detected six codon sites subject to positive selection. Four of these sites (614, 617, 680, and 811) are in the structural proteins, one (139) in Lpro, and one (1806) in 3Cpro. All three methods identified numerous sites under negative selection pressure in BRVs with SLAC, FEL, and IFEL detecting 511, 1143, and 843 sites subject to negative selection, respectively. The MEME method, which represents an advancement to the FEL method and detects positive selection sites that would be detected negative in the FEL method, detected 6 codon sites under positive selection among Aphthoviruses: three in the 3Dpol region at positions 2054, 2331, and 2394; one at position 49 in the Lpro region; and sites 735 and 884 in the P1 structural region (Table 1). Among BRVs, the MEME method detected 22 sites subject to positive purifying selection and distributed evenly along the ORF, clearly suggesting the existence of strong positive selection within BRV species (Table 2).

Table 2

Distribution of positively or negatively selected sites among Bovine rhinitis viruses.

ANALYSIS	NEGATIVE SITES				POSITIVE SITES
	NO.	AVG. dN/dS	P-VALUE	NO.	CODON SITE	AVG. dN/dS	P-VALUE
SLAC	511	0.116349	0.005–0.05	0	N.A.	N.A.	≥0.05
FEL	1143	0.016	≥0.05	1	680	Infinite	0.034
IFEL	843	0.00539	0.002–0.05	6	139, 614, 617, 680, 811, 1806	Infinite	0.008–0.05
MEME	N.A.	N.A.	N.A.	22	3, 5, 95, 448, 617, 630, 667, 684, 733, 789, 798, 811, 961, 985, 1008, 1010, 1249, 1412, 1479, 1490, 1745, 1901	N.A.	0.001–0.05

Determination of potential recombination events

During recombination, two molecules of DNA or RNA that carry matching sites (homologous sequences) exchange their segments to yield novel combinations. In fact, recombination is considered to contribute significantly to RNA virus evolution. When analyzing a group of nucleotide sequences for the probability of recombination events, the most common first step is to search for so-called “recombination breakpoints” in the existing sequences. Multiple algorithms have been designed to dissect an alignment of several nucleotide sequences for local regions that exhibit the hallmarks of a recombination breakpoint. The RDP software package uses multiple algorithms, such as RDP, GENECONV, CHIMAERA, MAXCHI, BOOTSCAN, PHYLPRO, LARD, SISCAN, and 3SEQ. In its default mode, RDP calculates even the least possible event. However, a recombination event should only be considered significant if it is evidenced by multiple methods. In this study, we considered an event to be significant only when evidence was provided by four or more methods. In this way, we could take into account the results produced by four out of nine methods employed. Such a strategy has been applied to interpret the recombination events detected using RDP.58 Table 3 shows the results from an analysis of BRAV RS3X and three other viruses: FMDV A24 Cruzeiro, BRBV, and ERAV. Due to several breaks produced in the ORF due to inclusion of both the 5′ and 3′ UTR, these regions were excluded from the alignment and only single ORFs from the Lpro to 3Dpol regions of these viruses were aligned and included for recombination detection. RDP that recognizes even the remotest possible event, recognized six potential recombination hot spots in the ORF of these viruses, and BRAV RS3X was involved in five of them. However, none of these were confirmed by more than four methods set up as cut-off, and hence we concluded that the possibility of inter-species recombination between different Aphthoviruses analyzed herein is negligible.

Table 3

Detection of potential recombination events between BRAV RS3X and closely related Aphthoviruses using ORF only.

RENO.	BREAKPOINT POSITIONS IN RECOMBINANT SEQUENCE							DETECTION METHOD
	BEGIN	END	RECOMBINANT SEQ (S)	MINOR PARENT SEQ (S)	MAJOR PARENT SEQ (S)	RDP	GENE-CONV	BOOTSCAN	MAXCHI	CHIMAERA	SISSCAN	PHYLPRO	LARD	3SEQ
1	3127	3371	BRBV	ERAV	FMDV A24	2.07E-03	NS	NS	NS	NS	NS	NS	NS	NS
2	6360	6564	BRAV RS3X	Unknown FMDV A24	BRBV EC11	2.47E-03	NS	1.76E-02	NS	4.17E-02	NS	NS	NS	NS
3	5750*	5824	FMDV A	ERAV	BRAV RS3X	1.81E-02	NS	NS	NS	NS	NS	NS	NS	NS
4	2238	2322	FMDV A24	ERAV	BRAV RS3X	2.21E-02	NS	NS	NS	NS	NS	NS	NS	NS
5	965	1122	BRAV RS3X	FMDV A24	Unknown BRBV EC11	2.67E-02	NS	NS	NS	NS	NS	NS	NS	NS
6	6640	54	BRAV RS3X	FMDV A24	Unknown BRBV EC11	0.044604	NS	NS	NS	NS	NS	NS	NS	NS

Notes:

The actual breakpoint position is undetermined (it was most likely overprinted by a subsequent recombination event). Minor Parent = Parent contributing the smaller fraction of sequence. Major Parent = Parent contributing the larger fraction of sequence. Unknown = Only one parent and a recombinant need be in the alignment for a recombination event to be detectable. The sequence listed as unknown was used to infer the existence of a missing parental sequence. NS = No significant. P-value was recorded for this recombination event using this method.

We then analyzed different BRVs for potential recombination events using RDP. As shown in Table 4, in contrast with inter-species recombination being insignificant within species, BRVs carry multiple recombination hot spots distributed throughout the ORF. No less than six of these were reproducible by multiple methods highlighting that these points are real recombination break points.

Table 4

Detection of potential recombination events among bovine rhinitis viruses using ORF only.

RENO.	BREAKPOINT POSITIONS IN RECOMBINANT SEQUENCE					DETECTION METHOD
	BEGIN	END	RECOMBINANT SEQ (S)	MINOR PARENT SEQ (S)	MAJOR PARENT SEQ (S)	RDP	GENECONV	BOOTSCAN	MAXCHI	CHIMAERA	SISSCAN	PHYLPRO	LARD	3SEQ
1	270	2930	BRAV BSRI4	BRAV 140032-1	BRAV SD-1	1.19E-48	6.66E-39	6.22E-58	4.26E-38	8.94E-33	9.16E-51	NS	NS	2.32E-74
2	6422	2963	BRBV BSRI1	BRBV 140032-2	BRBV EC-11	4.09E-28	1.36E-38	1.63E-35	1.15E-09	2.19E-15	2.86E-32	NS	NS	1.35E-11
3	2962	704	BRAV RS3X	BRAV_Sd-1	BRAV 140032-1 BRAV H-1	5.11E-15	NS	6.50E-09	4.18E-14	9.34E-17	1.33E-36	NS	NS	NS
4	542	6289	BRBV 140032-2	Unknown BRBV BSRI3	BRBV EC-11	1.32E-14	NS	3.62E-14	3.30E-11	1.48E-05	9.78E-05	NS	NS	NS
5	3415	2963*	BRBV BSRI1	BRBV EC-11	BRBV 140032-2	3.68E-09	3.23E-04	NS	3.60E-09	NS	3.33E-06	NS	NS	3.30E-05
6	3069	2899	BRBV BSRI3	Unknown BRAV SD-1 Unknown BRAV 140032-1	BRBV_140032-2	3.16E-07	NS	NS	2.18E-07	1.19E-03	NS	NS	NS	NS
7	138	6407	BRBV BSRI1	BRBV_BSRI3	BRBV_EC-11	3.24E-02	NS	NS	NS	NS	NS	NS	NS	NS

Notes:

The actual breakpoint position is undetermined (it was most likely overprinted by a subsequent recombination event). Minor Parent = Parent contributing the smaller fraction of sequence. Major Parent = Parent contributing the larger fraction of sequence. Unknown = Only one parent and a recombinant need be in the alignment for a recombination event to be detectable. The sequence listed as unknown was used to infer the existence of a missing parental sequence. NS = No significant. P-value was recorded for this recombination event using this method. Sequence id is in the original sequence not in the alignment or in relation to any parent.

The findings of these algorithms (Tables 3 and 4) for the assessment of recombination potential are consistent with other published studies regarding the potential of Picornavirus genomes to undergo recombination events.34 We concluded based on these findings and what was previously described in the existing literature that the likelihood of detecting recombination events between FMDV and BRV could be described as highly remote. We also inferred that given the absence of prior exchanges between the genomes with two different sample pools, the chance of future exchanges would be novel and highly unlikely.

Discussion

In this study, a novel strain RS3X of BRAV was sequenced to its near full-length. Annotation of features and comparison of sequence alignment of BRAV RS3X to related Aphthoviruses and BRVs revealed several unique features in its genome. As one would expect, the BRAV RS3X sequence exhibits conserved features of the BRAV species. Major differences between the prototypic Aphthovirus FMDV and BRAV RS3X were observed in the architecture of capsid and non-structural proteins 2B and 3A. Only the capsid protein VP1 of FMDV displays a cellular integrin-binding arginine-glycine-aspertate (RGD) tri-peptide motif in its G–H loop.46–48 The remaining three Aphthoviruses (BRAV, BRBV, and ERAV) lack this essential cellular receptor-binding site as revealed by structural modeling and analysis of sequence alignment, suggesting a different mechanism of cell attachment for these viruses. In fact, the ERAV crystal structure revealed that it binds the cell surface via sialic acid.59 The 2B viroporin of BRAV RS3X is significantly smaller than FMDV, and lack of a C-terminal trans-membrane helix of the former indicates a diversion in its topology or function from FMDV. Surprisingly, the modeled structure of the 2B protein from FMDV and BRBV suggest that the latter preserves the membrane pore-lining histidine residue, whereas FMDV 2B lacks this characteristic feature of a viroporin molecule.60 Conservation of active sites and other functionally critical residues in Lpro, 2Apro, 3Cpro, 3Dpol, and membrane anchorage residue in the 3A protein of BRAV RS3X when compared with Aphthoviruses (FMDV, BRBV, and ERAV) and within BRV species reinforces the functional retention of the molecular biology of Aphthoviruses in BRAV RS3X.

Neighbor-joining phylogenetic trees constructed with the most diverse P1 region and most conserved 3Dpol region parallel each other, confirming the accuracy of phylogeny tree inference. In both analyses, RS3X clustered with BRAV isolates, confirming its classification in that species. The close proximity of BRAV RS3X with isolates BSRI4 and SD-1 suggests a closer evolutionary relationship between these isolates, thus confirming their classification in the BRAV-1 serotype (formerly known as bovine rhinovirus type 1) to which it is antigenically related. Isolates H-1 and 140032-1 form a distinct cluster which has been serologically classified as a second serotype, BRAV-2 (formerly known as bovine rhinovirus type 3).

Potential recombination breakpoints in the genomes were detected with sites identified flanking the P1 region containing the structural protein-coding sequences and others further downstream in a region containing the nonstructural protein-coding sequences. The evaluation of potential recombination breakpoints was performed using distinct sets of algorithms dependent upon which software package was employed. The similarity in results obtained using both algorithmic sets strengthens the final conclusion that recombination events are unlikely between the FMDV and BRV genomes within the regions encoding the structural proteins. Based on examination of the recombination analysis data of Aphthoviruses and BRVs, we conclude that the inter-species recombination events involving FMDV, BRAV, BRBV, and ERAV do not seem likely. However, within each of the two BRV species, the occurrence of multiple recombination breakpoints confirms an underlying phenomenon of homologous recombination-mediated generation of diversity in Picornaviridae.61–63 Despite being the earliest discovered animal virus, there is no report of FMDV recombining with any other virus. On the other hand, there is a significant probability for intra-species recombination within each of the two BRV species. In fact, multiple studies have shown intra-species recombination in the prototypic Aphtho virus FMDV.58,64–67

Selection pressure analysis is a significant methodology for depicting a common lineage for rapidly evolving organisms, and has been employed extensively in inferring the evolutionary information of viral populations.32,58,68–71 Given that multiple analytical methods were employed in this study, there is a high degree of confidence added to the interpretation that positive selection is not operating between BRV and other Aphthoviruses. The observed higher number of positive selection breakpoints among BRVs is not a surprise and reinforces their common ancestry. As expected, negligible evidence of single-point positive selection sites proves a parallel but distinct lineage of the different Aphthoviruses included in this study.

In conclusion, the data from this study provide valuable information on Aphthoviruses and, more precisely, BRAV to serve as the genetic basis for future studies. Detailed knowledge of the evolution and divergence of Aphthoviruses at the molecular level could aid in the design of BRV-based molecular diagnostic tools and new bio-therapeutics.