Virus-virus interactions and host ecology are associated with RNA virome structure in wild birds.

Little is known about the factors that shape the ecology of RNA viruses in nature. Wild birds are an important case in point, as other than influenza A virus, avian samples are rarely tested for viruses, especially in the absence of overt disease. Using bulk RNA-sequencing ("meta-transcriptomics"), we revealed the viral diversity present in Australian wild birds through the lens of the ecological factors that may determine virome structure and abundance. A meta-transcriptomic analysis of four Anseriformes (waterfowl) and Charadriiformes (shorebird) species sampled in temperate and arid Australia revealed the presence of 27 RNA virus genomes, 18 of which represent newly described species. The viruses identified included a previously described gammacoronavirus and influenza A viruses. Additionally, we identified novel virus species from the families Astroviridae, Caliciviridae, Reoviridae, Rhabdoviridae, Picobirnaviridae and Picornaviridae. We noted differences in virome structure that reflected underlying differences in location and influenza A infection status. Red-necked Avocets (Recurvirostra novaehollandiae) from Australia's arid interior possessed the greatest viral diversity and abundance, markedly higher than individuals sampled in temperate Australia. In Ruddy Turnstones (Arenaria interpres) and dabbling ducks (Anas spp.), viral abundance and diversity were higher and more similar in hosts that were positive for influenza A infection compared to those that were negative for this virus, despite samples being collected on the same day and from the same location. This study highlights the extent and diversity of RNA viruses in wild birds and lays the foundation for understanding the factors that determine virome structure in wild populations.

INTRODUCTION

Wild birds are ubiquitous, found on every continent, and a massive biomass of these animals moves across the globe on annual cycles of migration creating a truly interconnected planet (Bauer & Hoye, 2014). In addition to natural environments, birds can be found in our cities, using wastewater treatment plants, landfills and our drinking water reservoirs. Beyond wild birds, it is estimated that one out of every 7–14 birds on earth are raised for human consumption (i.e., chickens; Barrowclough, Cracraft, Klicka, & Zink, 2016; Food and Agriculture Organization of the United Nations, 2012), which may act as important amplifiers of potentially zoonotic avian viruses, such as influenza A virus (IAV; Gao et al., 2013; Wan, 2012; Yoon et al., 2015). Despite our important relationship with birds, we have only a limited understanding of the diversity of avian viruses. Indeed, most studies of avian viruses have focused on those that cause mass mortality in wild birds (e.g., Wellfleet Bay virus; Allison et al., 2014; Ballard et al., 2017), result in substantial economic losses in food production birds (e.g., avian avulavirus type 1; Alexander, 2000; Leighton & Heckert, 2007; Tolf et al., 2013) or are zoonotic (e.g., IAV; Gao et al., 2013; Wan, 2012; Yoon et al., 2015).

Avian viruses have a rich and complex ecology (van Dijk, Verhagen, Wille, & Waldenström, 2018), with patterns of prevalence affected by seasonality (Latorre‐Margalef et al., 2014), host species (Munster et al., 2007), host age (van Dijk et al., 2014), latitude (Lisovski, Hoye, & Klaassen, 2017) and urbanization (Wille, Lindqvist, Muradrasoli, Olsen, & Jarhult, 2017). However, although most studies of virus ecology and evolution have implicitly assumed a “one‐host, one‐virus” model of host–pathogen interactions, such as the Mallard (Anas platyrhynchos)–IAV model (van Dijk et al., 2014), both hosts and their viruses exist in communities, and it is likely that these communities are the drivers of viral ecology. For example, despite intensive focus on the Mallard–IAV system (Latorre‐Margalef et al., 2014; van Dijk et al., 2014), IAV is in reality a multihost virus detected in over 100 species of wild birds, with different avian species playing different roles in virus ecology and evolution (Olsen et al., 2006). For example, gulls are reservoirs for evolutionary distinct IAV subtypes (Wille et al., 2011), and rare subtypes may be maintained in members of the Anseriformes and Charadriiformes that are infrequently sampled (Wille et al., 2018). In turn, numerous viruses have been detected in wild bird populations, and Mallards may be co‐infected with at least three different RNA viruses simultaneously in the absence of overt signs of disease (Wille et al., 2015, 2017). These other avian RNA viruses—avian coronavirus and avian avulavirus type 1—have seasonal prevalence patterns that generally mirror that of IAV, and it is therefore possible that they may also share similar host taxonomic or geographic differences in viral community structure (Wille et al., 2017). Virus co‐infection may also be an important driver of viral prevalence, as virus–virus interactions may enhance or interfere with infection (Diaz‐Munoz, 2017; Elena & Sanjuan, 2005; Henle, 1950; Jolly & Narayan, 1989). As a case in point, a higher prevalence of avian coronavirus was found in a bird population experiencing IAV infection, suggesting that the latter might play a role in structuring avian viromes in general (Wille et al., 2015).

Although one in 10 bird species is found in Australia, we know little of the accompanying viral diversity and abundance in these animals, nor of the large‐scale ecological factors that determine virome composition. We used a recently developed unbiased meta‐transcriptomic pipeline based on bulk RNA‐sequencing (Shi, Zhang, & Holmes, 2018) to reveal the viromes of four Australian avian species and to evaluate how the structure of entire viral communities is impacted by a variety of ecological correlates. In particular, we assessed the role of host taxonomy, location and co‐infection with IAV on virome structure. This study illustrates the extent of RNA viral diversity in wild birds, and the importance of the expansion of traditional host–pathogen systems beyond simple one‐host, one‐virus systems to disentangle ecological processes in viral presence and abundance.

MATERIALS AND METHODS

Ethics statement

This research was conducted under approval of the Deakin University Animal Ethics Committee (permit numbers A113‐2010 and B37‐2013). Banding was performed under Australian Bird Banding Scheme permit (banding authority numbers 2915 and 2703). Research permits were approved by the Department of Environment, Land, Water and Planning, Victoria (permit numbers 10006663 and 10005726), Department of Environment, Water and Natural Resources, South Australia (research permit numbers M25919‐1,2,3,4,5), and the Department of Primary Industries, Parks, Water and Environment, Tasmania (permit number FA 13032).

Sample selection

Samples were collected as part of a long‐term IAV surveillance study (Ferenczi, 2016; Ferenczi et al., 2016). Birds were captured using baited funnel walk‐in traps, cannon nets or mist nets as described previously (Ferenczi, 2016). Importantly, none of the birds in this study exhibited any signs of disease. Samples were collected from (a) three temperate locations in Australia—the Western Treatment Plant near Melbourne (37°59′11.62′′S, 144°39′38.66′′E), Yallock Creek (38°13′51.6″S 145°28′43.9″E) and King Island (39°55′52″S 143°51′02″E), and (b) an interior arid location—Innamincka Regional Reserve (27°32′28″S 140°35′47″E). Species selected for the study included both those known to be important in IAV ecology (Anas ducks and Ruddy Turnstone) and those in which IAV has not been described (Australian Shelduck and Red‐necked Avocet; Table 1).

Table 1

Eight libraries were sequenced reflecting different avian taxonomy, location within Australia and influenza A infection status

Taxonomy	Host species	Location	Sampling location	Influenza A status	Sample type	Year
Anseriformes, Anatidae, Tadorninae	Australian Shelduck (Tadorna tadornoides)	Temperate	Western Treatment Plant, Victoria	Negative	Cloacal	2012
Anseriformes, Anatidae, Anatinae	Dabbling duck (Anas sp.)	Temperate	Western Treatment Plant, Victoria	Positive	Cloacal	2013
	Dabbling duck	Temperate	Western Treatment Plant, Victoria	Negative	Cloacal	2013
	Dabbling duck	Interior	Innamincka Regional Reserve, South Australia	Negative	Cloacal	2013
Charadriiformes, Recurvirostridae	Red‐necked Avocet (Recurvirostra novaehollandiae)	Temperate	Yallock Creek, Victoria	Negative	Cloacal	2013
	Red‐necked Avocet	Interior	Innamincka Regional Reserve, South Australia	Negative	Cloacal	2013
Charadriiformes, Scolopacidae	Ruddy Turnstone (Arenaria interpres)	Temperate	King Island, Tasmania	Negative	Combined oropharyngeal/cloacal	2014
	Ruddy Turnstone	Temperate	King Island, Tasmania	Positive	Combined oropharyngeal/cloacal	2014

Cloacal samples (from 2012 to 2013) or a combination of oropharyngeal and cloacal samples (from 2014) were collected using a sterile‐tipped swab and was placed in viral transport media (VTM, Brain‐heart infusion broth containing 2 × 106 U/L penicillin, 0.2 mg/ml streptomycin, 0.5 mg/ml gentamicin, 500 U/ml amphotericin B, Sigma). All samples were assayed for IAV as described previously (Ferenczi et al., 2016).

RNA library construction and sequencing

RNA was extracted using the MagMax mirVana™ Total RNA Isolation Kit (Thermo Fisher Scientific) on the KingFisher™ Flex Purification System (Thermo Fisher Scientific). All extracted samples were assessed for RNA quality using the TapeStation 2200 and High Sensitivity RNA Reagents (Agilent Genomics, Integrated Sciences), and 10 samples with the highest concentration were pooled (based on species, location and IAV infection status) using equal concentrations and subsequently concentrated using the RNeasy MinElute Cleanup Kit (Qiagen). Libraries were constructed using the TruSeq total RNA library preparation protocol (Illumina), and host rRNA was removed using the Ribo‐Zero Gold Kit (Illumina). Paired‐end sequencing (100 bp) of the RNA library was performed on the HiSeq2500 platform. All library preparation and sequencing were carried out at the Australian Genome Research Facility (AGRF, Melbourne).

RNA virus discovery

Sequence reads were demultiplexed and trimmed with Trimmomatic followed by de novo assembly using Trinity (Grabherr et al., 2011). No filtering of host/bacterial reads was performed before assembly. All assembled contigs were compared to the entire nonredundant nucleotide (nt) and protein (nr) database using blastn and diamond blast (Buchfink, Xie, & Huson, 2015), respectively, setting an e‐value threshold of 1 × 10−10 to remove potential false positives. Abundance estimates for all contigs were determined using the RSEM algorithm implemented in Trinity. All contigs that returned blast hits with paired abundance estimates were filtered to remove plants and invertebrate reads that likely correspond to the host diet, as well as fungal, bacterial and host sequences. Blast results were used to initially classify viruses to their appropriate family and genus level, and the virus list was further filtered to remove viruses with invertebrate (Shi et al., 2016), plant or bacterial host associations using the Virus‐Host database (http://www.genome.jp/virushostdb/).

To compare relative viral abundance across libraries, three host reference genes were mined from the data using a custom blast database. As not all host reference genes are stably expressed, we utilized three genes that are stably expressed in the Mallard (Anas playrhynchos) lower gastrointestinal tract (Chapman et al., 2016): ribosomal protein L4 (RPL4), ribosomal protein S13 (RPS13) and NADH dehydrogenase 1 alpha subcomplex (NDUFA) from both Mallard (taxid: 8839), chicken (Gallus gallus; taxid: 9031) and Zebra Finch (Taeniopygia guttata; taxid: 59729).

Virus genome annotation and phylogenetic analysis

Contigs >1,000 bp in length were inspected using Geneious R10 (Biomatters, New Zealand), and open reading frames corresponding to predicted genome architectures based on the closest reference genomes were interrogated using the conserved domain database (CDD, https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) with an e‐value threshold of 1 × 10−5. Reads were subsequently mapped back to viral contigs to identify mis‐assembly using the Geneious mapping function. Viruses with full‐length genomes, or incomplete genomes but that possess the full‐length RNA‐dependent RNA polymerase (RdRp) gene, were used for phylogenetic analysis. Briefly, amino acid sequences of the polyprotein or gene encoding for the RdRp were aligned using mafft (Katoh & Standley, 2013), and gaps and ambiguously aligned regions were stripped using trimAL (Capella‐Gutierrez, Silla‐Martinez, & Gabaldon, 2009). Final alignment lengths are presented in Supporting Information Table S2. The most appropriate amino acid substitution model was then determined for each data set, and maximum‐likelihood trees were estimated using PhyML 3.0 (Guindon et al., 2010) with 1000 bootstrap replicates using the ATGC server (http://www.atgc-montpellier.fr/phyml/execution.php). For IAV and gammacoronavirus, phylogenies were also estimated using the nucleotide sequences of complete or partial reference genome sequences to better place viruses in context of currently described avian viral diversity. Similarly, the best‐fit model of nucleotide substitution was selected, and maximum‐likelihood trees were estimated using PhyML 3.0 with 1,000 bootstrap replications. Novel viral species were identified as those that had <90% RdRp protein identity, or <80% genome identity to previously described viruses.

Diversity and abundance across libraries

Relative virus abundance was estimated as the proportion of the total viral reads in each library (excluding rRNA). All ecological measures were calculated using the data set comprising “higher” vertebrate‐associated viruses (i.e., those associated with birds and mammals), albeit with all retroviruses and retrovirus‐like elements removed (hereafter, avian virus data set). Analyses were performed using r version 3.4.0 integrated into RStudio version 1.0.143, and plotted using ggplot2. Specifically, both the observed richness and Shannon effective [alpha diversity] were calculated for each library at the family and genus levels using the Rhea script sets (Lagkouvardos, Fischer, Kumar, & Clavel, 2017). Beta diversity was calculated using the Bray–Curtis dissimilarity matrix using the vegan package (Oksanen et al., 2007) at both family and genus levels, and presented as a network using the igraph (Csardi & Nepisz, 2006) and ggnet (https://github.com/briatte/ggnet) packages. Nonmetric multidimensional scaling (NMDS) ordination was additionally calculated based on Bray–Curtis dissimilarity and Adonis tests (PERMANOVA) applied using the phyloseq package (McMurdie & Holmes, 2013).

RESULTS

RNA‐Seq as a means to identify avian viruses

We characterized the total transcriptome of eight avian pools, representing two tribes in the order Anseriformes (waterfowl) and two families in the order Charadriiformes (shorebirds). These pools were designed to answer specific questions on the determinants of virome structure including the impact of bird taxonomy (within and between Anseriformes and Charadriiformes), location (temperate vs. arid sampling sites) and effect of IAV infection (Table 1). Each library comprised swab samples collected from 10 individuals at the same time point and location, therefore increasing the chances of finding viruses at lower prevalence.

RNA‐sequencing of rRNA‐depleted libraries resulted in a median of 44,345,130 (range 39,267,372–47,650,666) reads per pool, which were assembled into a median of 175,559 contigs (range 135,254–357,869). An assessment of the host reference gene RPS13, a proxy for sequencing depth of libraries, revealed similar abundances (0.000102%–0.000342% of reads), suggesting similar host sequencing depth across the libraries (Supporting Information Figure S1). These eight libraries had marked differences in the abundance of avian viral reads; Ruddy Turnstones (Arenaria interpres) and Anas ducks that were IAV positive (0.21% and 0.1% of reads), and Red‐necked Avocets (Recurvirostra novaehollandiae) from the interior (0.26% of reads) had relatively high abundances of avian viruses, while lower abundance levels were observed in Ruddy Turnstones and Anas ducks that were IAV negative (0.00006% and 0.00051% of reads, respectively; Figure 1).

Figure 1

Abundance of viruses in each library. (a) Abundance of all exogenous viruses including those from avian, invertebrate, lower vertebrate, plant, fungi or bacterial hosts. (b) Abundance of all viruses that are associated with birds. (c) Abundance of retroviruses or retrovirus‐like elements that have avian or mammalian signature. (d). Host reference gene RPS13

Blast analysis and virus characterization revealed the genomes of 27 RNA viruses, of which 18 were newly described species—that is, of sufficient phylogenetic distinction to represent new virus species—but that belonged to existing families and were most closely related to other avian viruses (Supporting Information Table S1). The viruses identified comprised double‐stranded RNA viruses (Reoviridae, Rotavirus; Picobirnaviridae), positive‐sense single‐stranded RNA viruses (Caliciviridae; Astroviridae; Picornaviridae, genus Avihepatovirus, Megrivirus and Unassigned genera; Coronaviridae, genus Gammacoronavirus) and negative‐sense single‐stranded RNA viruses (Rhabdoviridae, genus Tupavirus; Orthomyxoviridae, genus Influenza A virus). Members of the family Paramyxoviridae, known to circulate in wild birds (Ramey et al., 2013; Wille et al., 2015), were notably absent. No DNA viruses (i.e., the RNA transcripts of DNA viruses) were detected, potentially because cloacal and oropharyngeal samples are a richer source of shed viruses rather than those actively replicating within cells of the gastrointestinal tract. An array of retroviruses or retrovirus‐like elements was also detected, but due to the challenge in differentiating between endogenous and exogenous retroviruses they will not be discussed here (Figure 1).

Substantial undescribed diversity of RNA viruses in wild birds

An average of 80% of virus species in each library was novel (range 50%–100%), and in three libraries, all viruses were novel (Supporting Information Table S1), illustrating the large undiscovered viral diversity in wild birds. Numerous new viral species were identified from viral families that are not frequently screened for in wild birds, including rhabdoviruses, caliciviruses, picornaviruses and rotaviruses. In the case of some viral families, complete viruses were only found in a single species, such as the picobirnaviruses detected in Australian Shelducks (Tadorna tadornoides; Figure 2a). Other viruses, such as the caliciviruses, were highly abundant across all avian hosts, and full genomes were found in all avian species included in this study (Figure 4, Supporting Information Figure S9).

Figure 2

Phylogenies of double‐stranded RNA viruses. These trees show (a) segment 2 (RdRp) of picobirnaviruses and (b) the VP1 segment (RdRp) of the rotaviruses described in this study. All phylogenetic trees were midpoint rooted for clarity only. The scale bar indicates the number of amino acid substitutions per site. Bootstrap values >70% are shown for key nodes. Viruses described in this study are marked with a filled circle [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Across all the RNA‐dependent RNA polymerase (RdRp) phylogenies, the viral species from wild birds generated in this study were in the most part similar to previously described avian viruses, often forming an apparent “avian” clade within each group of viruses and suggestive of a relatively long‐term association with birds (Figures 2, 3, 4, 5, 6). Although there were some exceptions, such as the Shelduck picobirnaviruses which were most closely related to those viruses sampled in swine viruses, it is possible that this simply reflects poor sampling.

Figure 3

Partial RpRp phylogeny of members of the genus Avastrovirus. The tree is rooted between the avian and mammalian astroviruses. The scale bar indicates the number of nucleotide substitutions per site. Bootstrap values >70% are shown for key nodes. Viruses described in this study are marked with a filled circle. The phylogeny of the full‐length polyprotein is presented in Supporting Information Figure S6 [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Figure 4

Phylogenetic tree of the polyprotein, containing the RdRp, of members of the Caliciviridae. The most divergent calicivirus, Atlantic Salmon calicivirus, was used as an outgroup to root the tree. The scale bar indicates the number of amino acid substitutions per site. Bootstrap values >70% are shown for key nodes. Viruses described in this study are marked with a filled circle [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Figure 5

Phylogeny of the L gene (RdRp) of members of the Rhabdoviridae. Almendraviruses were set as the outgroup, and representative viruses for each genus (as per Walker et al., 2015) were also included in the analysis. The scale bar indicates the number of amino acid substitutions per site. Bootstrap values >70% are shown for key nodes. Viruses described in this study are marked with a filled circle [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Figure 6

Phylogeny of the polyprotein, containing the RdRp, of the Picornaviridae. Reference viruses are those from Boros et al. (2016), and the tree was midpoint rooted for clarity only. The scale bar indicates the number of amino acid substitutions per site. Bootstrap values >70% are shown for key nodes. Viruses described in this study are marked with a filled circle. The tree is midpoint rooted for clarity only [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Wild birds have previously been described as hosts for coronaviruses, astroviruses and IAVs, and the viruses identified in this study belonged to “wild bird” clades to which sequences from poultry fell as outgroups (Supporting Information Figures S3–S8). Specifically, four different IAV HA‐NA subtypes were found in this study; H12N5 and H9N3 viruses from Anas ducks, and H6N8 and H10N8 from Ruddy Turnstones (Supporting Information Figures S3–S5). In addition to H12N5 and H9N3 in Anas ducks, one short contig from the HA of H10 influenza virus was also identified. These three subtypes, all identified in Anas ducks, had different abundances in this library: H12N5 comprised 77% of all avian viral reads (0.081% of all reads in the library), compared to H9N3 that comprised 0.5% of avian viral reads in the library (0.0005% of all reads), and H10 that represented only 0.0068% of avian reads. This is in comparison with the H6N8 and H10N8 viruses identified in the Ruddy Turnstones that had similar abundances (25%–35% of avian viral reads, 0.05%–0.07% of all reads in the library).

Broadly, the IAVs from Australian birds described here were most similar to viruses sampled from Eurasian wild birds, which is expected given that Australian migratory birds use the East Asian–Australian flyway. However, while the H12 virus was more similar to viruses from Eurasia, this virus was phylogenetically distinct from currently circulating viruses, suggesting the presence of a potential “Australia‐specific” clade. Additionally, the Ruddy Turnstone H10 sequence fell into the North American clade rather than the Eurasian clade, in contrast to the N8 segment which fell into a Eurasian clade. Such a phylogenetic pattern is indicative of the intercontinental reassortment of these Ruddy Turnstone viruses. The gammacoronavirus in this study, identified in Red‐necked Avocets, was related to circulating wild bird gammacoronaviruses from waterfowl from both Eurasia and the United States (Supporting Information Figures S7 and S8).

Novel dsRNA viruses

Two complete picobirnavirus genomes were found in wild bird samples (Figure 2a). These viruses, all from Australian Shelducks, clustered together on the phylogenetic tree, although a more divergent partial virus was also found. Australian Shelducks are particularly rich in picobirnavirus diversity, with 21% of all avian viral reads in the library derived from picobirnaviruses (0.002% of total reads), although this virus family was also found at low abundance in Anas ducks from the interior. This potentially suggests a preference for the Anseriformes, although this will need to be confirmed with more data (Figure 2a, Supporting Information Figure S9). In addition, rotaviruses were found in almost all libraries and in all host groups (Supporting Information Figure S9) with three different subtypes revealed (Figure 2b). Specifically, we found previously described rotaviruses D and G in apparently healthy ducks and avocets, respectively, even though they are known to cause enteritis in poultry. Unlike other phylogenies in which wild bird viruses fell in clades that are distinct from poultry‐associated viruses, wild bird and poultry rotaviruses were phylogenetically similar and hence clustered on the tree, although sample size was limited. Ruddy Turnstones, however, carried a highly abundant (53% of avian viral reads, 0.11% of total reads in the library) and highly divergent rotavirus, characterized by a long branch, that fell as a sister species to rotavirus G (Figure 2b).

Novel ssRNA viruses

Two novel avastroviruses were identified in Red‐necked Avocets, both falling as outgroups to Group 2 viruses, the archetype of which is avian nephritis virus (Figure 3, Supporting Information Figure S6). These viruses share 60% and 40% pairwise amino acid identity to avian nephritis virus, respectively, suggesting that there is a large undescribed diversity of wild bird avastroviruses. Calicivirus reads were identified in all libraries, with the exception of Anas ducks from the interior of Australia. Furthermore, full genomes of five caliciviruses were identified in four libraries at high abundance (Australian Shelducks 0.005% of total reads, Ruddy Turnstone IAV positive 0.006% of total reads, Anas duck IAV positive 0.002% of total reads and Avocet interior 0.03% of total reads), and all these viruses belonged to the same clade as currently described poultry viruses within an unassigned genus (Figure 4). Two novel rhabdoviruses from Anseriformes found in this study fell as a divergent group within the genus Tupaviruses, within which Durham virus is the only avian virus previously described (Figure 5). Specifically, Shelduck rhabdovirus and duck rhabdovirus fell as relatively distantly related sister species, as illustrated by long branch lengths on the phylogeny, and thereby potentially represent a novel clade of wild bird viruses.

Our virome sampling also revealed a great diversity of picornaviruses (Picornaviridae), almost all of which were found in Red‐necked Avocets (Figure 6). Megriviruses were remarkably abundant in avocets, comprising 30% and 65% of all avian viral reads (0.004% and 0.169% of total reads) from birds sampled in two locations, representing two locations. Furthermore, megriviruses represent the only virus (other than IAV) found in more than one library (Figures 6 and 7). The library generated from avocets in the interior also contained (a) Avocet picornavirus B‐like A that formed a sister group to pigeon picornavirus B; (b) Avocet picornavirus B‐like B, a sister group to a clade containing both pigeon picornavirus B and Avocet picornavirus B‐like A; and (c) Avocet picornavirus, a highly divergent virus that likely represents a novel genus. Finally, in Anas ducks, we observed a divergent sister group to duck hepatitis A 1 and 3 (Wild Duck avihepatovirus‐like; Figure 6).

Figure 7

Composition of avian viromes. (a) Bipartite network illustrating species for which complete viral genomes found in each library. Each library is represented as a central node, with a pictogram of the avian species, surrounded by each viral species. Where no complete viral genomes were revealed, the pictogram is shown with no viruses. Where two libraries share a virus species, the networks between the two libraries are linked, and the edges are thicker for aesthetic purposes. Placement of libraries is arranged by influenza A status on the y‐axis and location on the x‐axis. Virus colour corresponds to virus taxonomy. A list of viruses from each library is presented in Supporting Information Table S1, and phylogenetic trees for each virus family and species can be found in (Figures 2, 3, 4, 5, 6, Supporting Information Figures S3–S8). (b) Nonmetric multidimensional scaling (NMDS) plot (applying the Bray–Curtis dissimilarity matrix) for viral abundance and virus family diversity. Colour, shape and fill correspond to host species, influenza A infection status and location, respectively. For increased clarity, influenza A‐positive libraries are indicated. Additional NMDS plots where data are analysed with and without IAV reads at both the viral genus and family level are found in Supporting Information Figures S12 and S13 [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Factors affecting the structure and abundance of avian viromes

One of the most important results of our study was that IAV status and location, but not host taxonomy, were associated with differences in viral abundance and diversity. Because of the potential impact of host phylogeny, we expected that virome structure would be similar within the Anseriformes and Charadriiformes but differ between these orders. However, across all the libraries and controlling for location and IAV status, libraries from members of the Anseriformes were no more similar to each other than they were to those from the Charadriiformes. This lack of taxonomic distinction was apparent whether the analysis was performed at the level of viral species, genus or family (Figure 7, Supporting Information Figures S10–S12, viral family, R 2 = 0.142, p = 0.353; viral genus R 2 = 0.153, p = 0.251), although the sample size was relatively small. To better control for other variables in our sampling scheme, we compared the libraries from avocets and dabbling ducks as these were all IAV negative and sampled from the same locations; in this case, the relationship remains statistically insignificant (viral family, R 2 = 0.25, p = 0.4; viral genus R 2 = 0.28, p = 0.3), although we lose statistical power due to the small sample size. Furthermore, the three libraries from Anas ducks had a different abundance and viral composition, and these three libraries represent differences in IAV infection status and location (Figure 7, Supporting Information Figures S10–S12). No viral species were shared within host species or family, with the exception of IAV and a megrivirus found in both avocet libraries (Figure 7). Finally, there were no differences in viral family distribution at the level of host species and order; all viral families were found in both Anseriformes and Charadriiformes, with the exception of picobirnaviruses which only occurred in the Anseriformes (Supporting Information Figure S12).

It might also be expected that birds in temperate locations would have a higher viral abundance and diversity than birds from the arid interior of Australia given greater prevalence of IAV in temperate latitudes (Lisovski et al., 2017). However, the library from Red‐necked Avocets from the arid interior had a higher viral abundance and diversity than individuals from temperate Australia and also had the highest viral abundance across all libraries (0.26% total reads; Figures 1, 7 and 8, Supporting Information Figures S10 and S11). There was also a clear virome difference between Anas ducks across locations: a higher viral abundance and species diversity were found in the ducks from the interior compared to the temperate ducks that were negative for IAV. However, temperate ducks that were IAV positive had a higher viral diversity and abundance compared to ducks sampled from the interior (Figure 8a,b). Overall, incorporating all libraries and controlling for IAV infection status and host species, location did not predict higher similarity between the libraries, as libraries from the same location were no more similar to each other than those from different locations (viral family, R 2 = 0.1, p = 0.554, viral genus R 2 = 0.093, p = 0.8; Figures 7 and 8, Supporting Information Figures S10–S12). When comparing only dabbling ducks that were IAV negative and avocets from arid and temperate locations, the relationship remained statistically insignificant (viral family R 2 = 0.199, p = 0.9, viral genus R 2 = 0.2, p = 1); however, there was limited statistical power due to small sample size.

Figure 8

Influenza A status and location are associated with differences in viral abundance and diversity. (a, c, e) correspond to the influenza A virus infection status, while (b, d, f) correspond to location. (a, b) Avian viral abundance in libraries in grey, and in (a), abundance of IAV is indicated in black. (c, d) Abundance of host reference gene RSP13. (e, f) Heatmap illustrating viral diversity, at the genus level in each library, with colour corresponding to viral abundance. Blue and purple correspond to viruses with high abundance, and pink corresponds to viruses with low abundance. Asterisks indicate cases in which at least one complete viral genome was obtained [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Finally, we expected that, in accord with previous studies (Wille et al., 2015, 2017), libraries containing IAV would have higher viral diversity compared to those that were negative for IAV. To address this, samples from IAV positive and negative birds were selected from the same location during the same sampling expedition to remove any potential bias. Libraries from both Ruddy Turnstone and Anas ducks that were positive for IAV indeed had a higher viral abundance (Figure 8a) [0.21%, 0.1% compared to 0.000061%, 0.0005% viral reads] and virus diversity (Figures 7 and 8c, Supporting Information Figures S10 and S11), at the family, genus and species levels (Figures 7 and 8c, Supporting Information Figures S10 and S11). This trend remained when abundance or diversity attributable to IAV was removed from the analysis. Furthermore, the two libraries containing IAV were statistically significantly more similar to each other in abundance and composition compared to all other libraries sequenced (viral family including IAV reads, R 2 = 0.24, p = 0.008; viral genus including IAV reads, R 2 = 0.24, p = 0.014; viral family excluding IAV reads, R 2 = 0.17, p = 0.04; viral genus excluding IAV reads R 2 = 0.24, p = 0.017; Figures 7 and 8, Supporting Information Figures S11 and S12). Only including the turnstone and Anas duck libraries negatively affects statistical power due to small sample size.

DISCUSSION

Despite the ubiquitous nature and economic importance of birds, we have a poor understanding of the natural viral diversity in this major animal phylum. To this end, we employed an unbiased metagenomics approach to reveal avian viromes, comprising 27 novel and previously described viral species, in a framework of ecological hypothesis testing.

Given the long‐term association between hosts and viruses, it was not unexpected that the viruses revealed in this study were most closely related to other avian viruses (Shi, Lin, et al., 2018), especially virulent poultry viruses that have been an important focus in virus characterization efforts (Boros et al., 2016; Day, Ballard, Duke, Scheffler, & Zsak, 2010). Based on previous studies, we anticipated finding low pathogenic forms of coronaviruses, astroviruses and avulaviruses (Wille et al., 2015, 2017). While we did detect most of these viruses, the absence of avian avulavirus was surprising. Avian avulavirus type 1 is present in wild birds in Australia (Hoque, Burgess, Greenhil, Hedlefs, & Skerratt, 2012; Hoque, Burgess, Karo‐Karo, Cheam, & Skerratt, 2012; Hore, Campbell, & Turner, 1973; Mackenzie, Edwards, Holmes, & Hinshaw, 1984; Peroulis & O'Riley, 2004) and globally (Alexander, 2000; Austin & Webster, 1993; Hanson et al., 2005; Ramey et al., 2013; Tolf et al., 2013; Wille et al., 2015), and has been detected in co‐infection studies (Wille et al., 2015, 2017). We also genomically described 18 new viral species which belonged to previously identified avian clades (genera) predominantly comprised of poultry viruses. For example, we identified five caliciviruses from wild birds that belonged to a previously described avian clade in the Caliciviridae comprised of chicken, turkey and goose caliciviruses (Liao, Wang, Wang, & Zhang, 2014; Wolf, Reetz, & Otto, 2011; Wolf et al., 2012).

While virus species and genotypes that were sister groups to poultry viruses were revealed, it is important to note all samples in this study were collected from birds with no clinical signs of disease. This raises two important issues. First, with the sequencing of more wild birds, those clades formally dominated by poultry will likely expand to include many viral species and genotypes from wild birds. This will be central to a better understanding the movement of avian viruses between wild bird reservoirs and poultry populations and hence of disease emergence in general. Indeed, poultry production has rapidly expanded in the last century (Kaleta & Rulke, 2008), to the extent that ~70% of avian biomass on the planet are now poultry (Bar‐On, Phillips, & Milo, 2018), creating a relatively new, but large niche for viruses. Furthermore, unlike wild birds, in viruses adapted to poultry such as Marek's disease virus (a double‐strand DNA virus), there has likely been selection for high transmissibility and high virulence (Rozins & Day, 2017). As such, sequencing wild bird viruses is imperative in understanding the evolutionary processes involved.

Second, these data raise the issue of how wild birds are able to tolerate such high levels of virus diversity and abundance seemingly in the absence of overt disease (Medzhitov, Schneider, & Soares, 2012; Råberg, 2014). In particular, cloacal swabs of Red‐necked Avocets from the interior had a high viral abundance (0.26% of reads were from avian viruses) and these 10 birds shed 13 viral genera and eight viral species for which full genomes were revealed. Ruddy Turnstone and Anas ducks that were infected with IAV similarly shed avian viruses at high abundance (0.21% and 0.1%, respectively), albeit with lower viral diversity. Although this must impose some physiological effect on the host, there continues to be conflicting data on the physiological effect of IAV infection in isolation (Kuiken, 2013), let alone the viral abundance described in this study. A large viral diversity in healthy, individual wild birds (Fawaz et al., 2016; Wille et al., 2015) and poultry (Day, Oakley, Seal, & Zsak, 2015; Day et al., 2010; Lima et al., 2017) has been previously described, in contrast to chickens in which viral diversity was described in diseased animals (e.g., diarrhoea; Boros et al., 2016). The leading hypothesis reflects a long history of host–pathogen co‐evolution: chickens are a relatively new host niche for IAV, and following spillover from wild birds highly pathogenic phenotypes evolve (e.g., H7; Seekings et al., 2018), sometimes resulting in catastrophic mortality. This is in contrast to wild birds that have likely been co‐evolving with IAV over long time periods, with natural selection perhaps favouring lower virulence; as a consequence, the highly pathogenic IAV found in wild birds is most likely due to spillover from poultry (Barber, Aldridge, Webster, & Magor, 2008; Little, Shuker, Colegrave, Day, & Graham, 2010; van Dijk, Fouchier, Klaassen, & Matson, 2015). Indeed, a muted inflammatory response translating to immunological tolerance to viral infections may allow some hosts, such as bats, to harbour a variety of viruses (Brook & Dobson, 2015; Xie et al., 2018). Similarly, Pekin ducks (Anas platyrhynchos domesticus) and wild strain Mallard ducks appear to have a controlled innate immune response against both low pathogenic IAV (Helin et al., 2018) and highly pathogenic IAV (Saito et al., 2018), with upregulation of the innate immune system occurring on day 1 postinfection and no evidence of hypercytokinaemia, or “cytokine storms.” Given the high (>0.1% of reads) viral abundance in some libraries, it is possible that some of the viruses described here do not cause disease in the absence of some other physiological or environmental stressor, although this is clearly an issue that needs to be explored in more detail.

Viral co‐infection is likely to be the rule rather than the exception and is shaped by both host ecology and virus–virus interactions (Diaz‐Munoz, 2017). To date, much of what we know about viruses in wild birds is derived from many years of research in IAV (Ferenczi et al., 2016; Latorre‐Margalef et al., 2014; Munster et al., 2007; Olsen et al., 2006) and those RNA viruses that have similar patterns of host preference and seasonality (Wille et al., 2015, 2017). Members of the Anseriformes and Charadriiformes have proven to be excellent model species for this study, and we detected viruses previously described in these hosts (Chu et al., 2011; Muradrasoli et al., 2010; Wille, Muradrasoli, Nilsson, & Jarhult, 2016; Wille et al., 2017) as well as a suite of novel viruses. One of the key observations of our study is that avian taxonomy did not drive virome structure; that is, there was no specific clustering of libraries in the NMDS plots by avian order (Anseriformes vs. Charadriiformes), and the three Anas duck libraries were different, although this analysis had limited statistical power. Given that IAV is prevalent in both these avian orders, we suggest that host ecology may play a more important role than host taxonomy in shaping virome diversity. For example, it is possible that waterbirds share viral families, genera and species as shallow water bodies facilitating virus transmission between individuals, as with IAV (Hoye, Fouchier, & Klaassen, 2012; van Dijk et al., 2018). In support of this, there was a difference between birds sampled (Red‐necked Avocets in particular) in lakes of the arid interior as compared to temperate coastlines.

The ecological factors assessed here are not mutually exclusive, as shown by the complex relationship between the three Anas duck libraries which had different structures given different conditions (location and IAV status). Samples from Ruddy Turnstones were collected from the same beaches on the same sampling trip, yet the 10 birds positive and negative for IAV had different viromes (total abundance, species abundances and viral diversity). Furthermore, the viromes of Anas ducks and Ruddy Turnstones that were positive for IAV were more similar to each other than to all other libraries. Anecdotally, while we were able to successfully culture H6N8 and H10N8 virus as part of an ongoing IAV surveillance project, the viruses from the Anas ducks were not successfully isolated, demonstrating the power of the meta‐transcriptomic approach used here. These trends, however, may be biased due to sample pooling, such that the patterns may be due to only a few individuals in the pool. Ultimately, therefore, the validity of the patterns observed here needs to be re‐assessed on the basis of individual transcriptomes, although such work will obviously be both costly and time‐consuming. In addition, it was previously shown (Wille et al., 2015) that co‐infection with IAV was important in modulating the prevalence of other RNA viruses. In particular, IAV modulates the interferon response of the host, changing the antiviral state (Garcia‐Sastre, 2001, 2011; Hale, Randall, Ortin, & Jackson, 2008), and this may promote co‐infection or prevent viral clearance of certain viruses. We have a surprisingly poor understanding of virus–virus interactions, although viruses do have mechanisms to mediate infection by other viruses. For example, viruses may have synergistic (enhancing) or antagonistic (inhibiting) interactions, and this may occur within and/or across viral species (Diaz‐Munoz, 2017; Elena & Sanjuan, 2005; Henle, 1950; Jolly & Narayan, 1989). Regardless, virus–virus interactions are important drivers of co‐infection and may be further affected by virus and host ecology (Diaz‐Munoz, 2017).

In sum, we have expanded our understanding of the diversity of avian viruses and laid the foundation for future hypothesis testing of the factors associated with virome structure in wild birds using high‐throughput meta‐transcriptomics. This study focused on avian orders known to be central in the ecology of IAV, but also a number of other RNA viruses. Although this study is of a relatively limited scale from an ecological perspective, we have successfully described viral diversity in samples collected from different sites, times and avian species and found evidence for differences across these factors. Finally, we demonstrate several potential applications of viral community analyses and anticipate a rapid expansion of viral ecology to move beyond the one‐host, one‐virus system and to consider both viruses and hosts as complex ecological communities.

AUTHOR CONTRIBUTIONS

The research was designed by M.W., M.K., A.C.H. and E.C.H. The research was performed by M.W. The data were analysed by M.W., J.S.E. and M.S. Reagents and analytical tools were contributed by M.K., A.C.H. and E.C.H. The manuscript was written by M.W., with input from all co‐authors. All authors gave final approval for publication.

Click here for additional data file.