Development of the Intestinal RNA Virus Community of Healthy Broiler Chickens.

Several RNA viruses such as astrovirus, rotavirus, reovirus and parvovirus have been detected in both healthy and diseased commercial poultry flocks. The aim of this study was to characterize (a) the development of the RNA viral community in the small intestines of healthy broiler chickens from hatch through 6 weeks of age (market age) and (b) the contribution of the breeder source vs. bird age in development of the community structure. Intestinal tissue samples were harvested from breeders and their progeny, processed for viral RNA extraction and sequenced using Illumina Hiseq sequencing technology resulting in 100 bp PE reads. The results from this study indicated that the breeder source influenced the RNA viral community only at hatch but later environment i.e. bird age had the more significant effect. The most abundant RNA viral family detected at 2, 4 and 6 weeks of age was Astroviridae, which decreased in abundance with age while the abundance of Picornaviridae increased with age.


Introduction
Vertebrates are host to large communities of bacteria, fungi, archaea and viruses collectively known at the microbiome. Microbial communities, especially the bacterial constituent, are known to contribute to the healthy states of their hosts. In contrast, eukaryotic viruses that infect the host and inflict cellular damage, have commonly been studied as pathogens rather than as members of a healthy microbiota. But, as evidence emerges that there are diverse viral communities in healthy vertebrate hosts, that paradigm is shifting [1]. With new methods to study the virome [2], the picture that is beginning to emerge is that of eukaryotic viruses as part of a healthy microbiota and part of a host's genetic identity [3].
Although viruses are an essential component of the normal microbiota of the gut, there are also virulent pathogens that can cause enteric disease. Most of these viruses have RNA genomes which infect the host and cause lesions. Rotaviral diarrheal diseases, for example, are a major cause of illness among children, calves, foals, pigs, chicks and poults [4,5]. In addition to rotaviruses, other RNA viral families have been linked to enteritis including Astroviridae, Reoviridae, Coronaviridae and Parvoviridae [6][7][8][9]. Despite links to disease, these same viral families have been detected in the gastrointestinal tracts of healthy animals of the same species [10][11][12][13] suggesting that the relationship between host and virus likely includes both symbiosis and disease, and some balance between the two states.
Although eukaryotic viruses may not be in an entirely symbiotic relationship with their hosts, it is evident that not all infections are symptomatic. To investigate the dynamics of the members of the RNA virome, we examined the RNA viruses in the gastrointestinal tract of healthy chickens. As a host, chickens can be infected with several enteric virus families and are susceptible to clinical enteritis. In our previous study, we developed a method to process the whole intestine, both tissue and contents, using high-throughput next generation sequencing to characterize the RNA viral community [2]. In this study, using the established method, we followed broiler chickens from hatch to market as well as their dams to understand how RNA viral communities develop and change as the host ages and how they relate to enteric health.

Study design
Tissues were collected post mortem from poultry raised on commercial farms (source: GNP company, St. Cloud, MN). At each sample collection timepoint except for the hatch timepoint, 30 birds were randomly selected, necropsied and examined for grossly evident lesions. The small intestines of six broiler breeder hens from a single flock were collected post mortem and numbered A1-A6. Small intestines from the progeny of the same breeder flock were collected at hatch (H1-H6). The remaining cohort of chicks from the sampled breeder flock was placed in the lower story of a two-story commercial broiler house (hereafter called the monitored flock). A second flock was placed at the same time on the second story (hereafter called the control flock). The control flock was composed of chicks from five breeder sources. The monitored and control flocks received water from the same source and were given the same feed according to standard practice. Intestines were collected from six birds post mortem from the monitored and control flocks each when the birds were 2 weeks (monitored: 2W_1-2W_6 and control: 2W_C1-2W_C6), 4 weeks (monitored: 4W_1-4W_6 and control: 4W_C1-4W_C6) and 6 weeks (monitored: 6W_1-6W_6 and control: 6W_C1-6W_C6) of age.

Sample collection, processing and RNA extraction
The jejunum and ileum along with their contents were harvested from each carcass and processed as described previously [2]. RNA extraction from each viral pellet was performed as described previously [2] with one modification. The viral pellet was treated with 2 units of DNAaseI (New England Bio Labs, Ipswich, MA) per 100 μl of reaction and incubated at 37°C for 10 min, before RNA extraction. The quality of the total RNA was measured with a Nanodrop 2000 (ThermoScientific, Wilmington, DE).

Sequencing
Total RNA from each sample was converted to Illumina sequencing libraries using Illumina's TruSeq RNA Sample Prep Kit v2 according to manufacturer's instructions. Final library size distribution was validated and indexed, and the libraries from paired ends were pooled and size-selected at 320 bp +/-5% using Caliper's XT instrument giving an average insert size of approximately 200 bp. The pooled libraries were sequenced on a Hiseq 2000 with 100 bp PE run. De-multiplexed FASTQ files were used for subsequent analysis. The data have been deposited with links to BioProject accession number PRJNA310100 in the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/).

BLAST and taxonomy mapping
The quality of sequences from each of the samples was evaluated using FastQC. The paired reads from each sample were concatenated and one unique copy of each fragment was retained. The joined fragments were split to regenerate the forward (split R1) and reverse (split R2) reads. AbokiaBLAST, a commercial parallel implementation of NCBI BLAST was used for all blast searches against the NT database downloaded from NCBI in August 2013. The paired reads were blasted independently against the NT database using the program's default parameters. Due to computational resource limitations and the large dataset size, blast could not be completed on the samples A3, H3, H4, 4W_C2 and 6W_C2 and hence these samples were excluded from further analysis. There were numerous indications of sample saturation allowing the blast output to be subsampled at the 50% level to improve analysis speed. The assignment of taxonomic ID associated with each sequence and reconciliation of information from paired ends was performed using custom python scripts as described previously [2]. The brief schematics of the steps involved are shown in Fig 1. The reconciled output was used to distribute weight evenly to leaf nodes and the leaf nodes with weights below the cutoff were pruned [2]. MEGAN 4 software was used to visualize taxonomic weights and extract the information at family, genus and species levels.
Each sample had a unique sequencing depth and library size resulting in different numbers of reads with hits to RNA viruses for each sample. Table 1 lists the number of total and unique reads, and the reads with hits to viral clades and RNA viral families for half of the blast output. The H samples had the lowest reads with hits to RNA viruses. Thus, the rarefaction analysis for the H samples was done at a frequency of 1,000 sequences while for all other samples it was done at a frequency of 10,000 sequences.
To compare the alpha diversity across samples, the number of reads with hits to RNA viruses for each sample was standardized to the number of reads in the smallest library i.e. sample H1 ( Table 1). The alpha diversity measures, Chao I index, Sobs (number of observed genera) and the invsimpson index were computed (average of 1000 replicates) on 137,668 randomly sampled sequences from each sample. However, in all samples (except H), this level  of sub-sampling failed to saturate the number of genera present in the datasets likely leading to an underestimation of alpha diversity (Fig 2). To compare the alpha diversity of the remaining samples, the number of reads with hits to RNA viruses was sub-sampled to the number of reads in the 6W_C3 library i.e. 270,349 (Table 1). Sub-sampling at 137,668 sequences will be called small sub-sampling and sub-sampling at 270,349 sequences will be called large sub-sampling. The Chao I, Sobs and the invsimpson indices for the non-H samples were computed on the large sub-sample. Two-way ANOVA on the invsimpson indices was done using the commercially available GraphPad Prism 5 software to evaluate the respective impacts of bird age vs. breeder source (statistical model of 'bird age + breeder source + breeder source Ã age') on progeny viral diversity.
To measure beta diversity, the Bray-Curtis dissimilarity index [15] was computed on the small and large sub-samples from each sample. These dissimilarity matrices were used as input for ordination analyses and non-metric multidimensional scaling (NMDS) plots were created. To further determine if the spatial separation between the groups observed was statistically significant, analysis of molecular variance (AMOVA) was performed and p < 0.05 considered significant.

Inferential statistics
Differential abundance analysis at the family level was computed using edgeR package in R v.3.1.2 [16] with the relative log expression (RLE) method to normalize the data as described previously [17]. The statistical model of 'bird age + breeder source' was used. The RLE normalized data was also used for abundance profiles. The balanced design considered for the twoway ANOVA and differential abundance analysis with edgeR involved test birds from breeder A at 2 (2W), 4 (4W) and 6 (6W) weeks of age, and control birds from breeder A at 2 (2W_C), 4 (4W_C) and 6 (6W_C) weeks of age.

Flock health
The breeder, monitored and control flocks were within normal production parameters for the source company. There were no gross lesions observed in the 30 bird samples collected from the breeder hens, nor from the birds in the monitored and control flocks at 2, 4 or 6 weeks of age.

Diversity of RNA viruses in individual birds
For all the samples, increased sampling caused the rarefaction curves to plateau indicating that the majority of RNA viral diversity had been measured (Fig 2). The Chao I index, the richness estimator, was similar to Sobs, the number of observed genera (Fig 3A and 3B). The number of genera detected as well as the diversity increased with age in birds from both test and control flocks (Fig 3C). Further, a two-way analysis of variance performed on the invsimpson indices confirmed that age had a significant effect with a p-value of 0.004 whereas breeder source and the interaction (breeder source Ã age) had no significant effect indicating that the major driver of viral diversity was age (Fig 3D).

Comparison of composition and structure of RNA viral communities of all birds
The Bray-Curtis dissimilarity index matrices were used as input for ordination analyses. Nonmetric multidimensional scaling (NMDS) plots were created resulting in R-squared values of  (Fig 4). The ordination patterns with the small sub-sample showed that samples with respect to bird age were separated on axis 2 while the A and H samples were separated from all other samples on axis 1 (Fig 4A). The ordination patterns with the large sub-sample revealed that samples were separated on axis 1 with respect to age; and the A samples were separated from all other samples on axis 2 (Fig 4B). Irrespective of the number of sequences sub-sampled, a clear spatial separation based on bird age, but not breeder source was observed.
The results of AMOVA (Table 2) revealed that (a) the A group was significantly different from all other groups except the H group, (b) the H group was significantly different from all groups except the A group, (c) the 2W, 4W and 6W groups from the monitored flock were significantly different from each other, (d) for the control flock, the 4W_C group was significantly different from 6W_C group, and (e) between the monitored and control flocks, only the 4W group was significantly different from 4W_C group. These results confirmed that the A and H groups were different from all the other groups as observed in NMDS plots (Fig 4A). These results also indicated that age, and not the breeder source, had a significant effect in spatially separating the groups on NMDS plots (Fig 4A and 4B).

Differential abundance
A total of 333 genera belonging to 80 RNA viral families were detected across all the samples. Differential abundance analysis at the family level showed 21 viral families based on bird age ( Fig 5A) and two families based on breeder source (Fig 5B) were differentially abundant at FDR of 10% (q < 0.1). The RNA virus family Reoviridae was differentially abundant due to both age and the breeder source effects and was highly abundant at 2 weeks of age. The abundance of families such as Narnaviridae, Virgaviridae, Iflaviridae, Leviviridae, Togaviridae and unclassified RNA viruses, with very low abundance at 2 weeks of age, increased with age ( Fig  5A). These results reaffirmed our AMOVA and NMDS results.

Abundance profiles
The median of the abundance values was first calculated for each viral family followed by a calculation of the median of all values in the family-group matrix. Of 80 viral families detected, 45 had at least one value equal to or greater than median of all the values in the family-group matrix, and these were considered to be abundant (Fig 6). The abundance distribution of the families in the A and H groups were similar to each other but very different than the birds  from the test and control flocks at 2, 4 and 6 weeks of age. The abundance distribution of families in test and control birds at 2, 4 and 6 weeks of age were similar to each other. Many of the RNA viral families were present in low abundance in progeny birds at 2 weeks of age from either breeder source but their abundance increased with age. The top four highly abundant RNA viral families detected across all groups were Picornaviridae, Retroviridae, Flaviviridae and Astroviridae (Fig 6). The median abundance values for all birds were calculated for each virus family. Those which were present at greater than 1% in an age group are represented in Fig 7. The RNA viral families present in greatest abundance in the A and H groups were Picornaviridae and Retroviridae respectively; whereas the RNA viral family most abundant in 2, 4 and 6 week old birds was Astroviridae (Fig 7). Paramyxoviridae and Coronaviridae were also present in H group samples and their presence was detected most likely because of the vaccination for Newcastle disease and Infectious bronchitis viruses respectively in the hatchery [18]. The abundance of Astroviridae decreased with age while the abundance of Picornaviridae increased with age. The Retroviridae first appeared at 4 weeks of age and its abundance increased with age. Overall, total RNA viral diversity increased with age.

Discussion
It is easy to understand how bacteria can be linked to animal health and are part of a normal healthy state. RNA viruses, on the other hand, infect the host, presumably to cause disease. The Our findings demonstrate that not only are RNA viruses present in healthy animals but also that their community composition changes and diversity increases as birds age. The number of RNA viral families with > 1% abundance increased from three at 2 weeks of age to nine at 6 weeks of age, while the birds remained healthy. The Astroviridae, Picornaviridae and Retroviridae appeared to be the core members of the broiler gut RNA virome across all ages. The detection of Astroviridae family is in line with the findings of other investigators reporting the widespread prevalence of astroviruses in commercial chicken and turkey flocks from 2 to 6 weeks of age [12,13].
Enteritis is an age related disease in poultry, most often occurring between 2 and 4 weeks of age. In this study we found that at these ages, healthy broiler chicks already have an abundance Intestinal RNA Virus Community of Healthy Broiler Chickens of Astroviridae. So, although Astroviruses have been identified as etiologic agents of enteritis [8], our data suggest additional factors might be at play in development of clinical symptoms as other studies have suggested [12]. These findings may also suggest that the period of enteritis is coincidental with a period of low viral diversity in the gut which may be an indicator of host susceptibility. However, the temporal progression of RNA viral families in the poultry gut suggest that many of the paradigms of viral enteritis need to be reconsidered in light of this new information. Additional studies in which the progression of viral communities in flocks with viral enteritis and healthy flocks are compared are warranted.
The question of maternal influence on the bacterial community composition of progeny has been studied in mammals [19] but impacts on the wider microbial community, which includes viruses, has not been studied. Of course, mammals differ from birds in their relationship to their progeny and the link, for example, between bacterial community composition and type of birth [19], is not a factor in birds. In commercial settings like the one in this study, broiler eggs are removed from the breeder flock and taken to a hatchery where they are incubated and hatched and then to a third location where the progeny are raised. Given this process, it is likely that maternal influences on bacterial communities and also viral communities are different than in mammals. In this study, the RNA viral abundance profile, community distribution and composition of the samples from chicks on the day of hatch (H) was similar to that of their dams (A) but at 2, 4 and 6 weeks of age, it was more similar to their environmental neighbors i.e. the control flock. Based on these results we conclude that the breeder source influences the composition and structure of the RNA viral community at hatch but that the development of the RNA viral community after that is more strongly influenced by the bird's environment. In addition, while the hatch samples were most similar to the breeder source, they also contained viral RNA from families Paramyxoviridae and Coronaviridae, which were not detected at any other age. We attribute the presence of these RNA viruses to vaccination at the hatchery.
Studies on RNA viromes in animal hosts are new and will require many iterations before a full picture of the biology of the interactions between viral families and their hosts will emerge. This study provides a baseline looks at RNA viral communities and their progression in the broiler chicken host.