Lineage diversification, homo- and heterologous reassortment and recombination shape the evolution of chicken orthoreoviruses

The near complete genome sequences of ten field avian orthoreovirus (ARV) strains collected from young chicken between 2002 and 2011 in Hungary have been determined in order to explore the genetic diversity and evolutionary mechanisms affecting ARVs in this region. Sequence analyses and phylogenetic calculations revealed similar geographic distribution and distinct evolution in case of eight studied strains that were closely related to the recently described Hungarian strain T1781. The remaining two strains showed the highest similarity with the US origin AVS-B. The topology of the phylogenetic trees of certain segments was affected by several potential homologous reassortment events between strains of Hungarian, Chinese and US origin. Analyzing the μB gene a possible heterologous reassortment event was identified in three Hungarian strains. Recombination events were detected in as much as a dozen cases implying that beside point mutations and reassorment this mechanism also plays an important role in the diversification of ARVs. All these mechanisms in concert may explain the reduced effectiveness of immunization using commercial vaccine strains.

Although the number of complete genome sequences available in GenBank is still limited, sequence data collected in the last few years suggest that the genetic material of ARV strains circulating in domestic poultry is continuously changing due to the well-known evolutionary mechanisms of reoviruses, resulting in genetically and antigenically heterogeneous strains. As observed in other RNA viruses, point mutation which is the primary mechanism of genetic and antigenic drift is relatively common. Antigenic shift, another evolutionary strategy, occurs when cognate genome segments reassort, and this mechanism is thought to be the driving force of the emergence of novel strains causing disease in poultry 1 . As a consequence of diverse evolutionary mechanisms, the commercially available ARV vaccines are not able to provide adequate protection against newly emerging field strains leading to outbreaks of disease in vaccinated poultry flocks [11][12][13] .
In this study we determined the near complete genome sequences of ten ARV field strains collected from young chicken with different clinical signs over a 10-year period in Hungary (Table 1) in order to gain more genetic data about the circulating strains, and to determine their relatedness to contemporary, predominantly Asian and US origin ARV strains, the genomic sequences of which are accessible in GenBank. Identifying intra-segmental recombination and heterologous reassortment events permits new insight into the evolutionary mechanisms of ARV strains commonly detected in chicken. European ARV strains have been underrepresented in previous studies thus this survey provides important new information about the geographic distribution and evolutionary origin of chicken ARVs.

Results and Discussion
Genomic organization and coding potential. The near complete genome sequences of ten, randomly selected Hungarian chicken ARV isolates had been determined by a modified sequence-independent, single-primer amplification method coupled with next-generation sequencing. The genomic organization of these ARV strains was similar and corresponded with that of the previously described reference ARV strains, AVS-B and S1133. With the exception of the S1, all genomic segments were found to encode a single protein and the homologues of the following orthoreoviral proteins were identified by open reading frame prediction and sequence comparison: λ A, λ B, λ C, μ A, μ B, μ NS, σ A, σ B, σ C, σ NS, p10 and p17. Where available, the 5′ and 3′ terminal sequences (GCUUUU[U/C] and UCAUC) were conserved among all segments as observed in other ARVs described so far 5,14 . Sequence homology and phylogenetic relationship among strains. Sequence comparison of the Hungarian strains revealed 71.1-100/83.0-100% nucleotide sequence identity and aa sequence similarity values with the exception of the μ B and σ C genes where higher genetic variability was observed, and these values ranged between 62.7-99.8% nt and 71.3-99.7% aa for μ B and 52.9-100% nt and 46.2-100% aa for σ C, respectively ( Fig. 1, Table 2).
Eight out of ten Hungarian strains (284-V-06, 924-Bi-05, 3211-V-02, 3457-M-11, 4599-V-04, 16821-M-06, 17203-M-06, 17227-M-10) showed the highest homology with T1781, another Hungarian ARV strain we characterized recently 8 . The common evolutionary origin of these strains is clearly reflected by the results of the analysis performed with the nt and aa sequences of the λ A gene. When comparing these strains with T1781, the nt identity and aa similarity values were high (falling between 91.1-95.1% and 98.8-99.3%, respectively). Lower values were obtained in the calculations performed with the sequences of US origin reference strains, S1133 and AVS-B, isolated in 1971 and 2006, respectively. In the case of these strains the nt identity values ranged between 83.1-84.1% and the aa similarity values were between 82.5-83.3% and 97-97.6%, respectively. Sequence analyses revealed that the remaining two strains 875-Bi-05 and 878-Bi-05 were most closely related to each other showing high, 98.8-100% nt sequence identity and 99.7-100% aa sequence similarity, and with the exception of the σ C   Phylogenetic analyses of the 10 genomic segments resulted in different tree topologies suggesting that certain segments may have been constituted by distinct evolutionary mechanisms; several hypothetical reassortment events between strains originating from different geographic regions were observed indicating that this phenomenon is relatively common in the Avian orthoreovirus species (Fig. 2). With the exception of the μ B gene in all phylogenies host specific evolution of gallinaceous and waterfowl origin ARVs was evident 15 . The geographic origin of different ARV strains based on their genetic relatedness could also be observed especially in the case of the more conserved core proteins. Our calculations suggested that numerous reassortment events occurred  recently between Hungarian, US, and Asian strains affecting mainly the outer capsid and non-structural protein coding genes, therefore the geographic origin of strains in these particular cases was less apparent. Phylogenetic analysis of the μ B and σ C outer capsid proteins, responsible for virus entry and transcriptase activation, and antigenic properties, respectively, revealed greater genetic distances between the different ARV strains. In detail, most of the studied Hungarian strains  Table 3). Calculated nt identity values between the NBV and TVAV, Pycno-1 and SSRV, not classified into the same orthoreovirus species, were even higher 65.5-66.1%. At the aa level the similarity values calculated with S1133, AVS-B and T1781 ranged between 71.7-73.4%, while these were lower, 64.4-69.1%, when calculations were performed with the ARV-like viruses  Table 3), suggesting that some sort of adaptation process to the gallinaceous host and/or fine-tuning among interacting virion components of true chicken origin ARVs has already begun. In some positions these strains had a unique aa composition characteristic only to these three Hungarian strains when compared with other known members of the genus Orthoreovirus (data not shown). These Hungarian reassortant strains were collected in 2005, 2006 and 2011 in Hungary and strains possessing a closely related gene have never been detected before and after in any other countries. The first seven nt of the 5′ terminal sequences (GCUUUUU) of μ B of the three strains showed high level of conservation with ARVs, implicating that this conserved terminal sequence may have been a prerequisite of the acquisition of this gene by reassortment. Apart from this conserved motif the 5′ untranslated region (UTR) sequences showed considerable variability and were 4 nt shorter (25 nt) when compared with the μ B encoding genome segment of other ARV strains supporting the theory of the heterologous reassortment. However, alignment of the 5′ UTR containing orthoreoviruses belonging to other genera and unclassified strains revealed the highest similarity with ARV-like viruses, especially with NBV and Pulau virus (data not shown). Three distinct lineages of orthoreoviruses have been described in previous studies, suggesting a common evolutionary history of each of the following clades: (i) the reptilian orthoreovirus strains with Baboon orthoreovirus and Broome virus, (ii) the classical mammalian orthoreoviruses, (iii) the avian orthoreoviruses along with NBV, TVAV and SSRV. It is possible that reassortment might occur between species belonging to the monophyletic group of the third clade, but the exact genetic breakpoint remains to be discovered. The genetic and phylogenetic distance is considerably smaller between representative members of the different species in the ARV-like clade, than that of the two other clades. Although it is documented that reassortment occurs only between members of a given orthoreovirus species 1 , to the best of our knowledge co-infection studies have not been systemically performed to rule out the possibility of an inter-reovirus-species reassortment event.
Intrasegmental recombination events. To uncover whether recombination involving cognate genome segments occurs in ARV, particularly the Hungarian isolates, we performed recombination analysis with the studied Hungarian strains and sequences of ARVs (chicken and turkey) available in GenBank. Analysis of the coding region of each genomic segment revealed twelve putative recombination events in the λ A (2), λ C (3), μ A (3), μ B (1), μ NS (2), and σ A (1) genes. Details of the analyses are provided in Table 4 and Fig. 4. Among the few identified cases, the Hungarian strains were involved in one possible recombination event, namely, the isolate 601 G was identified as a possible recombinant that may have arisen from a strain closely related to 16821-M-06 and another strain related to GuangxiR2. In the data set no recombination event was detected within the σ C coding gene, suggesting that serotype diversity is not derived from this mechanism.

Conclusions
Our results provide insight into the molecular characteristics of ARV strains randomly selected from a major strain collection. Genome sequencing and bioinformatics analyses of these strains provided new insight into the evolutionary mechanisms of ARV strains detected over a decade in Hungary. We conclude that lineage diversification, recombination involving cognate genome segments, reassortment among homologous strains and even among heterologous strains drive the evolution and genetic diversity of chicken orthoreovirus strains in Hungary. It remains to be studied whether similar mechanisms affect the genetic diversity of ARV strains in other parts of the world. A better understanding of the relative importance of these mechanisms will likely affect current vaccination strategies and the control and prevention of reovirus infections.

Materials and Methods
ARV isolates. Samples processed in this study originated from young chicken found dead following a course of disease with different manifestations in commercial poultry flocks between 2002 and 2011 in Hungary and submitted for post mortem evaluation to the Veterinary Diagnostic Directorate, National Food Chain Safety Office, Debrecen. Chicken ARV strains had been isolated from the collected specimens and the tissue culture supernatant of each isolate had been lyophilized for later use (Table 1). Prior to next generation sequencing, the lyophilized samples of ten randomly chosen strains had been dissolved in sterile double-distilled water. The routine diagnostic investigations were performed according to the current directives specified in law 1998./XXVIII. and 2008./XLVI., and regulation 41/1997. (V. 28.).
Whole genome sequencing. Sequencing was performed using the protocol described in detail elsewhere 16 .
In brief, RNA was extracted from diluted samples using TRI Reagent (Sigma-Aldrich, Saint Louis, MO, USA)  Computer analyses. Contigs were aligned with Sanger sequencing reads using MultAlin online software 17 and were edited in GeneDoc 18 . Phylogenetic analysis was performed with the MEGA6 program package 19 based on multiple sequence alignments generated by the TranslatorX online platform 20 ; the best-fit substitution models were selected for each gene-specific dataset based on the Bayesian information criterion. Maximum-likelihood trees were generated and tree topologies were validated by bootstrap analysis (100 replicates). Pairwise comparison of nucleotide sequences and amino acid p-distances of the Hungarian ARV strains and reference ARV strains, S1133, AVS-B and a Hungarian strain T1781 were calculated by applying the BioEdit 21 and MEGA6 program package, respectively. The accession numbers of orthoreovirus strains used in sequence and phylogenetic analyses are listed in Supplementary Table S1. The aligned sequences of the Hungarian isolates and ARV sequences available from GenBank were subjected to recombination analysis. The Recombination Detection Program (RDP) package v.3.44 22 was used for identification of recombinant sequences and breakpoints using the default parameters for the methods GENECONV, Bootscan, Chimaera, MaxChi, SiScan, 3Seq and RDP. Only those recombination events were taken into considerations that were supported by at least 4 methods to avoid misidentification using a single methodology. The best signals for recombination are associated with the lowest P-values; the highest acceptable P-value was set to 0.05. Recombination events detected with RDP v.3.44 were confirmed and visualized with SimPlot Version 3.5.1 23 . Analysis and visualization of aligned concatenated whole genome sequences was performed with mVISTA online platform 24 .