https://orcid.org/0000-0003-2244-7412
ABSTRACT
The rumen is known to contain DNA-based viruses, although it is not known whether RNA-based viruses that infect fungi (mycoviruses) are also present. Analysis of publicly available rumen metatranscriptome sequence data from sheep rumen samples (n = 20) was used to assess whether RNA-based viruses exist within the ovine rumen. A total of 2466 unique RNA viral contigs were identified that had homology to nine viral families. The Partitiviridae was the most consistently observed mycoviral family. High variation in the abundance of each detected mycovirus suggests that rumen mycoviral populations vary greatly between individual sheep. Functional analysis of the genes within the assembled mycoviral contigs suggests that the mycoviruses detected had simple genomes, often only carrying the machinery required for replication. The fungal population of the ovine rumen was also assessed using metagenomics data from the same samples, and was consistently dominated by the phyla Ascomycota and Basidomycota. The strictly anaerobic phyla Neocallimastigomycota were also present in all samples but at a low abundance. This preliminary investigation has provided clear evidence that mycoviruses with RNA genomes exist in the rumen, with further in-depth studies now required to characterise this mycoviral community and determine its role in the rumen.
INTRODUCTION
It has been well established that DNA viruses (phages) of the viral order Caudovirales are endemic in the rumen, infecting and replicating within the bacterial populations responsible for the enteric fermentation of plant material (Gilbert and Klieve 2015). In contrast, very little is known about the viruses infecting rumen eukaryotes (i.e. protozoa and anaerobic fungi). To date, the most commonly found viruses infecting fungi (mycoviruses) have dsRNA, non-segmented genomes and icosahedral particles (Kotta-Loizou and Coutts 2017). Mycoviruses with ssDNA genomes, negative sense and positive sense ssRNA genomes and multi-segmented genomes have also been described (Ghabrial et al. 2015; Sato, Caston and Suzuki 2018). Several mycoviral families also include viruses that may be found in both fungi and plants (Roossinck 2019).
Mycoviruses tend to only encode a limited number of genes; for example, the genus Mitovirus (family Narnaviridae) co-exists with its fungal host and consists of a particle-free, ssRNA genome encoding only a single RNA-dependant RNA polymerase gene (Schoebel et al. 2018). Despite their relative genetic and structural simplicity, mycoviruses can significantly affect the viability of their fungal host. Research on mycoviruses has, therefore, focused on their ability to infect and reduce the pathogenicity of commercially important pathogenic plant fungi (Coetzee et al. 2010; Marzano et al. 2016). For example, mycoviruses of white mould (Sclerotinia sclerotiorum) can decrease the growth of the fungus, providing a potential biocontrol treatment (Ghabrial et al. 2015; Mu et al. 2018).
Anaerobic fungi (phylum Neocallimastigomycota) are important members of the rumen microbiome due to their key role in degrading fibre through the enzymatic and physical disruption of plant structural barriers (Gordon and Phillips 1998; Edwards et al. 2017; Rabee et al. 2019). As reviewed by Puniya et al. (2015), a wide range of studies have shown that this translates into clear benefits for ruminants, including improved feed efficiency, growth rates, voluntary feed intake and milk production. The presence of anaerobic fungi suggests that mycoviruses may also exist within the rumen, co-existing and potentially having a detrimental effect on the activity and abundance of ruminal anaerobic fungi.
Recent advances in RNA sequencing and transcriptomics have greatly facilitated the detection of mycoviruses in environmental samples (Marzano et al. 2016; Nerva et al. 2019), and whole viral genomes have been sequenced and taxonomically classified using this approach (Kondo et al. 2013; Svanella-Dumas et al. 2018). Therefore, the aim of this preliminary study was to determine if mycoviruses with RNA genomes were present within the rumen using a metatranscriptomic approach. This was done by utilising a published dataset (Shi et al. 2014) that was generated by metatranscriptome sequencing of RNA isolated from the rumen of 10 rams fed a lucerne pellet diet on two different days. Furthermore, a corresponding metagenome dataset generated from the same samples was used to assess the fungal populations present, giving insight into possible RNA virus–fungal host interactions.
MATERIALS AND METHODS
Metatranscriptome and metagenome datasets
The datasets utilised in this study were previously published by Shi et al. (2014). Therefore, only a brief background to the dataset is given here. Sequence data was generated from 10 rams kept at the Woodlands Research Station, which were fed a lucerne pellet diet (dietary chemical analysis detailed in Shi et al. (2014) in Table S7, Supporting Information). Rumen fluid samples were collected by stomach intubation 4 h after morning feeding, and immediately snap frozen in liquid nitrogen. Samples were collected from each animal on two different days, therefore providing a total of 20 samples. From each rumen sample, DNA for metagenomic analysis was extracted using a non-mechanical lysis DNA extraction method (Rosewarne et al. 2011). In addition, for each rumen fluid sample, RNA for metatranscriptome analysis was extracted using a hot lysis-acid phenol extraction method (Shi et al. 2014). Sequencing of RNA and DNA from the 20 samples was undertaken using an Illumina HiSeq 2000 platform (Illumina, San Diego, USA), and generated a total of 4.5 × 108 metatranscriptome reads (2 × 150 bp) and 3.4 × 109 metagenome reads (2 × 150 bp), which were downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database website (http://www.ncbi.nlm.nih.gov/sra/Accession no. SRA075938).
Sequence data assembly
Metagenome sequence data generated from DNA extracted from each of the 20 rumen fluid samples was combined into a single file for assembly using the Spherical pipeline (Hitch and Creevey 2018). This pipeline used Velvet (v1.2.10; Namiki et al. 2012) with a kmer size of 51 and the ‘-exp_cov auto’ command. Bowtie2 (v2.2.3, Langmead and Salzberg 2012) was used to align the raw reads to the assembly during the iterative process. Similarly, metatranscriptome data for each respective RNA sample was combined into a single file for assembly. The combined metatranscriptome dataset was assembled using Velvet with a kmer size of 51 and the ‘-exp_cov auto’ command was selected. Bowtie2 was also used to align the raw reads to the assembly during the iterative process.
Assembly annotation
The metatranscriptome assembly was annotated against the mycovirus subset of UNIPROT (downloaded June 2016) and the metagenome assembly against the fungal subset of UNIPROT (Bateman et al. 2015) using DIAMOND (v0.7.9; Buchfink, Xie and Huson 2014) with default settings and a bitscore cut-off of 40. The annotation files were then converted into general feature format (GFF) using the ‘blast2gff’ command in MGKIT (Rubino et al. 2014). The GFF files were then filtered to remove overlapping annotations using MGKITs ‘filter-gff overlap’ command with a maximum allowed overlap of 100 bp. The GFFs were edited to include the taxonomic and functional information (Gene Ontology (GO)) for each annotation using the ‘add-gff-info’ command of MGKIT.
Mycoviral taxonomy was assigned to contigs using the Viral Taxonomic Lineage Classification (VTLC) dataset (https://github.com/thh32/VTLC) and in the absence of any known rumen mycoviruses in the VTLC dataset, the names of the closest taxonomic matches were used to designate provisional taxonomy. A taxonomy assignment confirmation step was included prior to combining taxonomic counts for determination of relative viral abundance (Supplementary methods 1, Supporting Information). The confirmation step takes into account the known level of protein conservation between viruses of each family (based on the UNIPROT database annotations, minimum percentage identity of 45.13% for family level and 58% for virus level taxonomy). For example, contigs failing to reach the threshold percentage identity for classification at virus level were classified instead at family level.
Contigs that were initially annotated as being mycoviral in origin were also compared to the NCBI virus reference sequence database (update 18 March 2019) using blastn (threshold for homology e-value < 10−5) (BLAST+ v 2.2.31; Camacho et al. 2009).
Abundance calculations for viral contigs
Bowtie2 was used for the alignment of each sample's metagenome and metatranscriptome reads against their corresponding assembly. Each sample's alignment file was compressed from SAM format into BAM format using ‘samtools view –Sb’ (Samtools v1.2.1, Li et al. 2009). Samtools was then used to sort and index each sample using the ‘sort’ and ‘index’ commands.
Samtools idxstats was used to calculate the number of reads assigned to each viral annotated contig within each sample. The numbers of reads assigned to all contigs with the same taxonomic lineage were combined at both the virus and family levels. Taxonomic groups with fewer than 1000 reads across the 20 samples were filtered out of the study. The counts were then scaled to the sample with the fewest mycoviral-assigned reads (6398 reads) to prevent sequencing depth bias. The scaled values were multiplied by the average read length (200 bp), and divided by the total number of bases within all contigs assigned to the designated taxa. This provided a prediction of abundance for each taxonomic group identified within each sample. Taxonomic groups were filtered out of the dataset if not present within ≥5 samples with a minimum abundance of 100.
As well as taxonomy-based abundance analysis, HTSeq (v0.6; Anders, Pyl and Huber 2015) was used to identify the number of reads aligning (minimum alignment quality ≥ 8) to functionally annotated mycoviral genes in the metatranscriptome assembly. This was done for each metatranscriptome sample using the intersection-nonempty overlap resolution mode. The geneID option was used to ensure counts for multiple assembled copies of a gene were combined. Each GO term was required to have 200 reads aligning to its genes across the 20 samples to be included in the analysis.
Abundance calculations for fungal contigs
Bowtie2 was used for the alignment of each sample's metagenome reads against the metagenome assembly. Each sample's alignment file was compressed from SAM format into BAM format using ‘samtools view –Sb’. Samtools was then used to sort and index each sample using the ‘sort’ and ‘index’ commands.
Samtools idxstats was used to calculate the number of reads assigned to each fungal annotated contig within each sample. The numbers of reads assigned to all contigs with the same assigned lineage were combined at the phylum level, and within the Neocallimastigomycota phylum also to the genus level. Taxonomic groups with fewer than 1000 reads across the 20 samples were filtered out of the study. The counts were then scaled to the sample with the fewest fungal-assigned reads (1 892 175 reads) to prevent sequencing depth bias. The scaled values were multiplied by the average read length (200 bp) and divided by the total number of bases within all contigs assigned to the designated taxa. This provided a prediction of relative abundance for each taxonomic group identified within each sample. Taxonomic groups were filtered out of the dataset if not present within ≥5 samples with a minimum abundance of 100. All contigs with the same lineage were combined at the phylum level, and within the Neocallimastigomycota phylum, also to genus level.
Co-abundance network analysis
Co-abundance correlation analysis was undertaken using the abundances calculated for the 20 combined rumen fluid samples using SparCC (Friedman and Alm 2012) to provide pseudo P-values for each correlation. SparCC calculated the correlations over 1000 iterations and resampled 1000 times to produce bootstrap pseudo P-values. These correlations were then filtered to retain only significant results (pseudo P ≤ 0.05, correlation co-efficient ≥ 0.5). In this way, co-abundance analysis was undertaken to assess potential virus:virus interactions and fungal:virus interactions.
Data availability
The assembly, annotation and count files used in this study are provided at https://github.com/thh32/Rumen_mycoviruses. The viral contig sequences are also deposited at NCBI (Bioproject number PRJNA543390).
RESULTS
Mycovirus relative abundance
A total of 11 438 698 contigs, containing 2 129 415 091 bp, were assembled from the metatranscriptome. Of these, 2466 contigs, representing 1 673 296 bp, were assigned to a total of 30 different mycoviruses based on provisional taxonomic annotations. The contigs identified as being related to known mycoviruses accounted for 0.025% ± 0.008 of the combined metatranscriptome reads.
From the 30 mycoviruses, nine viral families were identified within the combined metatranscriptome dataset. The most abundant contigs were classified as belonging to the Partitiviridae family of dsRNA mycoviruses and plant viruses with an average abundance of 3000 per sample (Fig. 1A). Relative abundance analysis of the mycoviral population indicated that many of the 30 mycoviruses were present at low abundance. The majority (9 out of 10) of the most abundant mycoviruses were designated as dsRNA viruses (highlighted in Fig. 1B.), on the basis of homology to known dsRNA sequences present in the UNIPROT database. The variability of the mycoviral population between individual rumen samples was also tested, with the deviation of each sample from the average abundance being used to determine the variability of each taxonomic group identified (Fig. 1). All the mycoviruses were identified as having highly variable abundance across the samples, deviating more than 10% from their median relative abundance value.
Relative abundance in rumen samples of the viral contigs identified in the metatranscriptome assembly. Relative abundances are presented grouped at the level of family (A) and specific viruses (B) based on the closest match obtained during annotation of the contigs using UNIPROT. The relative abundance of each grouping within each rumen sample (n = 20) was used to generate boxplots that show the median, interquartile range and maximum and minimum abundances. The colouring of the family boxplot (A) is used in the specific virus boxplot (B) to indicate the associated family for the specific virus where it was known. Where it is not known, boxes are coloured based on viral genome type: black boxes for dsRNA genomes and grey boxes for ssRNA genomes. All dsRNA taxa in both graphs are indicated by an asterisk at the end of their name.
Comparison of the mycoviral contigs against the NCBI virus reference sequence database identified high-homology matches (>90% identity) for only 45 of the 2466 contigs. This suggests that only a small subset of the mycoviral contigs are related to the sequences of known viruses present in the current databases. The majority of these matches were to known viral families, including the Partitiviridae (15 matches), Alphaflexiviridae (17 matches) and Betaflexiviridae (6 matches), all of which are known to contain mycoviruses (Marzano et al. 2016). Many of the other viral families identified using the NCBI virus reference sequence database have an association with fungi. Potyviridae (1 match) can be transmitted via fungi (Wagh et al. 2016) whilst both Picobirnaviridae (3 matches) and Amalgaviridae (1 match) share close evolutionary histories with Partitivirdae (Malik et al. 2014; Krupovic, Dolja and Koonin 2015).
Functions
GO terms were used to functionally group the mycoviral genes within the assembly together, and to study their abundance in the different rumen samples. Seventeen different GO terms were identified (Fig. 2). Many of the identified GO terms are related to each other, due to the hierarchical nature of GO terms. For example, the proteins involved in viral capsid formation may be classified within three different GO terms. All identified GO terms focus on viral structure (GO:00 19028) and replication (GO:0 003968), particularly RNA-based replication as shown by the most abundant GO term being RNA-directed RNA polymerase activity.
Relative abundance in rumen samples of function gene ontology (GO) terms within the viral contigs identified using UNIPROT. Relative abundance within each rumen sample (n = 20) was defined as reads aligning to genes assigned to each GO term. The data was then used to generate a boxplot that shows the median, interquartile range and maximum and minimum abundances.
Fungal relative abundance
Across the 20 rumen samples, fungi accounted for 1.6% ± 0.6 of the combined metagenome reads. The presence of the anaerobic fungal phylum Neocallimastigomycota was confirmed in all the rumen samples, along with five other fungal phyla (Fig. 3). The predominant fungal phyla were identified as Ascomycota and Basidiomycota (Fig. S1, Supporting Information).
Relative abundance in rumen samples of fungal contigs identified from the metagenome assembly using UNIPROT. The relative abundance of each fungal contig was grouped at the phylum level for each sample (n = 20) and the data used to generate a boxplot showing the median, interquartile range and maximum and minimum abundance for each phylum.
Co-abundance network analysis
A co-abundance network (pseudo P ≤ 0.05, correlation ≥ 0.5), was developed to explore any potential associations between the different mycoviruses identified. Across the 30 mycoviruses, 25 significant correlations were detected, of which 56% (14) were positive and 44% (11) were negative (Fig. 4). Contigs that were most related to the Amasya cherry disease associated mycovirus were identified as the most highly connected, correlating with 11 other mycoviruses. This was followed by contigs most related to the Oyster mushroom isometric virus ii, which correlated with 8 other mycoviruses.
Co-abundance analysis of the mycoviruses identified within the metatranscriptome assembly. The size of each node is proportional to the number of significant correlations with the node. The interactions are colour-based—green for positive and red for negative correlations. The width of each interaction is proportional to the strength of its correlation. Only significant interactions are shown (P < 0.05).
Co-abundance networks between mycoviruses and fungi were also compared, with fungi classified to phylum level and the genera detected within the phylum Neocallimastigomycota (Tables 1 and S1, Supporting Information). The contigs most related to the Oyster mushroom spherical virus were found to be co-abundant with all of the fungal phyla within the dataset, except for Neocallimastigomycota. Within the Neocallimastigomycota phylum, Piromyces and Neocallimastix were identified at the genus level and accounted for 32.6 ± 5.9% and 37.2 ± 5.5%, respectively, of the metagenome reads assigned to the phylum. As with the phylum level analysis, neither of these genera correlated with any of the mycoviruses in the co-abundance analysis.
Significant correlations (pseudo P ≤ 0.05, correlation ≥ 0.5) identified from co-abundance analysis of detected mycoviruses and fungal phyla using SparCC.
| Mycovirus | Fungal phyla | Correlation | P-value |
|---|---|---|---|
| Oyster mushroom spherical virus | Ascomycota | 0.54 | 0.014 |
| Oyster mushroom spherical virus | Basidiomycota | 0.54 | 0.009 |
| Oyster mushroom spherical virus | Chytridiomycota | 0.54 | 0.015 |
| Oyster mushroom spherical virus | Glomeromycota | 0.58 | 0.010 |
| Oyster mushroom spherical virus | Microsporidia | 0.54 | 0.014 |
| Mycovirus | Fungal phyla | Correlation | P-value |
|---|---|---|---|
| Oyster mushroom spherical virus | Ascomycota | 0.54 | 0.014 |
| Oyster mushroom spherical virus | Basidiomycota | 0.54 | 0.009 |
| Oyster mushroom spherical virus | Chytridiomycota | 0.54 | 0.015 |
| Oyster mushroom spherical virus | Glomeromycota | 0.58 | 0.010 |
| Oyster mushroom spherical virus | Microsporidia | 0.54 | 0.014 |
Significant correlations (pseudo P ≤ 0.05, correlation ≥ 0.5) identified from co-abundance analysis of detected mycoviruses and fungal phyla using SparCC.
| Mycovirus | Fungal phyla | Correlation | P-value |
|---|---|---|---|
| Oyster mushroom spherical virus | Ascomycota | 0.54 | 0.014 |
| Oyster mushroom spherical virus | Basidiomycota | 0.54 | 0.009 |
| Oyster mushroom spherical virus | Chytridiomycota | 0.54 | 0.015 |
| Oyster mushroom spherical virus | Glomeromycota | 0.58 | 0.010 |
| Oyster mushroom spherical virus | Microsporidia | 0.54 | 0.014 |
| Mycovirus | Fungal phyla | Correlation | P-value |
|---|---|---|---|
| Oyster mushroom spherical virus | Ascomycota | 0.54 | 0.014 |
| Oyster mushroom spherical virus | Basidiomycota | 0.54 | 0.009 |
| Oyster mushroom spherical virus | Chytridiomycota | 0.54 | 0.015 |
| Oyster mushroom spherical virus | Glomeromycota | 0.58 | 0.010 |
| Oyster mushroom spherical virus | Microsporidia | 0.54 | 0.014 |
DISCUSSION
This study provides the first description of mycoviruses within the ovine rumen. Using metatranscriptome sequencing, 30 mycoviruses, belonging to nine viral families, were detected. Overall, the mycoviral populations were of very low abundance, with homologous sequences accounting for only 0.025% of the metatranscriptome. The mycoviral population was also found to be highly variable (>10%) between individual samples, despite the animals being similarly managed and fed. The high variability observed could suggest that the abundance of this community is highly variable between individual sheep. Alternatively, the non-targeted sequencing approach used within this study may have resulted in sequencing depths that limited the capacity to fully capture the mycoviral populations present, particularly for individuals where mycoviruses were present at very low concentrations. The detection of a relatively small proportion of mycovirus-related sequences within the metatranscriptome suggests future studies should physically concentrate or enrich this viral community. This approach has been undertaken in previous studies of rumen DNA virus populations (Berg Miller et al. 2012; Anderson, Sullivan and Fernando 2017). However, studies of RNA viruses will require additional steps to purify and transcribe the RNA and construct sequence libraries prior to sequencing (Marzano et al. 2016; Mu et al. 2018).
Of the top 10 most abundant mycoviruses detected, 8 were identified as having dsRNA genomes. Mycoviruses with dsRNA genomes are more commonly found in nature (Kotta-Loizou and Coutts 2017). In accordance with this, the Partitiviridae family of dsRNA mycoviruses was found to be the most abundant mycovirus type present in the ovine rumen. Co-abundance analysis indicated that the mycoviruses were not associated with Neocallimastigomycota within the rumen. Therefore, it is speculated that the dsRNA mycoviruses detected in the rumen samples may have been mainly associated with the sac fungi (Ascomycota and Basidiomycota), which were the most dominant fungal phyla detected in the metagenome (Fig. S1, Supporting Information). These fungi were also previously reported to be detected in the DNA-based analysis of fungi present in the bovine rumen (Vaidya et al. 2018). These phyla have been shown to be more dominant in the rumen of animals fed high-grain diets (Ishaq et al. 2017), which was not the diet used in the Shi et al. (2014) study.
Whether the sac fungi are normal commensals of the rumen and can contribute to rumen fermentation remains to be demonstrated, although dietary supplementation with Aspergillus oryzae extracts and yeast have been reported to have beneficial effects on ruminal fermentation (Puniya et al. 2015). However, eukaryotic metatranscriptome analysis of the rumen did not find transcripts associated with the Ascomycota and Basidiomycota phyla, only the Neocallimastigomycota (Qi et al. 2011). As such, we consider that it is only the strict anaerobic fungi of the Neocallimastigomycota phylum that are likely to actively proliferate within the rumen. Based on the current study, it is speculated that the Ascomycota and Basidiomycota were potentially introduced into the rumen via the feed and/or environment (soil, water). If this is true, it is not clear how long they would remain viable under ruminal conditions. Future studies examining mycovirus populations in the rumen should, therefore, also include the analysis of feed and environmental samples.
Strictly anaerobic fungi, whilst consistently present in all the samples, represented only 1.9% (± 0.8) of the fungal sequences identified within the metagenome. This low relative abundance is likely to be associated with the rumen fluid sampling approach, as anaerobic fungi are only transiently present in rumen liquid as zoospores (Gruninger et al. 2014). Previous studies have observed that in animals fed daily, the timing of peak zoospore densities in the rumen fluid depends on the genus, with Piromyces and Caecomyces zoospores peaking 1 h after feeding (Orpin 1977) and Neocallimastix after approximately 30 min (Orpin 1976). With the sheep used to generate the dataset employed in this study, feed was given twice daily at 9 am and then 6 h later at 3 pm and the samples were collected 4 h post feeding (Shi et al. 2014). We therefore suggest that in future studies the mycoviral portion of the rumen microbiome be sampled from both the rumen liquid and solid fractions, following a time course directly after feeding.
Functional analysis of the identified mycoviral contigs showed a high abundance of gene ontology terms with roles in either viral particle structure or genome replication. This is likely to be due to the small genome size and functional simplicity of RNA viruses preventing the acquisition of non-essential genes (Hou and Lin 2009). Furthermore, the absence of any previously published rumen-associated mycoviruses in sequence databases limited the taxonomic assignment of the assembled contigs. Therefore, the mycoviral taxonomic annotations used within this study should be interpreted cautiously until the rumen-associated mycoviruses are physically isolated, characterised and taxonomically classified.
Co-abundance network analysis was conducted using the relative abundances of the mycoviruses, and 25 significant interactions were identified. Further studies are needed to clarify the basis of these interactions. From the high number of positive correlations (56%) linking mycoviruses, it might be the case that multiple mycoviruses infect the same fungal host species. Conversely, the identified negative correlations may suggest that, as observed for virus:host interactions (Arribas et al. 2018), mycoviruses may use strategies to inhibit their competitors in order to improve their ability to infect the host fungi.
Co-abundance network analysis using the relative abundances of the mycoviruses and fungi indicated that only the fungal phyla that contained no strict anaerobic fungi had significant interactions with the mycoviruses. As mycoviruses are known to influence fungal populations associated with plants and other environments (Ghabrial et al. 2015; Son, Yu and Kim 2015), it could be speculated that the rumen-associated mycoviruses may modulate these other fungal populations rather than the beneficial plant-degrading strict anaerobic fungi. However, due to the relatively low abundance of both fungal and mycoviral contigs in the datasets, the findings of this analysis should be interpreted cautiously.
While this study has provided new insights into the presence and types of mycoviruses present in the ovine rumen, it has also raised many questions and highlighted the need for further studies investigating them. Only following further detailed studies can the impact that these mycoviruses have on rumen fungal populations and function be determined.
ACKNOWLEDGEMENTS
We thank the authors of Shi et al. (2014) for providing access to their dataset, Professor Chris Creevey for his advice and mentorship during this work and Dr Pat Blackall for a critical reading of the manuscript.
Conflicts of interest
None declared.
