Molecular Detection and Characterization of Mycoviruses in the Genus Armillaria in the Czech Republic

1. Introduction

Fungal pathogens are important members in forest ecosystems. They substantially affect diversity, structure and dynamics of forest communities [1,2]. The fungal genus Armillaria includes more than 40 described species [3], which are causal agents of root rot in a wide variety of trees, shrubs and some herbs [4], including economically important conifers (eg. Abies, Picea, Pinus) and agricultural crops (eg. Citrus, Juglans, Malus, Prunus, Vitis) [5]. Norway spruce (Picea abies) is a major host of Armillaria in the Czech Republic. Over the last years, Czech forests of Norway spruce seem to be suffering a generalized decline [6,7], which is enhanced by pathogenic fungi like Armillaria [8,9]. The decline is often caused by a combination of abiotic and biotic factors that limit tree growth, reduce foliage quality, and weaken root systems. These interactions can kill individual trees and entire stands [10,11]. Although there are relatively efficient methods for protection of agricultural crops against Armillaria root rod [12]. None of the established methods can completely eliminate the mycelium from the contaminated site [13], therefore other strategies have to be pushed towards the development of environmentally friendly control approaches [14].

Fungal viruses (mycoviruses) of plant-interacting fungi are particularly significant for crop protection because they can influence the phenotype of their host. An increasing number of mycoviruses with the ability to induce hypovirulence in their host were being reported [15,16,17]. There are important forest pathogenic fungi, such as Chryphonectria parasitica, Phytophthora spp., Ophiostoma spp., Gremmeniella abietina, Hymenoscyphus fraxineus, Heterobasidion annosum and Armillaria spp., which host mycoviruses with diverse fungus-virus relationships [18,19,20]. In Armillaria, early studies by Blattny (1973) [21] and Reaves et al. (1988) [22] described the presence of virus-like particles in A. mellea and A. ostoyae, but it was not until 2021, when viruses were confirmed by molecular characterization [13,23,24].

The main objectives of this study were (i) to confirm the occurrence of mycoviruses in a collection of isolates of the genus Armillaria in Central Europe, (ii) to describe their genomic features and phylogenetic relationships.

2. Materials and Methods

Sampling

Armillaria isolates were collected in the years 2019 – 2020 in the Czech Republic, mainly northeastern Moravia, the Moravian Karst and southeastern Bohemia (Table 1). Spores from fruiting bodies or rhizomorphs were transferred and cultured on ME agar (7.5 g bacteriological agar, 3.5 g corn agar, 10 g malt extract, 10 g glucose / 0.5 l; pH 5.5) at 25 °C for aproximately four weeks. A total of 13 samples, include Armillaria ostoye and cepistipes, were cultivated and analyzed.

Isolation of Double-Stranded (ds) RNA

Samples were first examined for the presence of potential viral dsRNA molecules. Mycelium from one Petri dish was lyophilized and collected in a 50 ml Falcon tube. Two stainless steel beads with a diameter of 10 mm and liquid nitrogen were added to the Falcon tube and mycellium was homogenized by vortexing. After homogenization, dsRNA was extracted using protocol [25] with some modifications as desribed in Tonka et al. (2021) [26]. Previously confirmed dsRNA-hosting isolates of Gremmeniella [27] and Phytophthora [28] were used as positive control of each isolation.

Isolation of Total RNA

Thirteen samples were selected for stranded total RNA-seq analysis (Table 1).

Approximately 50 - 100 mg of fresh mycelium was collected from agar from each sample and transferred to the 2 ml tubes with steel beads. The tubes were immersed in liquid nitrogen and vortexed until the mycelium was ground to a fine powder. Isolation of total RNA was performed using a SPLIT RNA Extraction Kit (Lexogen, Vienna, Austria) following the protocol provided by the manufacture, eluted in 30 µl of EB buffer. Total RNA was visualized on a agarose gel and stored at -80 °C.

Stranded Total RNA Sequencing

Isolated total RNA of 13 Armillaria samples was pooled and treated with TURBO DNA-free TM Kit (InvitrogenTM, Thermo Fisher Scientific). RNA quantity was measured by Qubit^® Fluorometer (Life Technologies, Thermo Fisher Scientific). Total RNA was sent to SEQme s.r.o. (Dobris, CR) for RNA library construction and RNA sequencing. Prior to the library preparation, ribosomal RNA (rRNA) was depleted using the NEBNext rRNA Depletion Kit (Human/Mouse/Rat). The library was prepared using an NEBNext Ultra II Directional RNA Library Prep Kit for Illumina. The quality of prepared library was checked using Agilent Bioanalyzer 2100 High sensitivity DNA Kit, Invitrogen Collibri Library Quantification Kit and Qubit 1X dsDNA High-Sensitivity Assay Kit. A KAPA Library Quantification Kit for Illumina platform was used for absolute, qPCR-based quantification of the Illumina libraries flanked by the P5 and P7 flow cell oligo sequences. Libraries underwent paired-end (PE) (2 × 150 nt) sequencing on a NovaSeq6000 (DS-150) (Illumina, San Diego, CA, USA) using a NovaSeq S4 v1.5 reagent kit. An “in-lane” PhiX control spike was included in each lane of the flow cell.

Bioinformatics

Data Preprocessing

The raw sequencing data was downloaded from the data repository of the sequencing company SEQme s. r. o. (Dobris, CR). The processing was carried out on the local server of the University of South Bohemia. In the first step, data quality was assessed using the FASTQC v.0.11.9 program [29]. This control step revealed the presence of a small amount of adapter sequences. The nucleotide sequences of adapters were obtained according to the used library kit and sequencer machine mentioned earlier [30]. Adapter trimming and quality trimming was performed using the program Cutadapt v3.4 [31] with requirements for a Phred score higher than 30 and a minimum length of the truncated sequence of 50 bases. The software is available at the code depository GitHub [32]. Subsequently, data quality of trimmed reads was again assessed in FASTQC program mentioned previously.

Host Reads Removal

After preprocessing data, it was necessary to remove reads that are belong of host fungi. The STAR v2.7.9a program was used for this purpose [33]. Assembly GCA_900157425.1 version 2 of strain C18/9 of A. ostoyae was used as a reference sequence. After mapping the reading frames to the Armillaria genome, only unmapped reads were left for further processing.

Discovery of Known Virus

To reveal already known viruses, Viral NCBI database [34] were used for this purpose. We downloaded only complete viral refseq nucleotide sequences. Our reads were aligned to this reference database by BWA v0.7.17.-r1188 program package [35]. In the next step we calculated the coverage for each viral genome using samtools kit v1.16.1 [32]. If the coverage is higher then 80 percent, we visualized the alignment for particular viral genome with IGV program v2.16.1 [36] and manually confirm or denied the presence of the published virus.

Discovery of Novel Viruses

In the first step De Novo assembly was performer with unmapped reads from STAR aligner. The toolkit SPAdes version 3.15.3 [37] with default settings for metagenomics was used for this purpose [32]. Assembled contigs shorter than 1000 bp were discarded. In the next step, the designed contigs were compared to several databases:

Viral UniProtKB [38], Viral NCBI [34], RVDB [39], Virus-Host DB [40]. Each of the database was downloaded to local server and were used to search for similarities between database sequences and assembled contigs. For this purpose, program BLASTx and BLASTn v2.12.0+ was used. The threshold for evalue was set to 1e-3. All contigs with at least one hit were further used to search for similarity with the NCBI database nr (RefSeq non-redundant proteins) or nt (RefSeq Nucleotide) with the specified Armillaria taxon. In the event that a positive hit was found, the evalue values were subsequently compared. If this value was higher for the results from viral databases, the contig was removed. In the last step, the remaining contigs were searched with the entire nr and nt databases and threshold for evalue was set 1e-5.

Potential protein encoding segments were detected with coding open reading frame (ORF) finder using Geneious^® v. 8.1.9. Depth of Coverage: For the calculation of the coverage depth, we used the following formula: (Total reads mapped to the final identified virus * average read length)/virus genome or contig length.

Ribozymes Detection

The identification of potential ribozymes followed a customized protocol adapted from the viroid-like sequence search pipeline outlined by Lee et al. in 2023 [41].The process began with the recognition of covalently closed circular RNAs (cccRNAs) using a modified version of the reference-free CIRIT algorithm, as introduced by Qin et al. in 2020 [42]. This algorithm entails the exploration of overlapping regions between the starting and ending points of contigs to pinpoint cccRNAs. Unlike the original algorithm, our modified version repeatedly attempted to dissect potential cccRNAs into individual units, with the condition of maintaining a minimum 95% similarity within repeated regions. Subsequently, the identified cccRNAs underwent a search for known self-cleaving ribozymes, employing Infernal, a tool introduced by Nawrocki and Eddy in 2013 [43]. Ribozymes surpassing Rfam’s curated gathering threshold or exhibiting E-values below 0.1 were identified in each polarity. In addition to this, a subset of sequences underwent further scrutiny using RNAmotif, a tool developed by Macke et al. in 2001 [44]. This supplementary analysis aimed to uncover more diverse ribozymes that may not have been initially detected by Infernal.

Retro Transcription (RT) PCR and Sanger Sequencing

Aproximately 50-100 mg of fresh mycelium was collected a homogenized in mortar. Crushed mycelium was mixed with 600 µl of lysis buffer (LB) and 6 µl 2-mercaptoethanol and transffered to the 1.5 ml sterile tube. Isolation of total RNA was proceed using PureLinkTM RNA Mini Kit (Invitrogen, USA). The RNA was eluted in 30 μl of elution buffer and 1 µl of RNAase inhibitor was added. The quality of total RNA was assessed by gel electrophoresis in 1.5% agarose gel and total RNA was stored in -80°C. LunaScript^® RT SuperMix Kit (New England Biolabs, UK) was used for synthesis of cDNA as by Tonka (2022) [45]. The success of cDNA synthesis was verified by amplification of the eukaryotic regiontranslation elongation factor 1-alpha – eEF1A (tefa) and checked by gel electrophoresis in 1.5% agarose gel [46]. If amplification was successful, the cDNA was used in a PCR reaction with virus-specific primers. The virus-specific primers were designed to partially amplify the ORFs of putative viruses base on sequencing date from RNA-seq using program Geneious v. 8.1.9 PCR amplification was performed with 12 μl PPP Master Mix (Top-Bio,Vestec, ČR) 1 μl of each 10 mM primer (Table S1), 4 μl of cDNA and PCR grade water in a total volume of 25 μl. PCR conditions were: 94°C 2 min, 25x (94°C 1 min., 60°C 1 min., 72°C 2 min.) 72°C 5 min. Due to the low concentration of putative viruses, the PCR amplification was repeated with fresh chemicals and as a template was used 4 μl of resulting product from the first PCR amplification, then PCR products were visualized by gel electrophoresis in 1.5% agarose gel. All resulting fragments with the appropriate size were cutted from the gel, cleaned using NucleoSpin^® Gel and PCR Clean-up kit-Macherey-Nagel (BioTech a.s, ČR), and ExoSAP-IT™ PCR Product Cleanup Reagent (Thermo Fisher Scientific, USA), then fragments were sent to SEQme (Dobříš, ČR) for Sanger sequencing.

Conserved Domains

NCBI CD-search tool [47] was used for search of putative conserved domains. In order to look for putative conserved motifs of RdRp regions of putative viruses, predicted amino acid (aa) sequences were aligned using MUSCLE 3.8.425 in Geneious v. 8.1.9 to aa sequences of related viruses retrieved from the GeneBank.

Phylogenetic Analyses

Amino acid (aa) sequences of RdRp regions of each virus were included in phylogenetic analyses. Sequences were aligned in Geneious v. 8.1.9 by MUSCLE 3.8.425 [48] together with known aa sequences of viruses considered to be related. Phylogenetic trees were built using the maximum likelihood method [49] in RAxML-HPC v.8 on XSEDE running in the CIPRES Science Gateway web portal [50]. Bootstrapping was performed by using the recommended parameters provided by the CIPRES Science Gateway portal. The trees were visualized in FIGTREE (V1.4.4).

3. Results and Discussion

In this study, we described the genomes of five putative single stranded (ss) RNA viruses hosted by pathogenic fungi belonging to the genus Armillaria from Czech Republic. These new viruses include member of family Tymoviridae and members of recently described virus group, tentatively named “Ambiviruses” [51,52,53].

dsRNA Screening

No dsRNA elements were detected in Czech isolates of Armillaria spp. whereas the corresponding dsRNA patterns of the positive controls were detected in every assay. As later demonstrated by high-throughput sequencing (HTS) of total RNA, no dsRNA viruses appear to occur in our collection of isolates. This result agrees with the study Linnakoski et al. (2021) [13], that dealt with the detection of mycoviruses in isolates of fungi of the genus Armillaria from Finland, Russia and North Africa and it also corresponds to the results of Dvořák in 2008 [54], which did not detect dsRNA viruses in Armillaria from the Czech Republic. However, Shamsi et al. (2004) in their work report two detected mycoviruses with a dsRNA genome (a partitivirus and a phlegivirus) in samples originating from Switzerland [24]. Dvořák in 2008 [55] and Blattný et al. in 1973 [21] studied the Armillaria fungi virome from the Czech Republic in their works. In his research, Dvořák (2008) [54] detected possible double stranded molecules in forty samples from the Czech Republic, but they were not confirmed to be viral, which coincides with the results of this study. Likewise, Blattný et al. (1973) [21] and Reaves et al. (1988) [22] detected virus-like particles in their works on Armillaria mellea and Armillaria ostoye, but they were no further verified molecularly. Blattný et al. (1973) [21] described these particles as rod-shaped (22–28×119 nm) or isometric (30 nm).

Identification of Final ssRNA Viral Sequences

Total RNA-sequencing generated 113 million paired-end (PE) reads. After quality trimming and de novo assembly a total of 8809 contigs longer than 1 kb were obtained. BLASTX comparison of the contigs revealed four viral contigs with sequence similarities to members of the recently described circular RNA virus group “Ambiviruses” and one related to members of order Tymovirales, in particular, the family Tymoviridae (Table 2).

Ambi-like viruses. Sutela et al. In 2020 [51] was the first study describing “ambiviruses” in the endomycorrhizal fungi Ceratobasidium sp. and Tulasnella sp. Then, they were discovered in Cryphonectria parasitica [52] and many agaricomycetes including Armillaria spp. [13,24], Heterobasidion spp. [20], Rhizoctonia spp. [13], and Phlebiopsis gigantea [55]. As in these studies, the genomic sequences of the putative ambi-like viruses detected in Czech isolates of Armillaria spp. were circa 4.5 kb long and contained typical RdRP conserved domains including the GDD motif (Figure S3A). Based on the pairwise sequence comparison (PASC) percentages (Table S2A, S2B) four different viruses were identified and not considered as variants of the same viruses following the same criteria as Sutela et al. in 2020 [51] and Forgia at al. in 2021 [52].

They were tentatively named considering their host identity (Table 2). AALV1, AoALV2 and AoALV3 contained two and AoALV4 contained three non-overlapping ambisense ORFs. AALV1 encodes a first positive (+) sense ORF (ORF A) corresponding to the RdRP (713 aa long) and a second negative (-) sense ORF (ORF B) corresponding to a hypothetical protein (HP) 410 aa long. AoALV2 encodes the (-) sense ORF B-HP 655 aa long, and the (+) sense ORF A-RdRP 710 aa long. AoALV3 contains first (+) sense ORF corresponding to a hypothetical protein (ORF B) 400 aa long and a second negative sense ORF (ORF A) corresponding to RdRp 682 aa long. This disposition shows that this contig was assembled in the 3′ – 5′ sense. AoALV4 encodes a first (-) sense ORF B-HP 405 aa long, a second negative ORF 220 aa long, and a third (+) sense ORF A-RdRP 713 aa long.

Ambivir genomes have the unique feature of having circular genomes encoding RdRPs and divergent ribozymes in various combinations in both sense and antisense orientation [53]. In our analyses, several hammerhead (HHRz) and hairpin (HpRz) ribozymes were detected. This is in agreement with Forgia et al. in 2023 [53], which describes HHRz and HpRz motifs and the predicted cleavage-sites in fungal known ambiviruses, but also in ambiviricot RNA-directed RNA polymerase (RdRp) palmprints of 439 distinct species-like operational taxonomic units (sOTUs) encountered in the GenBank. Although HHRz typically consist of 30 to 40 nucleotides, in Armillaria ambiviruses some of them seem to be a bit longer ranging from 55 to 78 nts. The representation of Armillaria ambiviruses and their ribozymes is illustrated in Figure 1a–d, and the positions in the particular ambivirus contigs are shown in Supplementary Table S3.

The phylogenetic relationships of Czech A. ostoyae ambi-like viruses with other “ambiviruses” (Figure 2) show that they group in two separated virus clusters, both of them including viruses from A. mellea and A. borealis. This result suggests a monophyletic origin of the ambiviricot sequences in Armillaria spp.

Armillaria tymovirus 1. A viral contig resembling features of those from the genus Tymoviridae was detected in Armillaria ostoyae. Armillaria ostoyae tymovirus 1 is circa 6.8 kb length and encodes one (+) sense ORF with 2172 aa. Based on sequence analysis of tymoviruses, conserved sequences were revealed (Figure 3), include conserved regions belonging to pfam01660 super family Vmethyltransf cl03298 in region 143-314 (e-value 6.72e-09). This methyltransferase domain has been detected in a wide range of hosts and is involved in mRNA capping. Furthermore, conserved sequences were detected belonging to pfam01443 superfamily Viral_helicase1 cl26263 located 1113-1322 (e-value 2.66e-14) and ps-ssRNAv_RdRp-like super family cl40470 located 1600-1799 (e-value 2.73e-25) [56].

The order Tymovirales was first described in 2004 and currently comprises five families: Alphaflexiviridae, Betaflexiviridae, Gammaflexiviridae, Deltaflexiviridae and Tymoviridae [57]. Members of the order Tymovirales have a 5.9 to 9.0 kb (+) ssRNA genome that is often polyadenylated. The largest protein, the replication-associated polyprotein (RP) is encoded by all members of the Tymovirales order. Tymovirales RP usually contains sets of conserved functional domains [58].

Viruses belonging to this order usually have a wide range of host organisms. Betaflexiviridae and Tymoviridae are usually considered plant viruses [58], a single member of Gammaflexiviridae has been detected in a filamentous fungus [57,59], and Alphaflexiviridae have been detected both in plants and fungi [60]. Although most viruses in the order Tymovirales are plant viruses, several members of this order are known to infect plant pathogenic fungi, including Botrytis cinerea, Fusarium boothii, Fusarium graminearum, R. solani, and S. Sclerotiorum [61,62,63,64,65]. Armillaria ostoyae tymovirus 1 is the second tymovirus described within the genus Armillaria after the study Shami et al. (2004) [24].

The phylogenetic relationships with other members from order Tymoviridae from GenBank are shown in Figure 3. The phylogenetically closest virus to the virus AoTV1 is Lentinula edodes tymo-like virus 1, which was detected in a Chinese sample of the fungus Lentinula edodes.

RT-PCR Screening

Some of the viruses described in Czech isolates of Armillaria (Figure 1) are hosted by more than one isolate, and, even by different species, suggesting interspecies transmission. AALV1 is present in isolate 1 (Armillaria cepistipes) and in isolate 7 (A. ostoye), AoALV2 was detected just in isolate 7 (A. ostoye), AoALV3 is present in isolates 7 and 13 (A. ostoye), AoALV4 in isolates 6, 7, 8 and 13 (A. ostoye), and AoTV was detected just in isolate 13 (A. ostoye). Interesting, isolates 1 and 7, between which virus interspecies transmission has occurred, were collected at locoliths approximately 50 km apart. These results suggest that ambiviruses are transmitted efficently in Armillaria in the Czech Republic, as it seems to happen in other regions of the Northern hemisphere [13,24], where ambi-like viruses have been found to be very common in Armillaria isolates and occur often in A. borealis and mellea from the Finland, Siberia and Switzerland. This efficiency could also be related to the fact that fungi of the genus Armillaria, but also other fungi causing tree root rot, such as Heterobasidion spp. or Rosellinia necatrix, usually form large clonal individuals that grow for decades. It can therefore be considered that virus accumulation and transmission. Interspecies transmission of mycoviruses is considered rather rare. However, laboratory experiments have shown that it is possible within the fungal genera Aspergillus [66], Sclerotinia [67] and Cryphonectria [68]. In vitro studies have shown mycovirus transmissions between somatically incompatible fungal strains of Heterobasidion [69,70,71], being relatively common in both the laboratory and in nature.

5. Conclusions

This study describes the tymovirus in an Armilaria sp. and confirms the occurrence of “ambiviruses” in Central European populations of A. ostoyae and A. cepistipes. The effect of the detected viruses on infected Armillaria hosts and their ability to be transmitted intra- and interspecies should be further investigated.