De novo transcriptomic assembly and mRNA expression patterns of Botryosphaeria dothidea infection with mycoviruses chrysovirus 1 (BdCV1) and partitivirus 1 (BdPV1)

Background Pear ring rot, caused by Botryosphaeria species, is responsible for substantial economic losses by causing severe recession of pear tree growth in China. Mycovirus-mediated hypovirulence in plant pathogenic fungi is a crucial biological method to control fungal diseases. Methods We conducted a large-scale and comprehensive transcriptome analysis to identify mRNA in B. dothidea in response to mycovirus. De novo sequencing technology from four constructed libraries of LW-C (Botryosphaeria dothidea chrysovirus 1, BdCV1), LW-P (Botryosphaeria dothidea partitivirus 1, BdPV1), LW-CP (LW-1 strain infection with BdCV1 and BdPV1), and Mock (free virus) was used to investigate and compare gene expression changes in B.dothidea strains infected with mycovirus. Results In total, 30,058 Unigenes with an average length of 2128 bp were obtained from 4 libraries of B. dothidea strains. These were annotated to specify their classified function. We demonstrate that mRNAs of B. dothidea strains in response to mycovirus are differentially expressed. In total, 5598 genes were up-regulated and 3298 were down-regulated in the LW-CP group, 4468 were up-regulated and 4291 down-regulated in the LW-C group, and 2590 were up-regulated and 2325 down-regulated in the LW-P group. RT-qPCR was used to validate the expression of 9 selected genes. The B. dothidea transcriptome was more affected by BdCV1 infection than BdPV1. We conducted GO enrichment analysis to characterize gene functions regulated by B. dothidea with mycovirus infection. These involved metabolic process, cellular process, catalytic activity, transporter activity, signaling, and other biological pathways. KEGG function analysis demonstrated that the enriched differentially expressed genes are involved in metabolism, transcription, signal transduction, and ABC transport. mRNA is therefore involved in the interaction between fungi and mycovirus. In addition, changes in differential accumulation levels of cp and RdRp of BdCV1 and BdPV1 in B. dothidea strains were evaluated, revealing that the accumulation of BdCV1 and BdPV1 is related to the phenotype and virulence of B. dothidea strain LW-1. Conclusions The identification and analysis of mRNAs from B. dothidea was first reported at the transcriptome level. Our analysis provides further insight into the interaction of B. dothidea strains infection with chrysovirus 1 (BdCV1) and partitivirus 1 (BdPV1) at the transcriptome level. Electronic supplementary material The online version of this article (10.1186/s12985-018-1033-4) contains supplementary material, which is available to authorized users.

Pear ring rot, caused by the destructive pathogen Botryosphaeria, is responsible for substantial economic losses through widespread fruit rots and stem canker. This caused a severe recession in the growth of pear fruit trees in China [1][2][3][4][5]. Due to the lack of disease-resistant varieties, many fungicides are recommended for treatment of Botryosphaeria canker, which induces drug resistance and environmental pollution [3]. It is urgent to find new, safe, and effective means to control pear ring rot disease. Mycovirus-mediated hypovirulence in plant pathogenic fungi is a powerful method to control fungal crop diseases such as the hypovirulent strain of C. parasitica (CHV1) to heal cankers of chestnut blight, and Sclerotinia sclerotiorum hypovirulence-associated DNA virus 1 (SsHADV-1) to control the rapeseed stem rot [6][7][8][9][10][11][12][13]. Recent studies of fruit tree fungal diseases have reported that the mycovirus Rosellinia necatrix megabirnavirus 1 (RnMRV1) has potential to control white root rot (Rosellinia necatrix) diseases of fruit trees [8]. In recent years, several dsRNA mycoviruses have been identified and sequenced in B. dothidea strains isolated from pear trees with ring rot and stem wart symptoms in China [14,15]. Chrysovirus BdCV1 and partitivirus BdPV1 mixed infection results in hypovirulent B. dothidea LW-1 strain in our study, revealing BdCV1 to be the first mycovirus found to be responsible for hypovirulence (reduced virulence), growth rate, and phenotypic sectorization of the phytopathogenic fungus B. dothidea [16]. BdCV1 is therefore a good candidate for biological control of apple and pear ring spot diseases, while BdPV1 does not cause any visible symptoms in virulent pathogenicity [5,16].
High throughput sequencing technology provided comprehensive information on gene expression. Measuring the expression patterns of mRNAs at the transcriptome level from fungi infection with mycovirus is vital to reveal the mechanism of mycovirus-mediated hypovirulence in pathogenic fungi [6,17,18]. Only a limited number of studies, using high throughput sequencing technologies and bioinformatics, have demonstrated transcriptional or translational changes in fungi infection with mycovirus, such as Cryphonectria parasitica, Aspergillus fumigatus, Neurospora crassa, Sclerotinia sclerotiorum, and Fusarium graminearum [17][18][19][20][21][22][23][24]. For instance, 150 genes representing a wide spectrum of biological functions were down-regulated in the strain Ep-1PN by Sclerotinia Sclerotiorum Debilitation-associated RNA Virus (SsDRV), of which S. sclerotiorum integrin like gene (SSITL) was verified affect pathogenic fungus virulence and growth rate [20,23]. Genes associated with viral replication, maintenance of viral life cycle, transcription and translation machinery, and signal transduction were up-regulated while those including membrane-associated transporters and cellular transport systems were down-regulated based on a genome-wide transcriptome analysis of F. graminearum infected with Fusarium graminearum virus 1 (FgV1) [18,21,22,25].
The main objective of our study is to characterize the mRNAs from B. dothidea involved in biological processes associated with the host infecting with mycovirus. Currently, there are no studies published on the B. dothidea transcriptome. No direct and detailed functional genomics resources in public databases for understanding the molecular mechanism of B. dothidea exist except the draft genome sequences from two B. dothidea strains reported, the pathogen of Apple ring rot disease and isolated from grapevine host, respectively [26,27]. Therefore, large-scale transcriptome sequencing using Illumina sequencing technology was performed to explore the potential mRNAs expression in B. dothidea related to pathogenic factors and hypovirulence determinants in response to two mycoviruses BdCV1 and BdPV1.

Botryosphaeria dothidea strains and infection with BdCV1 and BdPV1
Three B. dothidea strains, LW-1, LW-C, and LW-P, infected with hypovirulent and non-hypovirulent mycovirus were generated [5,16]. LW-1 (or designated as LW-CP) was infected with Botryosphaeria dothidea chrysovirus 1 (BdCV1) and Botryosphaeria dothidea partitivirus 1 (BdPV1), which was isolated and identified from the trunks of sand pear'Jinshuiyihao' cv., collected in Wuhan city, Hubei province of China [5]. The LW-P strain infection with BdPV1 and LW-C with BdCV1 were derived from LW-1 by single hyphae isolation, respectively, as previously described [5,16]. Meanwhile, mock derived from LW-C strain by protoplast isolation was virus-free as negative control. The fungal strains were grown at 25°C in darkness on potato dextrose agar (PDA) medium and stored on PDA plate at − 4°C or in glycerin at − 80°C.

Total RNA extraction from B. dothidea strains
The mycelia of the mycovirus-infected three strains and the virus-free strain of B. dothidea were harvested on 5 day of PDA culture. The mycelia were frozen and ground using liquid nitrogen for nucleic acid extraction. Total RNAs were isolated and obtained with TRIzol Reagent (Invitrogen, USA), phenol and chloroform to remove proteins and ethanol precipitation according to manufacturer's instructions. Total RNAs digested by DNase I (Takara, Dalian, China) were precipitated with ethanol and dissolved in DEPC treated water to be used as template for Next Genomics Sequencing (NGS). The RNA integrity Number (RIN > 7) and concentration of obtained tRNAs were further assessed using the Bioanalyzer Agilent 2100 and NanoDrop spectrophotometer (Agilent, USA). High quality samples were used to construct the four de novo libraries.

cDNA library construction and de novo sequencing
For each RNA-seq library, total RNA from three biological repetitions were pooled. Following the total RNA extraction and treatment with DNase I, mRNA was isolated and enriched by Oligo (dT). mRNA was cut into small fragments and cDNA synthesis was performed with random hexamer primers. PCR amplification was then prepared for library construction according to the Illumina RNA library protocol for transcriptome sequencing. Then, the quantification of the four cDNA libraries was verified by the Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System. Finally, the libraries were sequenced using the HiSeq 4000 instrument (Illumina, USA) at the Beijing Genomics Institute (BGI) (Shenzhen, Guangdong province, China).
Transcriptome data processing, de novo assembly, and functional annotation Internal software was used to filter raw reads to generate clean reads data via the following processes: 1) The removal of adaptor-polluted reads, 2) removal of reads containing a high content (> 5%) of unknown bases (N), 3) removal of low-quality reads. Following filtering, the remaining reads were used for downstream analyses. Trinity was used to perform de novo assembly with clean reads which PCR duplication removed, and Tgicl was used to cluster transcripts to Unigenes. The sequence by Trinity was called transcripts. Then, gene family clustering with Tgicl for each sample's Unigene to obtain final Unigenes was performed for downstream analyses. The Unigenes were divided into two classes, including clusters and singletons. The prefix used for clusters was CL with the cluster id behind it (in one cluster, there are several Unigenes with more than 70% similarity between them). The prefix used for singletons was Unigene. For unigene functional annotation analysis, NT (ftp:// ftp.ncbi.nlm.nih.gov/blast/db), NR (ftp://ftp.ncbi.nlm.nih.gov/blast/db), GO (http://geneontology.org), COG (http:// www.ncbi.nlm.nih.gov/COG), KEGG (http://www.genome.jp/kegg), SwissProt (http://ftp.ebi.ac.uk/pub/databases/ swissprot), and InterPro (http://www.ebi.ac.uk/interpro) were used as functional databases. Blasts aligned unigenes to NT (NCBI non-redundant nucleotide sequences), NR (NCBI non-redundant protein sequences), COG (Cluster of Orthologous Groups), KEGG (Kyoto Encyclopedia of Genes and Genome), and SwissProt to obtain the annotation. Blast2GO with NR annotation to obtain the GO (Gene Ontology) annotation, and InterProScan5 was used to get the InterPro annotation [28][29][30]. The functional databases with a priority order of NR, SwissProt, KEGG, and COG were selected to determine the sequence annotation of unigenes. Coding sequences (CDS) were extracted from sequences of unigenes that best mapped to functional databases. ESTScan was used to predict coding regions of unigenes if no hits were obtained in any database mentioned above [31].

Differential gene expression analysis of mRNA in response to mycovirus
De novo sequencing from four constructed cDNA libraries designated as LW-CP, LW-C, LW-P, and Mock was performed. To analyze the differentially expressed genes (DEGs) between two samples in response to mycovirus, clean reads to assembled unigenes using Bowtie2 were mapped [32]. Then, gene expression value was calculated with RSEM [33]. DEGs were identified with PossionDis based on the Poisson distribution, which was performed as described by Audic et al. [34]. The parameters of PossionDis was set as |Fold Change | > 2.0-fold (|log 2 FC| > 1) with a false discovery rate (FDR) < 0.001.

Clustering analysis of DEGs of B. dothidea strains infection with mycovirus
With DEGs, pheatmaps of hierarchical clustering were performed using the heatmap.2 package in R software [35]. For clustering more than two groups, we used the intersecting DEGs.

GO and KEGG pathway analysis of DEGs
Gene Ontology (GO) classification and functional enrichment was performed on DEGs. GO displayed three ontologies including molecular function, cellular component, and biological process. In addition, KEGG pathway functional enrichment was also analyzed using a corrected P-value less than 0.05 for significantly enriched DEGs.

RT-qPCR analysis of gene mRNA expression level
RT-qPCR was used to validate DEG expression levels obtained from sequencing and to assess the expression levels of cp and RdRp from BdCV1 and BdPV1. Primers used are listed in Additional file 8. Total RNA was extracted from the four groups: Mock, LW-CP, LW-C, and LW-P cultured for 5 d, 10 d and 15 d using TRIzol (Invitrogen, USA) according to manufacturer's instructions. Total RNA was digested by DNase I and used as template (Takara, Dalian, China). The cDNA was obtained to perform reverse transcription using PrimeScript™ RT reagent Kit with genomic DNA Eraser (Takara, Dalian, China). According to the manufacturer's instruction, qPCR was performed using 2.5 μl diluted cDNA, 10 μl of the SYBR Premix Ex Taq II PCR mixture (Tli RNaseH Plus) (Takara, Dalian, China), 1 μl of each 5 mM primer, and deionized water to a final volume 20 μl. All reactions were run in triplicate. The qPCR was performed on the Bio-Rad iQTM5 Real-time System machine (BIO-RAD, USA). The products were verified by melt curve analysis. The 18 s gene from B.dothidea strains was used as the internal reference for standardization of cDNA expression levels from samples. mRNA relative expression level changes were quantified by a comparative C T method (ΔΔC T ) using the formula 2 -ΔΔCt [36].

Differential accumulation of cp and RdRp in B. dothidea infection with BdCV1 and BdPV1
To determine the culture time point to analyze the B. dothidea transcriptome in response to mycovirus infection, the relative expression levels of cp and RdRp mRNA from BdCV1 and BdPV1 in LW-1, LW-C, and LW-P cultured for 5 d, 10 d, and 15 d were quantified by RT-qPCR analysis. The results demonstrate that the relative expression levels of cp and RdRp mRNA from BdCV1 and BdPV1 are expressed differentially in B. dothidea strains LW-C, LW-P, and LW-1 infected with single and mixed mycovirus. In general, BdCV1 viral accumulation increased, while BdPV1 viral accumulation decreased slowly in LW-1 strain cultured periods from 5 to 15 d (Fig. 1).

De novo sequencing and unigene assembly of B. dothidea strains
The mycelia of the mycovirus-infected three strains and the virus-free strain of B. dothidea at 5 d grown on PDA were collected for RNA extraction to be further used for De novo sequencing. The B. dothidea transcriptome was sequenced and assembled de novo. This is a valuable resource since there are no current genomic data available for B. dothidea except two B.dothidea strains reported [26,27]. Four libraries designated as LW-CP, LW-C, LW-P, and Mock were constructed to conduct the comprehensive genome-wide transcriptome analysis following Illumina Hiseq 4000 platform sequencing. In total, clean reads [47. were obtained after trimming, and then assembled using de novo assembly program of Trinity (Additional file 1: Table S1). The assembly of the clean reads resulted in 26,461 transcripts for LW-C, 30,707 for LW-CP, 25,808 for LW-P, and 25,298 for Mock. Tgicl was used to cluster transcripts to unigenes. The clustering quality metrics are shown, including of 23,056 unigenes for LW-C, 27,554 for LW-CP, 22,403 for LW-P, and 22,068 for Mock, which reveals the effect of mycovirus infection on a number of gene expression at the transcriptome level. In total, 30,058 unigenes with an average length of 2128 bp and N50 length of 3338 bp were obtained from 4 libraries of B. dothidea strains (Additional file 2: Table S2). The size distribution of these unigenes ranged from 300 bp with 3889 unigenes to higher than 3000 bp with 7593 unigenes as shown in Fig. 2. To validate the unigene sequence by de novo assembly, RT-PCR and sequencing were used to measure seven unigenes selected randomly. Results demonstrated that they were more than 98% identity comparable to de novo sequencing (Additional file 3: Table S3). This is the first report of B. dothidea transcriptome.
Functional annotation of predicted proteins from B. dothidea unigenes Thirty thousand fifty eight non-redundant sequences were submitted to a BLASTx search and annotated to classify function from seven public functional databases including NR, NT, GO, COG, KEGG, Swissprot, and Interpro. In total, the number of unigenes annotated to least one functional database was 24 Table 2).
Searching the KEGG database revealed that 19,722 unigenes (65.61%) matched to KEGG pathways. The assembled unigenes were assigned to five specific functional categories, including cellular processes, environmental information processing, genetic information processing, metabolism, and organism systems. The most represented KEGG function category was the "global an overview maps" (7138 unigenes) and the 'carbohydrate metabolism' (3927 unigenes), belonging to the metabolism cluster. 'Signal transduction' (776 unigenes) and 'membrane transport' (122 unigenes) were annotated. 'Drug resistance' (5 unigenes) was the smallest group (Fig. 3c). These results reveal that B. dothidea infection with mycovirus mainly involve the metabolism pathway.

Differential expression profiles of mRNAs from B. dothidea strains infection with mycovirus
Transcriptional changes in fungal genes following mycovirus infection have been little reported [17,18,22]. Here, we examined genome-wide transcriptional dothidea genes that show significant differential gene expression with functional annotation are summarized in Table 3. In addition, we detected the expression of genes involved in gene silencing pathways. These were up-regulated in response to mycovirus infection, such as argonaut-like genes (Ago), dicer-like genes (Dicer), and RNA-directed RNA polymerase (RdRp) genes (Additional file 5: Table S5). Genes annotated heat shock protein (Hsp)-related proteins involved in 'response to stimulus' pathway were also up-regulated (Additional file 6: Table S6). The above results reveal that the numbers of up-regulated transcripts were more prevalent than down-regulated   transcripts. Notably, more genes were differentially expressed under BdCV1 compared to BdPV1. A hierarchical clustering of the mycovirus-responsive transcriptome demonstrates the significant similarities and differences in gene expression profiles among the three groups (LW-CP, LW-C, and LW-P; Fig. 4b). We next compared DEGs to determine the numbers of mycovirus-specific or commonly expressed genes in the three groups (Fig. 4c). A total of 2483 genes were identified in the three groups, revealing a potentially critical biological process. The number of specifically expressed genes in each mycovirus-infected sample was 2203 (1876 up-regulated and 602 down-regulated) for LW-CP, 2057 (881 up-regulated and 1501 down-regulated) for LW-C, and 1177 (1059 up-regulated and 746 down-regulated) for LW-P. Although many DEGs were mycovirus-specific, 2483 DEGs were common to LW-CP, LW-C, and LW-P, suggesting that they could be involved in stress response. The B. dothidea transcriptome analysis revealed that many host genes responsive to viral infection were common and mycovirus-specific expression [17,38].
Increasing the stringency of differentially expressed genes progressively up to FDR < 0.001 and log 2 FC ≥ 2 (FC ≥ 4) also reveals a high number of significantly expressed transcripts (Additional file 7: Figure S7). DEGs were the highest in LW-CP (5313 genes) and the lowest in LW-P (2,738 genes). Furthermore, the hierarchical clustering and Venn

Validation of mRNA transcriptome analysis data by RT-qPCR
To confirm RNA-Seq results, we selected nine genes and purified total RNAs from different biological samples. The results of RNA-Seq were highly consistent with those of RT-qPCR using specific primers (Additional file 8: Table S8). For example, the expression of CL1217.Contig1 and CL1346.Contig4_All were strongly down-regulated, and CL2349Contig4_All, CL51.Contig10_All, CL1218.Conti-g1_All, Unigene3107_All, CL4042.Contig3_All, Unige-ne4082_All, and CL5019.Contig2_All were up-regulated in response to mycoviruses, which were all confirmed by RT-qPCR. This reveals that the transcript patterns by RT-qPCR are consistent with the RNA-Seq data analysis; however, the extent of differential expression differed ( Table 4).

Effect of mycovirus infection on expression of predicted transcription factors
Because many fungal genes are unknown, only a limited number of transcription factors (TFs) were annotated [17,18,22]. In this study, 1083 TFs belonging to 18 families, including Zn2Cys6 (684 TFs) and C2H2 zinc finger (165 TFs), were predicated from de novo transcriptome data, many of which were differentially expressed in response to chrysovirus 1 BdCV1 and partitivirus 1BdPV1 infection. In total, 11 putative TF families were differentially expressed in response to different mycovirus infections. One hundred twenty eight TFs were commonly identified in all three mycovirus-infected groups by Venn diagrams. Up-regulated TFs were more numerous than down-regulated TFs in the LW-CP and LW-P groups, while more numerous TFs were down-regulated in the LW-C group.   Table S9 and Additional file 10: Figure S10).

GO enrichment analysis of mycovirus-responsive unigenes from B. dothidea strains
To determine whether the mycovirus-responsive genes are involved in specific pathways, a functional classification of the DEGs was performed by GO enrichment analysis. By comparing LW-CP versus Mock, 4866 DEGs were functionally assigned to the relevant terms in three categories including biological processes, cellular component, and molecular function. The GO terms 'metabolic process' (GO: 0008152, 690 genes) in biological process, 'catalytic activity' (GO: 0003824, 817 genes) in molecular function, and 'membrane' (GO: 0016020, 235 genes) in the cellular component category were highly enriched (Fig. 5a).

KEGG pathway analysis of mycovirus-responsive genes from B. dothidea strains
To elucidate the biological processes, DEGs from B.dothidea strains regulated by mycovirus infection, and KEGG pathway were analyzed. The top 20 enriched pathways are Note: "/" indicated no significantly differential expression summarized from the 3 groups. In the LW-CP group, 2554 DEGs were significantly enriched (corrected p-value < 0.05) in 'metabolic pathways' (PATH: ko01100), and 830 DEGs were related to 'biosynthesis of secondary metabolites' (Fig. 6a). However, this is approximately 36 and 12% of the total genes involved in 'metabolic pathways' and 'biosynthesis of secondary metabolites' , respectively (Additional file 12: Table S12A). In the LW-C group, similar to the results of LW-CP, genes associated with 'metabolic pathways' (PATH: ko01100) and 'biosynthesis of secondary metabolites' (PATH: ko01110) (Fig. 6b), approximately 39 and 13% of the DEGs, respectively, were significantly enriched (Additional file 12: Table S12B). 'Stach and sucrose metabolism' (PATH: ko00500) related to carbohydrate metabolism, 'Glycerophospholipid metabolism' (PATH: ko00564) related to lipid metabolism, and 'Indole diterpene alkaloid biosynthesis' (PATH: ko00403) related to biosynthesis of other secondary metabolites, were all significantly enriched. Lastly, MAPK signaling pathway of signal transduction (PATH: ko04011) and ABC transporters of membrane transport (PATH: ko02010) were enriched. In the LW-P group, more DEGs were enriched in the 'metabolic pathways' and 'Starch and sucrose metabolism' (PATH: ko00500), with approximately 39 and 11% of the DEGs significantly enriched among the DEG, respectively (Fig. 6c). In addition, 'peroxisome' (PATH: ko04146), with functional  Table S12C). From the above functional analysis, the DEGs associated with metabolic pathways, transport and catabolism, membrane transport, transcription, and signal transduction are presumed to play a role in the interaction of fungal defense and mycovirus counter defense strategies.

Discussion
In this report, the transcriptome of B. dothidea strains was sequenced using the Illumina platform. A total of high-quality 18.87 Gb bases were obtained. We then assembled 4 fungi samples occurred in B. dothidea infection with mycovirus and obtained 30,058 Unigenes (≥200 bp).
The N50 and GC content of unigenes was 3338 bp and 56.32%, respectively. We then annotated the unigenes using 7 functional databases. 24,836 (82.63%) unigenes were annotated with at least one functional database (Tables 1 and 2 Previous research has demonstrated that mycovirus infections change the expression of a broad range of genes and cause hypovirulence or phenotypic alterations in the fungal host except plant [17-19, 22, 38, 39]. C. parasitica infected with Cryphonectria Hypovirus 1 (CHV1) is used as a model for studies on virus/host interaction. cDNA library analysis of differentially expressed genes involved viral replication, transmission, host growth, development, and defense. The pro1 and CpBir1 genes have important biological functions for hypovirus and chestnut blight fungus host [17,40,41]. In our study, RNA-Seq-based genome-wide expression profiling analysis in response to single or mixed mycovirus chrysovirus 1(BdCV1) and partitivirus 1 (BdPV1) infection was performed. The expression patterns of 9 putative genes involved in mycovirus stress measured by quantitative real-time PCR were consistent with their transcript changes as identified by RNA-seq (Table 4). It revealed specific and common alterations in host gene transcript accumulation following infection of B. dothidea by BdCV1 and BdPV1 [17]. As expected, more transcriptional changes occurred in response to the hypovirulent LW-CP and LW-C compared to the LW-P (Table 3; Fig. 4), revealing that chrysovirus 1(BdCV1) would have a stronger effect than partitivirus 1 (BdPV1) on the B. dothidea transcriptome, and gene expression changes in transcriptome caused by hypovirulent and non-hypovirulent mycoviruses were related to the observed host phenotypes and pathogenicity [17,18]. Meanwhile, changes in differential accumulation of BdCV1 and BdPV1 in B. dothidea strains were demonstrated that the accumulation level of BdPV1 decreased by co-infection with BdCV1 (Fig. 1). It revealed that the BdCV1 and hypovirulent determinants may inhibit the expression level of BdPV1, prediction that it is the main reason to induce to the phenotype and hypovirulence of B. dothidea LW-1 strain (Figs. 1 and 4), which need to be further verified.
It is known that transcription factors play important roles in fungi response to mycovirus. Mycovirus can activate TFs expression differentially [18]. In this study, 1083 transcription factors classified into 18 families were predicted. The above data revealed many TFs from B. dothidea were expressed differentially among the mycoviruses, especially Zn2Cys6 and C2H2 zinc finger families (Table 5). This corroborates the results from the F. graminearum transcriptome, whose differentially expressed TFs included fungal-specific and dominant Zn2Cys6 and the C2H2 zinc finger members involved in transcriptional regulation [18,22,42,43]. The results demonstrate that TFs are down-regulated in response to BdCV1 in LW-C, revealing that BdCV1 inhibits TFs expression in B. dothdiea (Table 5; Additional file 9: Table S9). In addition, a class of heat shock proteins (HSPs) was up-regulated in response to mycovirus infection (Additional file 6: Table S6), which supports the speculation that HSPs homologues are influenced by mycovirus [17].
To better survey the biological behavior of defense response, it is necessary to understand the functional distribution of these DEGs in B. dothidea following mycovirus infection based on the transcriptome and bioinformatics analysis [44][45][46][47]. Through the enriched GO terms analysis, it is also revealed that in biological processes, DEGs mainly distributed to metabolism. Differential expression of genes related to metabolism might be associated with the host phenotype [17,39,48]. In addition, GO terms for genes associated with signal transduction, including phosphatidylinositol signaling system and MAPK signaling pathway were more enriched in response to infection by BdCV1 compared to BdPV1. This reveals that it is necessary to determine the role of the MAPK signaling pathway in the regulation of mycovirus and host interactions [25]. Genes involved in membrane, oxidoreductase activity, RNA biosynthetic processing, and ribosomal assembly were enriched in B. dothidea following mycovirus infection ( Fig. 5; Additional file 11: Table S11), which indicates the diversity of genes affected by mycoviral infection. Furthermore, KEGG pathway analysis uncovered DEGs involved in important pathways. A mass of metabolism pathways, both primary and biosynthesis of secondary metabolism, are identified and significantly enriched in response to mycovirus infection ( Fig. 6; Additional file 12: Table S12). Similar metabolic pattern is exhibited in different fungi-virus combinations, which indicated these pathways play important role in fungi host in response to mycoviral infection [17,18,22]. In addition, 'carbohydrate metabolism' , 'lipid metabolism' , 'membrane transport system' , 'transport and catabolism' , 'translation' and 'signal transduction' which were highly enriched ( Fig.6; Additional file 12: Table S12). It also provides insight into the various biological pathways associated with mycoviral infection with plant pathogenic fungi.
RNAi machinery is involved in the regulation of fungi and mycovirus infection by controlling endogenous and exogenous RNA [49]. Indeed, the biological functions of the RNA silencing pathway have been characterized in the Neurospora crassa, Cryphonectria parasitaica, Rosellinia necatrix, and F. graminearum [50][51][52][53][54][55][56]. As reported, FgAGO-1 and FgDICER-2 are responsible for hairpin RNA-triggered RNA silencing and related small interfering RNA accumulation in F. graminearum [18,49,51]. In this study, RNA interference components, including dicer-like (Dicer), argonaute-like (Ago), and RdRp genes in B. dothidea were expressed differentially in response to mycovirus infection (Additional file 5: Table S5). Small RNA sequencing demonstrates that microRNA exists and is expressed differentially in B. dothidea infection with mycovirus (not published). This suggests that B. dothidea possesses RNA silencing pathways for antiviral defense and endogenous gene regulation. These data indicate that mycovirus may activate host antiviral defense. Furthermore, it is necessary to determine and analyze the RNA silencing component responsible for transcriptional regulation and antiviral defense mechanism in B. dothidea.

Conclusions
In this study, 30,058 unigenes were obtained from hypovirulent and non-hypovirulent B. dothidea strains infected with mycovirus by de novo assembly. To identify potential