Viral Interactions and Pathogenesis during Multiple Viral Infections in Agaricus bisporus

ABSTRACT Viral interactions during multiple viral infections were examined in Agaricus bisporus cultures harboring 9 viruses (comprising 18 distinct viral RNAs) by statistically analyzing their relative abundance in fruitbodies. Four clusters of viral RNA were identified that suggested synergism and coreplication. Pairwise correlations revealed negative and positive correlations between clusters, indicating further synergisms and an antagonism involving a group containing a putative hypovirus and four nonhost ORFan RNAs (RNAs with no similarity to known sequences) possibly acting as defective interfering RNAs. The disease phenotype was observed in 10 to 15% of the fruitbodies apparently randomly located among asymptomatic fruitbodies. The degree of symptom expression consistently correlated with the levels of the multipartite virus AbV16. Diseased fruitbodies contained very high levels of AbV16 and AbV6 RNA2; these levels were orders of magnitude higher than those in asymptomatic tissues and were shown statistically to be discretely higher populations of abundance, indicating an exponential shift in the replicative capacity of the virus. High levels of AbV16 replication were specific to the fruitbody and not found in the underlying mycelium. There appeared to be a stochastic element occurring in these viral interactions, as observed in the distribution of diseased symptoms across a culture, differences in variance between experiments, and a number of additional viruses undergoing the step-jump in levels between experiments. Possible mechanisms for these multiple and simultaneous viral interactions in single culture are discussed in relation to known virus-host regulatory mechanisms for viral replication and whether additional factors could be considered to account for the 1,000-fold increase in AbV16 and AbV6 RNA2 levels.

D uring viral infection, numerous and complex virus-host interactions occur, including appropriation of host machinery and induction of host antiviral defense mechanisms. Further interactions, virus-virus, can result from double and triple virus infections, such as synergistic effects on disease symptom development (1)(2)(3)(4). However, the complexity of interactions is massively increased in multiple viral infections where a host harbors large numbers of viruses, a condition that appears to be widespread in some perennial plants and fungi. For example, the grapevine Vitis vinifera harbors large numbers of mycoviruses and plant viruses (5,6), and there are numerous reports of multiple viral infections in fungi, such as 21 virus-like sequences in a hypovirulent strain of the plant pathogen Sclerotium rolfsii (7), 11 viruses in a single isolate of Beauveria bassiana (8), 16 diverse viruses in Fusarium poae (9), and 26 distinct viral RNAs in a single isolate of Agaricus bisporus (10). Furthermore, deep sequencing technology and bioinformatic analyses of publicly available transcript databases are likely to further increase the number of viral observations in fungi (11).
Multiple viral infections are often asymptomatic or benign or at least have a low fitness cost to the host. However, specific viruses within mixed infections have been associated with a phenotype as pathogenic, changing host behavior by reducing a host's virulence or stabilizing a coinfection, for instance, in Cryphonectria parasitica (4,12). Some plant diseases occur only upon coinfection, while singular infections are relatively benign, e.g., rice tungro (13), maize lethal necrosis (14), and cowpea stunt disease (15). Virus levels within multiple virus infections can be dynamic over time with variable viral profiles influenced by host, virus-virus, or environment factors, which all can affect pathogenicity (3,16). Synergistic and antagonistic viral interactions have been associated with changes in viral pathogenesis and their evolution and spatial separation in plant tissues. The outcome of multiple infections can be difficult to predict (17)(18)(19), for instance, with the activation or suppression of RNA silencing and the formation of aggregated multiple viral genome structures producing collective infectious units, both having the potential for synergistic and antagonistic virus-virus interactions (20,21).
The fungus Agaricus bisporus is cultivated commercially to produce white button mushroom fruitbodies, a high-value crop. Thirty-nine nonhost RNAs of probable viral origin have been identified in A. bisporus. Three viral diseases of A. bisporus are economically important due to yield loss or quality loss. La France disease has symptoms of slow mycelial growth and distorted, delayed, or massively reduced fruitbody emergence and is caused by the particulate virus AbV1, consisting of 6 RNAs and 3 associated RNAs (22). Brown cap mushroom disease (BCMD) has symptoms of uniform to variable fruitbody cap browning, a significant quality-loss feature (23), and mushroom virus X (MVX), or Patch disease, has symptoms similar to those of La France disease, with massive yield loss (as nonproductive patches) and distorted mushrooms but the absence of AbV1 virus (24). BCMD and Patch disease were first observed and described at a similar time in the 1990s but often at different locations. The name MVX was originally used to describe both of these emerging diseases, but they are now considered distinct. Fruitbodies of both diseases contain multiple viral infections, and 30 RNAs have been sequenced, consisting of 18 viruses and 8 nonhost ORFans (RNAs with no similarity to known sequences) (10). Sequence comparison suggests that many of these viruses and ORFans are common for both diseases despite initially disparate identification based on gel separation. Gel-separated band H1 is equivalent to the hypovirus-like AbV2, band H3 is equivalent to ORFan2, and bands 18,19,22, and 23 are equivalent to AbV16 RNAs 2, 1, 3, and 4, respectively (10,24). Mushroom bacilliform virus (MBV) has been characterized and identified in the deep-sequenced RNAs of mushrooms with BCMD symptoms and also occurs as a coinfection with AbV1 in La France disease (10,25).
Previous investigations into virus-virus interactions have largely involved binary or lower-order interactions and have not considered multiple virus infections. It is not known how large numbers of viruses in a tissue interact and compete for the host cells' resources and how adaptive strategies for low-level viral coexistence adopted by some member viruses are affected by the high replicative activity of others. In this study, viruses from an A. bisporus strain showing pathogenic symptoms were introduced to an established mycelium of A. bisporus harboring ubiquitous viruses, and the subsequent viral interactions were observed through the analysis of viral abundance levels. A variety of statistical analysis and modeling approaches were used to examine virus-virus interactions in fruitbodies displaying disease symptoms and in underlying mycelium. This is the first time the interactions among such a large number of coinfecting viruses have been studied.

RESULTS
Fully colonized mycelia of Agaricus bisporus strain (A15), growing on compost and harboring ubiquitous viruses and ORFans, were infected with additional viruses and ORFans from A. bisporus strain MVX-4569. This resulted in the production of brown symptomatic mushrooms (color values, 9.5 to 20.8) among a prevalence of white asymptomatic mushrooms (color values, 7.7 to 8.6). A preliminary experiment was carried out to determine how many of the 30 viral RNAs identified by Deakin et al. (10) were present in MVX-4569-infected mushrooms using quantitative PCR (qPCR) of a single brown mushroom and a single white mushroom from each experiment. This revealed the presence of 18 viral RNAs, comprising 9 distinct viruses, including the multipartite AbV6 (2 segments) and AbV16 (4 segments) and 5 nonhost RNAs, which, per Deakin et al. (10), we term ORFans. The abundances of these 18 viral RNAs were measured by qPCR in the mushroom samples (40 samples from experiment 1 and 49 from experiment 2). Viral clusters. Viral RNAs were clustered on the basis of similar relative abundancies in the 89 samples using hierarchical and k-medoids clustering techniques (Fig. 1). Four clusters were identified by the k-medoids technique, determined by both gap statistics and within sum of squares (see Fig. S1 in the supplemental material). Both hierarchical and k-medoids methods produced clusters with very similar members (Table 1), with only AbV2 differing between the two. A silhouette plot showed there was very little difference in overall cluster scores between the two clustering methods for AbV2 (Fig. S2). Cluster 1 consists of ORFans 2, 3, 5, and 7 and, more loosely, AbV2; cluster 2 contains AbSV, AbV10, AbV12, AbV6 RNA1, and AbV6 RNA2; cluster 3 consists of the four components of AbV16 and ORFan8 (statistically confirming previously described associations [10,23,24]); and cluster 4 consists of AbV14 and AbV9. MBV was located in cluster 1 using hierarchical and k-medoids clustering techniques.
Relationships within and between viral clusters: correlation and interactions. Pairwise correlation of the abundances of all RNAs showed positive correlations within clusters and both positive and negative correlations between clusters ( Fig. 2 and Fig. S3). Interaction between the groups was apparent; for instance, cluster 2 and cluster 3 had strong positive correlations among their member RNAs, while cluster 1 and cluster 4 member RNAs were negatively correlated (Fig. 2). MBV showed moderate correlation to cluster 3.
Population variance of viral clusters. Striking differences can be observed in the variance of cluster 2 between the two experiments (Fig. 3). These differences in variance were statistically significant (P , 0.001). The large interquartile ranges for cluster 2 in experiment 1 contrast with much smaller interquartile ranges in experiment 2 (   RNA1 3  3  3  3  AbV16 RNA2 3  3  3  3  AbV16 RNA3 3  3  3  3  AbV16 RNA4 3  3  3  3  ORFan8  3  3  3  3  AbV14  4  4  4  4  AbV9  4  4  4  4 a Hierarchical clusters were determined by visual inspection of Fig. 1a. evidence of bimodality by comparing the likelihood of fit to either a normal distribution or double normal distribution (indicating bimodality) using expectation maximization ( Fig. 4). AbV16 RNAs 1 to 4 from cluster 3 and AbV6 RNA2 from cluster 2 fit a double normal distribution in both experiments and, therefore, were present in two distinct populations. The difference between the two peaks of abundance in AbV16 RNAs of cluster 3 was greater than 1,000-fold (.10 on the 40 2 DC T scale, i.e., 2 10 ) between the median values of high and low abundance. In addition, in experiment 1, AbVSV, AbV10, AbV12, and AbV6 RNA1 from cluster 2 had a better fit for a double normal distribution, whereas in experiment 2, a different set of additional viral RNAs, ORFans 2, 3, 7, and 8, fit better the double normal distribution (Fig. S6).
Relationship between viral abundance and symptom expression. Pearson's pairwise correlations were calculated to identify possible correlations between symptom expression (mushroom cap color [log 10 DE]) and viral RNA abundance (40 2 DC T ) both collectively as clusters F5-6 ( Fig. 5) and as individual RNAs (Fig. 6). Significant correlations were established between the abundance of cluster 3 RNAs and cap color in both experiments 1 and 2, with the higher viral RNA abundance correlating with greater cap browning (Fig. 5). Cluster 1 abundance also significantly correlated with color but in experiment 2 only. Examination of correlations between individual RNA abundance and cap color enabled greater precision, and significant correlations were identified for the AbV16 components for both experiments and for ORFan8 and MBV in experiment 2 ( Fig. 6 and Table S3).
Comparison of AbV16 viral levels in fruitbodies and underlying mycelium. RNA levels of AbV16 RNA 1 and RNA 2 were compared between 16 fruitbodies (8 symptomatic brown and 8 white mushrooms) and the underlying mycelium directly below these 16 mushrooms, a subset of the total data. By blocking the variance due to the difference between the RNA types, a clear difference was identified in AbV16 abundance between brown and white fruitbodies (t = 5.149, df = 43, P , 0.001), but no difference was found in AbV16 abundance between the mycelia below brown and white fruitbodies (t = 0.488, df = 43, P = 0.628) (Fig. 7).

DISCUSSION
Agaricus bisporus mycelium (strain A15) harboring a relatively low number of viral RNAs was infected with additional viral RNAs from a small inoculum of MVX-4569 A. bisporus mycelium. The viral RNAs moved and replicated throughout the A15 mycelium, including the mycelium of the casing layer from where fruitbodies developed following environmental manipulation. Eighteen viral RNAs representing nine viruses were detected in the cultures, and their abundances were analyzed to elucidate their relationships and how they interact. The nine viruses in the cultures represent nine distinct RNA-dependent RNA polymerase (RdRp) enzymes with the potential for competing for resources required for viral replication. However, these RNAs could be grouped into four clusters using unsupervised statistical measures of their abundances. The grouping of AbV2 and MBV into cluster 1 was less strong. The identification of clusters suggests that virus levels had reached equilibrium in the interval between infection and fruitbody sample harvest (20 days).
The abundance of the four components of AbV16 correlated with the disease symptom, cap browning. The levels of AbV16 RNAs also interacted positively with viral levels of a different cluster, cluster 2, by pairwise correlation. In addition, the components of AbV16 and AbV6 RNA2 (cluster 2) exhibited two distinct populations of RNA levels, low and very high, orders of magnitude apart. We propose that the disease phenotype is caused by interaction between these clusters leading to the step-jump in viral RNA levels. However, we acknowledge the possibility of the inverse, that the changes in viral abundance are a response to the cellular pathology. There appears to be a fitness cost to the step-jump in RNA levels, as AbV16 abundance can crash in culture after exhibiting the disease culture (23). High levels of AbV16 have been found early during fruitbody development, and thereafter there is no evidence of transition between the two viral replicative states; levels remain either low or very high (23). Here, we show that before fruitbody formation, AbV16 levels were low in the mycelium. At the point of environmental stimulation for differentiation to form fruitbodies, 15% of the developing fruitbodies underwent an abrupt shift to 1,000 times higher levels of AbV16, and for the remaining fruitbodies the levels stayed low.
The clusters represent evidence of synergism and viral coreplication, probably resulting from different mechanisms. Cluster 3 RNAs have previously been proposed as a single virus, AbV16 and the tightly associated ORFan8 (10), and therefore is replicated by the single RdRp of the cluster. Cluster 1 (ORFans 2, 3, 5, and 7 and, more loosely, AbV2 and MBV), cluster 2 (AbSV, AbV10, AbV12, AbV6 RNA1, and AbV6 RNA2), and cluster 4 (AbV9 and AbV14) consist of phylogenetically diverse viruses with multiple RdRps, which, therefore, are unlikely to be replicated by a single replicase. The mechanisms by which these RNAs cluster, implying coordinated replication, are not fully understood  Table S2 in the supplemental material). but may involve complementary suppression or saturation of host antiviral defenses or forms of molecular interactions that are similar to the collective infectious units that enhance comovement and replication (18,21,28,29). In C. parasitica, the suppression of the host antiviral silencing mechanism by the hypovirus CHV4 enabled the stable  The replication of all positive-strand RNA viruses is generally associated with cellular membranes that are heavily modified or reorganized to act as efficient platforms for genome amplification (30). Major membrane modification is consistent with the cause of the browning reaction (by pathology or mechanical damage), i.e., the oxidation of phenols by tyrosinase, which are separated into different compartments by membranes (31). Thus, membrane modification and reorganization plausibly lead to leakiness, in turn promoting the reaction and the formation of the brown-colored products. The exponentially higher replicative capacity of the step-jump in viral levels may involve additional cellular factors, such as the centrosome-microtubule-binding domain gene, which has been reported to be upregulated in A. bisporus with high viral levels, or structural proteins encoded by the viruses, possibly acting as scaffolds for viral replication complexes (23). The two viral RNAs with known structural protein domains are AbV6 RNA2 and MBV, each of which encodes capsid-like domains (10,32). RNA silencing is the primary antiviral defense mechanism in fungi; however, it is unknown whether its suppression alone could account for the 1,000-fold increase in viral levels or whether an additional mechanism of high-throughput viral replication is involved. The step-jump in viral RNA replication appears to be generic and not sequence specific, as a number of additional RNAs underwent this change in one experiment only.
Correlations both positive (clusters 2 and 3) and negative (clusters 1 and 4) were identified between individual clusters (Fig. 2). A lack of correlations also was observed between cluster 1 and clusters 2 and 3 and between cluster 4 and clusters 2 and 3, indicating neutral or no interactions. The mechanism for the antagonism between clusters 1 and 4 could be negative dominance, whereby defective interfering RNAs (such as the ORFans of cluster 1) interfere with viral replication in a cell by faster but defective replication (21,33). RNA silencing can contribute to the production of defective interfering RNA (20). Cluster 1 consists of AbV2, MBV, and ORFans 2, 3, 5, and 7. AbV2 has been described as hypovirus-like based on its 28% sequence identity at the protein level to Cryphonectria hypovirus 2 (10). There is evidence that ORFan3 is subgenomic to AbV2, as they are joined in some sequence assemblies (10). Cluster 1 included many RNAs previously described as ubiquitous and asymptomatic, consistent with the observation that persistent viruses minimize the negative effects of parasitizing their host or confer an advantage on the host in compensation, the trade-off hypothesis (34,35). An alternative mechanism for the antagonism could be mutual exclusion, i.e., that viruses of each cluster are spatially separated in different cells (18). This hypothesis was not tested in these experiments, as the samples consisted of whole fruitbodies.
MBV did not consistently group with the clusters. The hierarchical and k-medoid clustering methods placed MBV with cluster 1 (Fig. 1), but it falls outside clusters from principal component analysis (see Fig. S3 in the supplemental material), and pairwise correlation showed MBV correlating with RNAs of cluster 3 (Fig. 2). Unlike AbV16, MBV did not have a bimodal distribution of abundance, although its abundance correlated with cap browning in experiment 2. MBV, a particulate virus, also has been found as a coinfection of AbV1, where it appears to exploit incidents of high viral replication (36).
There appears to be a stochastic nature to the viral interactions and the development of pathogenesis. For instance, (i) brown symptomatic fruitbodies were distributed apparently randomly in beds of growing white mushrooms, (ii) large differences were identified in variance of cluster 2 RNAs between experiments 1 and 2 (Fig. 3A), and (iii) the viruses with double normal distributions differed between experiments 1 and 2 (Fig. 4). That stochastic factors could play a role in the complexity of the multifactorial and dynamic state of a large number of viral RNAs is perhaps not surprising, an otherwise jostling equilibrium tipped over the edge by the microenvironment. The RNAs that consistently demonstrated the exponential step-jump comprise two multipartite viruses, AbV16 and the capsid-encoding AbV6 RNA2, and the stochastic event may result from a combination of these RNAs at the same cellular location at a specific phase in the cell (duplication) cycle during early fruitbody formation. In mycelial cells, AbV16 and AbV6 RNA2 have been shown to be separately located using fluorescence in situ hybridization (37).
Rather than competing for the cells' resources, the numerous viruses interacted synergistically as groups, as shown by their similar relative abundances. Further interactions between groups were identified as synergistic, antagonistic, and neutral. This statistical and modeling study has revealed simultaneous viral interactions in the same fungal culture. The study did not examine mechanisms of interaction, but the data are consistent with previously described mechanisms, such as suppression of antiviral defenses, defective interfering RNA, spatial separation, and replication by a single RdRp. It is unknown whether these mechanisms can collectively account for the 1,000fold increase in AbV16 and AbV6 RNA2 levels. Host gene expression changes, viral translocation, and effects of environmental stresses will be subject to further research. Multiple viral infections are likely to be prevalent in nature, although they may be undetected, as they are asymptomatic. However, the random expression of severe symptoms may be a factor in the uncertainties associated with biological control agents, such as fungi and viruses, in controlling fungal diseases.

MATERIALS AND METHODS
Biological material. Two sources of Agaricus bisporus were used, (i) commercial spawn of Sylvan strain A15 (Sylvan Inc., Kittanning, PA, USA) and (ii) the virus-infected strain 4569 (MVX-4569), which was isolated and cultured from a diseased fruitbody originally collected in 2005 from a commercial mushroom farm growing strain A15 (23) and maintained in the laboratory as mycelial subcultures per Grogan et al. (24). Experimental inoculum of MVX-4569 was produced as per Fleming-Archibald et al. (38). The inoculum was stored at 4°C for 5 months between the two experiments.
Experimental conditions. Two experiments were carried out under the same growth conditions but at different scales: experiment 1, large scale, where mushrooms were grown on a metal shelf system (0.24 m by 0.24 m by 4 m), reflecting common industry practice for growing mushrooms, and experiment 2, small scale, where mushrooms were grown in three polypropylene crates (0.6 m by 0.4 m by 0.24 m).
In both experiments, pasteurized mushroom compost, known commercially as phase II compost, was inoculated with 0.5% (wt/wt) commercial A15 spawn. For experiment 1, the metal shelf was filled with 160 kg of inoculated phase II mushroom compost. For experiment 2, three polypropylene crates were each filled with 15 kg of inoculated phase II mushroom compost. For each experiment after the filling process, the mycelium was allowed to colonize the compost for 17 days in a controlled environment at 25°C, .6,000 ppm CO 2 , and 90 to 95% relative humidity per standard mushroom industry practice.
The fully colonized compost was then inoculated with MVX-4569-infected compost fragments by incorporating them into the surface of the A15 colonized compost. In experiment 1, 2 g of MVX-4569infected compost fragments was incorporated into the first 0.5-m length of compost in the trough (approximately 20 kg). In experiment 2, 1.5 g was incorporated uniformly into the surface layer of the colonized compost in the polypropylene crates (15 kg). The effective rate of inoculation in both cases was 0.01%, as described by Fleming-Archibald et al. (38).
The compost for both experiments was then covered with 55 mm of casing soil (Harte Peat; Clones Co., Monaghan, Ireland), incubated for a further 7 days under the same conditions as detailed above. Mushroom fruitbody production was initiated by an environmental manipulation, reduction of CO 2 levels to 1,200 ppm and temperature to 18°C (38). Mushroom fruitbodies from the first flush were harvested for experimental use 20 days after casing application. In both experiments, the inoculation with MVX-4569 resulted in the production of symptomatic brown-colored mushrooms (10 to 15%) randomly distributed among white asymptomatic mushrooms largely in the zone where inoculation had taken place.
Fruitbody and compost sampling. In experiment 1, 40 mushroom fruitbodies were harvested from the growing shelf; 10 of these were visibly the brownest in color, and three white fruitbodies near each brown mushroom were also selected and harvested. For experiment 2, 49 fruitbody samples were harvested with various degrees of color symptoms, and 16 compost samples were taken directly below specific fruitbodies, 8 from below fruitbodies displaying the brown symptom and the remaining 8 from below asymptomatic fruitbodies, white in color. To sample the compost, the fruitbody was harvested, the casing layer was carefully removed, and approximately 5 g of the compost below was collected. All samples were harvested by hand. After color measurement, the samples were freeze-dried for approximately 24 h until completely dry and stored at 220°C.
Mushroom cap color determination. Mushroom cap color was measured using a CR-410 Chromameter (Konica Minolta, Basildon, UK) as three color parameters, L*, b*, and a*, which were then combined as a single parameter, DE (39).
RNA extraction from fruitbodies. RNA was extracted from mushroom fruitbodies using a modified version of our previously published methodology (10). Briefly, 80 mg dry weight of mushroom fruitbody was mixed with 4Â (wt/vol) lysis buffer (STE, 1% SDS) and 4Â (wt/vol) 5:1 phenol-chloroform (pH 4.5) (Fisher Scientific, Loughborough, UK), and cells were lysed using bead beating (Retsch, Hann, Germany). RNA was precipitated by adding 0.1 volumes of 3 M sodium acetate, pH 5.2, and 2 volumes of absolute ethanol. The RNA was then cleaned by running through 10% polyvinylpolypyrrolidone (PVPP) spin columns. If the flowthrough retained a strong color, it was applied to a further PVPP column and the procedure repeated. RNA quality and quantity were measured using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Loughborough, UK) and the samples stored at 280°C.
RNA extraction from compost. RNA was extracted from compost using a methodology similar to that for fruitbodies. A volume of 75 mg of freeze-dried compost was mixed with 1 ml TRIzol (Sigma-Aldrich, Haverhill, UK) and cells lysed by bead beating. RNA was then recovered by binding to 50 ml of silica gel in suspension. The silica was pelleted by pulse spinning on a microcentrifuge for 10 s and RNA cleaned by sequentially washing in (i) 1 ml wash buffer (100 ml 0.1 M Tris-HCl, pH 6.4, and 120 g guanidine thiocyanate), (ii) 1 ml wash buffer, (iii) 1 ml 70% ethanol, (iv) 1 ml 70% ethanol, and (v) 1 ml acetone with pulse spinning and discarding of the supernatant at each step. RNA was resuspended in 50 ml RNase-free water and stored at 280°C.
DNase treatment and reverse transcription. The resulting nucleic acid was treated with RNase-free DNase 1 (Fisher Scientific, Loughborough, UK) by adding 1 ml DNase 1 (1 U/ml), 1 ml 10Â reaction buffer, nucleic acid at 1,000 ng, and nuclease-free water up to 10 ml in a nuclease-free 0.2-ml PCR tube (Fisher Scientific, Loughborough, UK) on ice. These were incubated at 37°C for 60 min. To inactivate the DNase, 1 ml 50 mM EDTA was added and incubated at 65°C for 10 min. The resulting RNA was reverse transcribed to cDNA by adding to a nuclease-free 0.2-ml PCR tube on ice, 1.2 ml random hexamer primer (Thermo Scientific, Loughborough, UK), 2 ml deoxynucleotide triphosphates (5 mM each) (Fisher Bioreagents), 1 ml RNA, and 8.8 ml nuclease-free water. This was incubated at 65°C for 5 min. It was then placed on ice for at least 1 min. Four microliters of 5Â first-strand buffer, 1 ml 0.1 M dithiothreitol, 1 ml nuclease-free water, and 1 ml reverse transcriptase (Superscript III; 200 U/ml) (Invitrogen, Loughborough, UK) was added, and the tubes were incubated at 25°C for 5 min, 50°C for 60 min, and 70°C for 15 min.
qPCR. The qPCR assays were performed in duplicate using reaction mixtures containing 7.5 ml SYBR green master mix (Applied Biosystems, Loughborough, UK), 1 ml diluted cDNA, 15 pM forward and reverse primers, and nuclease-free water to 15 ml. In addition, nontemplate water control was included for each primer set on every plate. Detection and quantification were achieved with an Applied Biosystems 7500 real-time PCR system (Applied Biosystems, Loughborough, UK). The following amplification cycles were used: 2 min at 50°C, 10 min at 95°C, and 40 cycles of 15 s at 95°C and 60 s at 60°C. The cycle threshold (C T ) values for the viral RNAs were normalized against the 18S gene to determine the DC T values used for analysis.
Statistical analysis. Several statistical approaches were taken to identify novel associations between viruses and, by inference, their interactions. All statistical analyses were carried out in R 3.4.0 (40). Two unsupervised clustering methods, hierarchical and k-medoids (41), and principal component analysis (PCA) were applied to the abundance data of all 89 fruitbody samples to identify viral clusters. Properties of viral clusters were investigated using further techniques. Pearson's pairwise correlations were used to identify (i) positive and negative interactions between viruses and their clusters, (ii) the relationship between viral abundance and fruitbody color, and (iii) the relationship between RNA abundance in fruitbodies and the underlying mycelium. Differences in the variance of abundance between (i) clusters, (ii) experiments, and (iii) their interaction were tested by fitting a linear model to the absolute deviance from the median value of each virus per experiment. Significance of each term in the model was calculated with analysis of variance (ANOVA). Estimated marginal means and post hoc comparison test were calculated with emmeans (42) with comparison P values adjusted for multiple testing using the Tukey method (43). To determine whether the sample abundances of each RNA fit a single or double normal distribution, single and double normal models were fit to the data using expectation maximization implemented in the mixtools package (44); significance was calculated using bootstrapping by producing 10,000 iterations of the likelihood ratio statistic for testing the null hypothesis-single normal fit versus the alternative hypothesis-double normal fit with significance for rejecting the null taken at the 0.05 level.
To examine any differences in the abundances of AbV16 RNA1 and RNA2 between underlying mycelia and fruitbodies, diseased and asymptomatic, data were fit to a linear model, and analysis of variance was used to identify differences between (i) the RNAs, (ii) fruitbody color, (iii) tissue type (fruitbody or mycelium), and (iv) the interactions (RNA-color-tissue) while controlling for paired sample (fruitbody and underlying mycelium). The significance of each term in the model was calculated with ANOVA. Estimated marginal means and post hoc comparison test were calculated with emmeans (42) with contrast P values controlled for familywise error rate using the Tukey method (43).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.