- Open Access
Comprehensive analysis of full genome sequence and Bd-milRNA/target mRNAs to discover the mechanism of hypovirulence in Botryosphaeria dothidea strains on pear infection with BdCV1 and BdPV1
IMA Fungus volume 10, Article number: 3 (2019)
Pear ring rot disease, mainly caused by Botryosphaeria dothidea, is widespread in most pear and apple-growing regions. Mycoviruses are used for biocontrol, especially in fruit tree disease. BdCV1 (Botryosphaeria dothidea chrysovirus 1) and BdPV1 (Botryosphaeria dothidea partitivirus 1) influence the biological characteristics of B. dothidea strains. BdCV1 is a potential candidate for the control of fungal disease. Therefore, it is vital to explore interactions between B. dothidea and mycovirus to clarify the pathogenic mechanisms of B. dothidea and hypovirulence of B. dothidea in pear. A high-quality full-length genome sequence of the B. dothidea LW-Hubei isolate was obtained using Single Molecule Real-Time sequencing. It has high repeat sequence with 9.3% and DNA methylation existence in the genome. The 46.34 Mb genomes contained 14,091 predicted genes, which of 13,135 were annotated. B. dothidea was predicted to express 3833 secreted proteins. In bioinformatics analysis, 351 CAZy members, 552 transporters, 128 kinases, and 1096 proteins associated with plant-host interaction (PHI) were identified. RNA-silencing components including two endoribonuclease Dicer, four argonaute (Ago) and three RNA-dependent RNA polymerase (RdRp) molecules were identified and expressed in response to mycovirus infection. Horizontal transfer of the LW-C and LW-P strains indicated that BdCV1 induced host gene silencing in LW-C to suppress BdPV1 transmission. To investigate the role of RNA-silencing in B. dothidea defense, we constructed four small RNA libraries and sequenced B. dothidea micro-like RNAs (Bd-milRNAs) produced in response to BdCV1 and BdPV1 infection. Among these, 167 conserved and 68 candidate novel Bd-milRNAs were identified, of which 161 conserved and 20 novel Bd-milRNA were differentially expressed. WEGO analysis revealed involvement of the differentially expressed Bd-milRNA-targeted genes in metabolic process, catalytic activity, cell process and response to stress or stimulus. BdCV1 had a greater effect on the phenotype, virulence, conidiomata, vertical and horizontal transmission ability, and mycelia cellular structure biological characteristics of B. dothidea strains than BdPV1 and virus-free strains. The results obtained in this study indicate that mycovirus regulates biological processes in B. dothidea through the combined interaction of antiviral defense mediated by RNA-silencing and milRNA-mediated regulation of target gene mRNA expression.
Pear (Pyrus) is an important fruit crop, with rich pear germplasm resources including thousands of cultivars belonging to five domesticated species of Pyrus pyrifoila, Pyrus communis, Pyrus bretschneideri, Pyrus ussuriensis and Kuerle pear, which are grown all over China (Wu et al. 2018). However, pear stem wart and stem canker diseases induced by Botryosphaeria species are widespread across all pear producing areas, inducing rotten fruit, branch warts, branch dieback, stem canker, and pear tree growth recession. B. paeva, B. rhodina, B. obstuse and B. dothidea are the dominant species that seriously affect pear crop yield and quality (Slippers et al. 2004; Slippers and Wingfield 2007; Garibaldi et al. 2012; Yan et al. 2013; Zhai et al. 2014; Wang et al. 2017b; Yang et al. 2017b). Botryosphaeria spp., belonging to Botryosphaeriaceae in the class of Dothideomycetes, destroy many types of fruit trees such as apple, pear, grapevine, peach, plum, apricot, begonia, chestnut and jujube (Slippers and Wingfield 2007; Garibaldi et al. 2012; Tang et al. 2012; Yan et al. 2013; Zhai et al. 2014). B. dothidea is the main destructive pathogen and dominant strain responsible for pear ring rot disease in pear-growing areas (Zhai et al. 2014). Despite the importance of B. dothidea, little is known about the genetics and pathogenic mechanism of this pathogen (Wang et al. 2018a, 2018b). To date, the main strategy used to reveal the mechanism of phytopathogenic infection has involved sequencing of the complete genome and annotation of gene function in fungal diseases such as rice sheath blight, rice blast, and wheat stripe rust pathogen fungus, Cenococcum geophilum and Cercospora sojina (Dean et al. 2005; Zheng et al. 2013a, 2013b; Dong et al. 2015; Peter et al. 2016; Luo et al. 2018; Plissonneau et al. 2018). Similarly, sequencing of the Ustilago esculenta genome and gene annotation revealed the molecular mechanism of the artificial selection and smut fungi-host interaction involved in RNA-silencing in host defense (Ye et al. 2017). Therefore, full genome sequencing and analysis of the B. dothidea LW-Hubei strains infecting pear hosts will provide important molecular information that can be used to explore the pathogenic mechanisms.
Although the application of numerous chemical fungicides and fruit bagging techniques have contributed to the prevention pear ring rot disease, the associated environmental pollution has serious detrimental effects on food safety and health (Guo et al. 2009; Zhao and Cao 2012; Zhou et al. 2012). Thus, new, safe and effective measures for the biocontrol of pear ring rot disease induced by B. dothidea are urgently needed. Mycovirus-mediated hypovirulence in plant pathogenic fungi has become an important strategy for the biocontrol of fungal diseases (Anagnostakis 1982; Nuss 2005; Chiba et al. 2009; Ghabrial and Suzuki 2009; Pearson et al. 2009; Yu et al. 2013; Xie and Jiang 2014). Recently, hypovirulence of B. dothidea LW-C strain infected with B. dothidea chrysovirus 1 (BdCV1) has been reported as a potential candidate for the biocontrol of fungal diseases (Wang et al. 2014; Wang et al. 2017b) Studies of the molecular basis of hypovirulence of BdCV1-infected B. dothidea have provided further clarification of the pathogenesis of B. dothidea and provided a theoretical basis for the prevention and control of disease caused by this fungus.
RNA-silencing is a conserved mechanism in eukaryotic organisms, in which small RNAs (sRNAs), particularly microRNAs (miRNAs), regulate transcriptional or post-transcriptional gene silencing (Kim et al. 2009; Huntzinger and Izaurralde 2011; Bologna and Voinnet 2014). sRNAs are involved in regulation of morphological and developmental processes in response to abiotic and biotic stresses (Dang et al. 2011; Wang et al. 2016a, 2016b). However, studies on sRNA-mediated gene silencing in fungi and mycoviruses are rare (Li et al. 2010; Nuss 2011; Chiba et al. 2013). RNA-silencing as a mechanism of sRNA-based antiviral defense requires the involvement of Dicer, Ago, RdRp genes, which regulate sRNA biosynthesis and mycovirus expression in fungi (Kadotani et al. 2004; Segers et al. 2007; Janbon et al. 2010; Campo et al. 2016). In fungi, sRNAs suppress plant immunity by hijacking host RNAi pathways, as reported for Botrytis cinerea, Rhizoctonia solani, Verticillium wilt and Puccinia striiformis (Wang et al. 2016a, 2016b; Weiberg et al. 2013, 2014; Wang et al. 2017a, 2017c; Cai et al. 2018). sRNAs can be transported between plants and fungi to silence virulence gene, indicating that resistance based on cross-kingdom sRNA transport has extensive application potential for the prevention and control of various plant diseases (Weiberg et al. 2013, 2014; Cai et al. 2018). High-throughput sequencing technologies and bioinformatics analysis have been used to identify microRNAs (miRNAs) in some plant pathogenic fungi, such as Rhizoctonia solani, Neurospora crassa, Puccinia striiformis, Cryptococcus neoformans, Trichoderma reesei, Sclerotinia sclerotiorum, and Fusarium graminearum, F. oxysporum, and Penicillium marneffei (Lee et al. 2010; Jiang et al. 2012; Zhou et al. 2012; Kang et al. 2013; Lau et al. 2013; Chen et al. 2014, 2015; Dahlmann and Kück 2015; Mueth et al. 2015; Lin et al. 2016). Mature miRNAs negatively regulate gene expression by complementary binding to the open reading frame (ORF) or untranslated region (UTR) of specific target genes, although the mechanism of miRNA-mediated regulation of target mRNAs remains poorly understood in most fungal species (Lee et al. 2010; Li et al. 2010; Dang et al. 2011; Huntzinger and Izaurralde 2011). RNAi components and sRNAs have been identified in P. striiformis. Pst-milR1 was identified as an important pathogenicity factor that impairs wheat resistance by regulation of the pathogenesis-related 2 gene (PR2) in wheat, and PsDCL was found to regulate Pst-milR1 production (Mueth et al. 2015; Wang et al. 2017a). Therefore, miRNA-mediated regulation of biological processes may be involved in different biological pathways in host fungi (Li et al. 2010; Huntzinger and Izaurralde 2011).
As the causative agent of pear ring rot, B. dothidea is associated with serious agricultural losses. Some studies of the genome from B. dothidea infecting with apple and grapevine hosts, and transcriptome of B. dothidea infecting pear hosts in response to mycovirus, have been reported; however, the sRNAs in B. dothidea remain to be identified (Garibaldi et al. 2012; Zhai et al. 2014; Liu et al. 2016; Wang et al. 2018a, 2018b). The existence of microRNA-like RNAs (milRNAs) in B. dothidea and the potential mechanism by which these molecules regulate their target genes to counter pathogen infections remains to be established. In our previous study, we used de novo transcriptome sequencing to show that the DCL, Ago and RdRp genes were expressed in B. dothidea strains in response to mycovirus infection (Wang et al. 2018b). RNA-silencing has been predicted to exist in B. dothidea strains, although this remains to be confirmed. Furthermore, the roles of RNA-silencing components in sRNA transcription and milRNA generation in B. dothidea infected with mycovirus remain to be clarified.
Mycoviruses are used for biocontrol, especially in fruit tree disease. BdCV1 (Botryosphaeria dothidea chrysovirus 1) and BdPV1 (Botryosphaeria dothidea partitivirus 1) belonging to dsRNA mycovirus, influence the biological characteristics of B. dothidea strains. BdCV1 is a potential candidate for the control of fungal disease. Therefore, it is vital to explore interactions between B. dothidea and mycovirus to clarify the pathogenic mechanisms of B. dothidea and hypovirulence of B. dothidea in pear. In this study, we sequenced the entire genome of B. dothidea and constructed four sRNA libraries of B. dothidea infected with mycovirus to examine the evidence for RNAi induction and milRNA expression in response to mycovirus in hypovirulent and non-hypovirulent B. dothidea strains. This information is important in providing insights into milRNAs as potential pathogenic factors and elucidating the molecular basis of the interaction of B. dothidea with mycovirus.
MATERIALS AND METHODS
B. dothidea strains, conidia and hyphal morphological analysis by light microscopy
The Botryosphaeria dothidea LW-1 (designated as LW-CP) strain, with a mixed infection of Botryosphaeria dothidea chrysovirus 1 (BdCV1) and Botryosphaeria dothidea partitivirus 1 (BdPV1), was isolated from sand pear (Pyrus pyrofolia) trunk in Hubei Province, China, as previously described (Wang et al. 2014). The virulent Botryosphaeria dothidea LW-P strain infected with BdPV1, the attenuated Botryosphaeria dothidea LW-C strain infected with BdCV1, and Botryosphaeria dothidea Mock strain (designated as LW-Hubei) as a virus-free control, were used in this study (Wang et al. 2017b). In addition, we used the HL-1 and HBWH-1 B. dothidea strains isolated from sand pear trunks collected from Hubei Province (Wang et al. 2014; Wang et al. 2017b).
Colonies of the four B. dothidea strains (LW-CP, LW-C, LW-P and Mock) were cultured for 3 days on potato dextrose agar (PDA) and Murashige Skoog (MS) at 25 °C in darkness, respectively. To induce sporulation, the mycelia were harvested and exposed to blank light for 11 d. Sporulation and conidial angle were observed under a stereo microscope (Olympus, SZX16; Japan). Conidia were stained with DAPI (Solarbio, Beijing of China) and conidial size was measured from images obtained under DAPI filter using a fluorescence microscope (Leica DM2500, Germany) (Tebeest et al. 1989). In addition, the mycelia were cultured for 3 d on agarose gel to observe the morphology and assess the effect of mycovirus on the biological characteristics of B. dothidea.
Evaluation of the cellular morphology of B. dothidea strain mycelia by transmission electron microscopy (TEM)
The four samples of mycelium from B. dothidea strains were cultured for 48 h before collection. The samples were then fixed with 2.5% glutaraldehyde for 2 h at room temperature (RT), rinsed three times with 0.1 M phosphate-buffered saline (pH 7.4) and subsequently fixed in 1% osmic acid for 2 h at 4 °C. The four samples were dehydrated in a graded ethanol series (50, 70, 80, 85, 90, 95 and 100%), permeated successively with 2:1 and 1:1 acetone and epoxy as a penetrating agent, and embedded in ethoxyline resin. Ultra-thin sections (thickness, 60–100 nm) were prepared with a diamond microtome (Leica, EM UC7, Germany) and collected onto Formvar-coated copper grids. The sections were then stained with lead and uranium for observation of the cellular morphology and damage to the phytopathogenic fungi caused by mycovirus infection using a TEM (America FBI company, Cat. Tecnai G2 20 TWIN) under a 200-kV accelerating voltage.
Evaluation of B. dothidea strain pathogenicity
The pathogenicity of B. dothidea strains was evaluated on apple and pear fruit, and pear branches as previously described with a little modification (Wang et al. 2014). The mycelial plugs from strains were cultured on PDA for 3 d to obtain fresh samples for inoculation. The inoculation was performed using the wounding method and the expansion of the lesion area was observed and measured.
Vertical transmission of BdCV1 and/or BdPV1 to B. dothidea LW-Hubei strain
Single conidia were obtained from the induced pycnidia of B. dothidea strains LW-C, LW-P and LW-CP generated as described in the previous section. The derivatives were screened for BdCV1 and/or BdPV1 dsRNAs to assess vertical transmission ability and stability of BdCV1 and BdPV1 in B. dothidea strains (Yang et al. 2017a, 2017b).
Horizontal transmission between BdCV1 and BdPV1 in B. dothidea strains
Transmission of hypovirulence traits of strain LW-CP, LW-C and LW-P was assessed according to a previous method with some modifications (Wang et al. 2017b). B. dothidea strains were co-cultured at 25 °C for 7 days to allow contact between the two colonies. The hypovirulent strains LW-CP and LW-C, and the virulent strain LW-P served as the donors, and the virus-free strain Mock served as the recipient. In addition the LW-C and LW-P strains were also used as both recipient and donor, respectively, to evaluate BdCV1 and BdPV1 transmission between the two strains. After incubation of the contact cultures, mycelia agar plugs from the colony margin of the Mock, LW-C and LW-P strains were transferred to a fresh PDA plate to obtain three derived isolates from each co-cultured recipient strain. The derivatives were screened for dsRNAs as described previously (Yang et al. 2017a). The parental strains LW-CP, LW-C and LW-P were included as controls.
B. dothidea strain information and genomic DNA isolation
The mycelia were cultured for 5 d at 25 °C on PDA containing 50 μg/ml ampicillin and streptomycin. The mycelia were then collected and ground in liquid nitrogen. Genomic DNA from the Mock strain (designated as the LW-Hubei isolate) were isolated using the standard cetyltrimethylammonium bromide (CTAB) method with some modifications (Murray and Thompson 1980). The gDNA titer and quality were assessed by Qubit Fluorometer, agrose gel electrophoresis and NanoDrop spectrophotometer prior to library construction.
Genome sequencing, assembly and annotation
The genomic DNA of the B. dothidea LW-Hubei isolate was sequenced using the Single Molecule Real-Time (SMRT) method at the Beijing Genomics Institute (BGI) in Shenzhen, China (Chin et al. 2013). A 270-bp paired-end library and a 10-kb mate-pair DNA library were constructed, and subsequently subjected to paired-end 150-bp sequencing using the Illumina HiSeq4000 platform. The treated sequencing data (filtered reads: 6.84 G, sequencing depth: 147×) were used to estimate the genome size, repeat content, and genome heterozygosis by K-mer analysis and then confirmed by alignment between assembled scaffolds and original reads. The 20-kb library was then constructed according to standard PacBio methods, including genomic DNA fragmentation, end-repair, adaptor ligation, and template purification. The 20-kb library was quantified using a 2100 Bioanalyzer (Agilent, USA) and sequenced by SMRT. The raw data generated by the PacBio platform were treated according to standard protocols. The reads were assembled using variety of software in combination with bioinformatics analysis of read assemblies, contig linkage to scaffolds and gap-filling (Kim et al. 2014; Badouin et al. 2015; Faino et al. 2015; Sit et al. 2015; Tsuji et al. 2015). The sequencing data (filtered reads: 3.72G, sequencing depth: 80×) were assembled by CANU (Version-1.2) with default parameters and Falcon (Version-0.3.0).
Protein-encoding genes were annotated by a combination of independent ab initio predictors based on homology (Johnson et al. 2008; Stanke et al. 2008; Ter-Hovhannisyan et al. 2008). The transcriptome data obtained from the four libraries of LW-Hubei strains were mapped to the genome using Trinity to assess the quality of B. dothidea genome assembly and annotation (Wang et al. 2018b).
Functional annotation of predicted genes
Functional annotation of all predicted genes was performed using multiple databases including Swiss-Prot, NR (NCBI non-redundant protein sequences), KEGG (Kyoto Encyclopedia of Genes and Genome), GO (Gene Ontology), COG (Cluster of Orthologous Groups) and KOG by BlastP with a cut-off E-value of 10− 5 (Ashburner et al. 2000; Jones et al. 2014; Galperin et al. 2015; Makarova et al. 2015; Kanehisa et al. 2016). In addition, potential virulence-related proteins were identified by searching against the pathogen-host interaction database (PHI) (Torto-Alalibo et al. 2009). The following databases were used to analyze genes from the B. dothidea LW-Hubei isolate: fungal cytochrome P450 database, carbohydrate-active enzymes (CAZy) database (http://www.cazy.org/), virulence factor database (VFDB), type III secretion system effector protein (T3SS), and transporters (Fischer et al. 2007; Vargas et al. 2012; Levasseur et al. 2013; Chen et al. 2016; Elbourne et al. 2017). Putative kinase proteins were predicted and annotated using the databases (http://ekpd.biocuckoo.org/download/HMM/kinase.hmm).
Repeat annotation, DNA methylation analysis and synteny annotation
A de novo repeat database of the B. dothidea LW-Hubei isolate was generated using Repbase, ProMask, the De novo database and the Tandem Repeats Finder (TRF) (Benson 1999). Transposable elements (TEs) were identified by aligning the assembly with the known transposon sequence database and the De novo database. The software package used for predication repeat annotation contains Repeat Masker (Version: 4–0-6), Repeat Protein Masker and the De novo database (Smit et al. 2014). Tandem repeats (TR) were predicted by TRF (Version: 4.04). Whole-genome DNA modification detection and motif analysis were performed using the PacBio SMRT software (version 2.3.0, http://www.pacb.com) (Fang et al. 2012; Davis et al. 2013; Manso et al. 2014; Adhikari and Curtis 2016; Cohen et al. 2016). The whole genome syntenic analysis of B. dothidea LW-Hubei isolate compared with Macrophomina phaseolina, Neofusicoccum parvum, Diplodia corticola and Diplodia seriata was performed with MUMmer (version 3.0) (Fu et al. 2012). Linear synteny graphs were constructed at the nucleic acid and amino acid levels (E-value ≤1e− 5, sequence identity ≥85% and the best hit for each protein were selected). In addition, a schematic representation of the synteny between the selected scaffold1 of the B. dothidea LW-Hubei isolate with those of the four strains was also generated with CIRCOS (http://mkweb.bcgsc.ca/circos/).
Comparative genomic analysis of gene families
The protein sequences of the B. dothidea LW-Hubei isolate genome was compared with those of four pathogenic species of Botryosphaeriaceae family (M. phaseolina, N. parvum, D. corticola and D. seriata). OrthoMCL (version 2.0) was used to analyze the orthologs, co-orthologs and in-paralog pairs. Gene families were constructed and analyzed based on the genes from B. dothidea LW-Hubei isolate and the other fungi as reference strains. In addition, single-copy gene family analysis was performed according to the BGI protocols, including alignment of the protein sequence in BLAST, gene family TreeFam clustering with Hcluster_sg software, and conversion of the alignment proteins into those of the multiple amino acid sequences in the CDS region using MUSCLE-3.8 (Edgar 2004a, 2004b).
Phylogenetic tree analysis
The phylogenetic relationships between B. dothidea LW-Hubei and the four Botryosphaeriaceae isolates (M. phaseolina, N. parvum, D. corticola, D. seriata) as well as four additional Dothideomycete isolates (Cenococcum geophilum, Dothistroma septosporum, Exserohilum turcica, Lepidopterella palustris) were analyzed at amino acid level based on the consistent phylogenetic core-pan gene and single-copy orthologous families. The sequence alignments were performed using MUSCLE-3.8(Edgar 2004b; Nandi et al. 2010).
The phylogenetic tree was generated by MEGA 7 using a TreeBeST method, with the bootstrap value set at 1000 replicates (Kumar et al. 2016).
RNA preparation and isolation
The mycelia of the mycovirus-infected three strains and the virus-free strain of B. dothidea were collected after culture on PDA for 5 days and then ground under liquid nitrogen. The total RNAs (tRNAs) were extracted using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions with some modifications. The concentration and integrity of the tRNAs was further assessed using a NanoDrop spectrophotometer and agarose gel electrophoresis. High quality tRNAs were used for high-throughput sequencing of sRNA, verification of milRNAs and mRNA expression analysis.
Small RNA library construction and deep sequencing
Small RNA libraries were constructed using 20 μg of tRNA and the kits of PAGE gel separation of RNA segments of different size to cut out between 18 and 30 nt strip (Illumina, USA) according to the manufacturer’s protocols for sRNA fraction isolation, adapter ligation, reverse transcription (RT) (Invitrogen) in combination with adapter-specific RT-primers, PCR amplification and purification by PAGE gel electrophoresis. The quality and yield of purified high-quality cDNA library were detected using Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System, and sequenced on an Illumina Genome AnalyserIIx (BGI, Shenzhen, China).
Identification of milRNAs and expression by bioinformatics analysis
High-quality reads obtained from the four B. dothidea strains by elimination of the low-quality reads, adaptors and other contaminants to obtain 18–30 nt clean reads. sRNAs with perfect genomic matches were used for further analysis. All trimmed sequences between 18 and 30 bp in length were used to align clean reads (tags) to the reference genome in the study and to other sRNA in databases to identify known milRNA by Bowtie (Langmead et al. 2009). Alignments with no gap and two mismatches between the query sequences and known milRNAs were used in subsequent analysis. The reads that did not match any known miRNAs were further processed to discover putative novel milRNAs with their predicted hairpin precursors using miRDeep2 software (Friedlander et al. 2008). To further confirm precursor stem-loop structure of these predicted novel milRNAs, the putative precursor sequence was subjected to folding analysis using RNAfold software. In addition, the sRNA expression levels were estimated in terms of Transcripts Per Kilobase Million (TPM). The milRNA gene expression levels were used directly to determine differential gene expression compared with that of the B. dothidea strains in response to mycovirus infection. To evaluate the significance of differential gene expression, FDR ≤0.001 and the absolute value of |log2FC| ≥ 1 were set as the default threshold.
Bd-milRNA and mRNA expression analysis in B. dothidea strains
The Bd-milRNAs were detected by stem-loop RT-PCR (Chen et al. 2005). Briefly, a stem-loop RT primer was used to reverse-transcribe mature milRNAs to obtain cDNAs. First, 500 ng tRNAs and 0.25 μM of each individual stem-loop RT primer were added to be 10 μl RT reaction mixture and incubated at 95 °C for 5 min before placing on ice immediately. Following the addition of 2 U RNase inhibitor (Yeasen, Shanghai city of China), 0.5 mM dNTPs (Takara, Dalian of China), cDNA was generated using the PrimeScript™ RT Kit with genomic DNA Eraser (Takara, Dalian of China) in a 20 μl RT reaction mixture, which was incubated at 42 °C for 60–90 min.
MilRNAs were amplified by polyadenylation-mediated PCR. sRNA was isolated using the miRcute miRNA Isolation Kit (Tiangen) according to the manufacturer’s protocol. The template was polyadenylated, and reverse transcribed to cDNA using the Poly (A) miRNA cDNA Synthesis Kit (Invitrogen) according to the manufacturer’s instructions. MilRNAs were amplified from the cDNA pool by using the mature miRNA-specific forward and reverse primers provided in the kit. The PCR products were then separated by agarose gel electrophoresis.
In addition, B. dothidea mRNA expression was evaluated by RT-qPCR analysis using 2.5 μl of diluted cDNA as a template and the primers described in Additional file 17: Table S1 and Additional file 18: Table S2. The reaction mixture contained 10 μl of the SYBR Premix Ex Taq II PCR mixture (Tli RNaseH Plus) (Takara, Dalian of China), 1 μl of each 5 mM forward and reverse primer, and deionized water to a final 20 μl volume. All reactions were run in triplicate on a CFX96 Real-time System (BIO-RAD, USA). The products were verified by melting curve analysis. The actin gene (GME8592_g) was included as an internal reference for each sample. The relative mRNA expression level changes were calculated using a comparative CT method (ΔΔCT) according to the formula 2-ΔΔCT (Livak and Schmittgen 2001).
Prediction of milRNA targets and function annotation
The potential target mRNAs of known and novel milRNAs were predicted using two types of software. The predicted milRNA target cleavage sites were identified by alignment of the mature milRNA sequences with reported B. dothidea RNA-seq data, which was reanalyzed based on the complete genome sequence of B. dothidea LW-Hubei using psRNATarget with default parameter values set as a cut-off threshold of 2.5 with the expectation of higher prediction coverage. In addition, TargetFinder software was used to predict the putative target mRNAs (Audic and Claverie 1997; Fahlgren and Carrington 2010; Dai and Zhao 2011). The KEGG pathway was also used to perform pathway enrichment analysis of the predicted biological function of the differentially expressed target genes (Kanehisa et al. 2008). In addition, the functional enrichment analysis was performed using WEGO online (http://wego.genomics.org.cn/) (Ye et al. 2006).
All experiments were performed in triplicates or more. Data were presented as mean ± SD and analyzed using SPSS (Version 21.0) software. Data were further evaluated using Duncan’s new multiple range test. P < 0.05 was considered to indicate statistical significance.
The interaction between dsRNA mycoviruses and B. dothidea strains
BdCV1 and BdPV1 influenced the biological characteristics of the B. dothidea LW-Hubei strain. In the current study, BdCV1 reduced the virulence of LW-C and LW-1 (designated as LW-CP) inoculation on ‘Malus domestica cv. Fuji’ apple fruit, Pyrus pyrifolia cv. ‘Hohsui’ pear fruit and one-year branch growth, whereas BdPV1 enhanced the virulence of LW-P. Virus-free Mock (designated as LW-Hubei strain), HL-1 and HBWH-1 exhibited virulence as controls (Additional file 1: Figure S1). These findings are in accordance with our previous results showing hypovirulence of LW-C and LW-1 in ‘Huangguang’, ‘Kuerle’, and Pyrus bretschneideri Rehd. cv. ‘Mili’ pear fruit and P. pyrifolia Nakai cv. ‘Hongxiangsu’ pear branches (Wang et al. 2014; Wang et al. 2017b). These results further demonstrated that BdCV1 is associated with hypovirulence of B. dothidea infecting different pear cultivars and apple hosts; thus, BdCV1 is implicated as a promising candidate for biological control of pear ring rot disease.
In addition, we analyzed the effects of BdCV1 and BdPV1 on conidial bodies, microscopic morphology and cell ultrastructure of the mycelium from B. dothidea LW-Hubei strains. The LW-P strain produced conidiomata and conidium angle more effectively than the LW-C and LW-CP strains when cultured on PDA and MS medium, respectively (Fig. 1a, Additional file 2: Figure S2). The LW-P and Mock strains produced solitary, globular pycnidia covered with mycelia, whereas the condial bodies of the LW-CP and LW-C strains formed a yellow mass. The conidia produced by LW-P were longer than those produced by LW-C, LW-CP and Mock (Fig. 1b, Additional file 19: Table S3). Observation of the conidial growth and morphology revealed that the four strains produced hyaline conidia, with septa rarely formed before germination, and predominantly containing eight nuclei (Fig. 1b, Additional file 3: Figure S3). Observation of the mycelia produced after 48 h in culture by transmission electron microscopy (TEM) showed that the cell structure of the LW-C and LW-CP strains was markedly different from that of the Mock and LW-P strains. The cells of the LW-C and LW-CP strains were small and deformed, with irregular thickening and some degradation of the cell wall, exosome formation and plasmolysis of the protoplasm, swollen mitochondria, and heavily vacuolization. In contrast, the cell wall and membrane in the LW-P and Mock strains were uniform and smooth, and the nucleus profile was clearly visible (Fig. 2). Microscopic observation of B. dothidea strains cultured on PDA for 3 d showed that the hyphae from LW-C and LW-CP were deformed and with more branches compared with those of LW-P and Mock (Additional file 4: Figure S4).
BdCV1 and BdPV1 transmission to B. dothidea hosts was assessed by the detection of mycoviral dsRNA (Yang et al. 2017a). The results showed that BdCV1 and/or BdPV1 exist stably in LW-1 (109 of 109), LW-C (84 of 84) and LW-P (94 of 94) strains following vertical transmission of the offspring conidia from B. dothidea strains (Additional file 5: Figure S5). In addition, the horizontal transmission in co-cultures of LW-C/BdCV1 + LW-P/BdPV1 was also analyzed by detection of mycoviral dsRNA in the LW-C and LW-P strains. The results demonstrated that BdCV1 was transmitted to LW-P, whereas BdPV1 was not transmitted to LW-C (Fig. 3). In contrast, in the LW-CP/Mock, LW-C/Mock, and LW-P/Mock control groups, BdCV1 and/or BdPV1 were effectively transmitted to the Mock strain at higher frequencies than that reported previously (Wang et al. 2014, Wang et al. 2017b).
These results revealed that BdCV1 has a marked effect on the growth, development and pathogenic properties of B. dothidea during pear infection; however, the mechanisms underlying these effects remain to be clarified. To further understand the interactions between B. dothidea and BdCV1 and BdPV1, it is important to obtain a complete knowledge of the genome of B. dothidea LW-Hubei strain infecting the pear host.
Sequencing and assembly of the B. dothidea LW-Hubei strain genome
To analyze the pathogenic mechanism by which B. dothidea infects pear, total genomic DNA from the mycelia of the B. dothidea LW-Hubei isolate was extracted and sequenced by the Single Molecule Real-Time (SMRT) technique in combination with the HiSeq platform. Three different sizes of insert libraries (270 bp, 10 kb and 20 kb) were sequenced to generate 6835 Mb trimmed data with a genome coverage of 147× for the Illumina and 3721 Mb (genome coverage of 80×) for the PacBio sequencing platform (Additional file 20: Table S4 and Additional file 21: Table S5). For the PacBio platform, subread distribution analyses confirmed the high quality of the 20-kb library (Additional file 6: Figure S6). Evaluation of the heterozygosity and distribution by k-mer revealed a volume peak at 30, indicating high homozygosity of the B. dothidea LW-Hubei isolate genome. In total, 68 contigs were generated with an N50 length of 2.72 Mb and 49 scaffolds with an N50 length of 3.27 Mb, including the 12 largest scaffolds with lengths > 2 Mb (40,021,009, accounting for 86.4% of the complete genome) displayed as a circos-plot (Table 1, Fig. 4). The total assembly size was approximately 46.34 Mb and 14,091 protein-coding genes were predicted, with a gene density of 304 genes per 1 Mb. In total, 465 non-coding RNAs (ncRNAs) were predicted, including 137 tRNAs (Table 1). Furthermore, 13,188 putative protein-coding genes were supported by mapping of the RNA-seq data to the obtained genome (Wang et al. 2018b), confirming the accuracy of the full genome assembly for the LW-Hubei strain (Additional file 22: Table S6). The average length of the predicted genes was 1684 bp, containing 3.1 exons (average length, 470 bp) and 2.1 introns (average length, 116 bp). The length of the coding regions of the predicted genes accounted for 43.89% of the total genome. The genomic features of the B. dothidea LW-Hubei isolates are shown in Table 1 and the genome structures are displayed in a circos-plot (Fig. 4). The complete genome has been deposited in the DDBJ/EMBL/GenBank databases under accession numbers SACU00000000. It is also worth noting that the complete genome of the B. dothidea LW-Hubei isolate is much larger than those of the previously reported isolates of B. dothidea LW030101 and B. dothidea Sxg01s that infect apple and grapevine, respectively, with genome sizes estimated at approximately 45.26 M and 43 M, respectively, based on assemblies of 125-bp libraries (Liu et al. 2016, Yan et al. 2018).
Repetitive sequences and potential methylation sites
Repetitive sequences and transposable elements (TEs) increase the genome size and play vital roles in the genetic evolution and classification of fungal species (Bodega and Orlando 2014; Peter et al. 2016; Krishnan et al. 2018). A total of 4,309,902 bp (accounting for 9.3% of the genome) of repetitive sequences were identified in the B. dothidea LW-Hubei isolate genome, including DNA transposons, tandem repeat (TRs) sequences, long terminal repeats (LTR), short interspersed elements (SINE) and other unclassified transposons by using four types of prediction software. Interestingly, 94.9% (4,089,081 of 4,309,902) of the repetitive sequences were TEs, whereas the TRs accounted for only 5.1%. Notably, LTR retrotransposons and SINEs accounted for 52.2% (2,132,909 of 4,089,081) and 29.7% (1,216,468 of 4,089,081) of all TEs, respectively. In addition, DNA transposons accounted for a relatively low percentage of the sequence at 5.3% (216,254 of 4,089,081), and unknown TEs comprised 14.8% (604,268 of 4,089,081) (Table 2).
It was recently reported that DNA methylation may exist in fungi (Luo et al. 2018). In total, 374,995 m4C (4-methyl-cytosine) residues, and 20,121 m6A (6-methyl-adenosine) were identified by SMRT analysis (Fig. 4). Moreover, most of the methylation sites (1,117,684) were unclassified. In accordance with the identified methylation positions, multiple motifs recognized specifically by transmethylases were detected (Additional file 23: Table S7).
Comparative genomic analysis of B. dothidea LW-Hubei strains in Botryosphaerace
Synteny analysis of the B. dothidea LW-Hubei isolate genome revealed that B. dothidea LW-Hubei exhibits a high degree of synteny with the other four Botryosphaeriaceae spp. genomes, Macrophomina phaseolina, Neofusicoccum parvum, Diplodia seriata and Diplodia corticola (Islam et al. 2012; Blanco-Ulate et al. 2013; van der Nest et al. 2014; Yan et al. 2018), at the nucleotide and amino acid levels, indicating sharing of the conserved core genes of Botryosphaeriaceae (Additional file 7: Figure S7). In detail, the scaffold1 sequences of the B. dothidea LW-Hubei isolate corresponded well with the corresponding scaffolds of M. phaseolina, N. parvum, D. seriata and D. corticola (Fig. 5a). Venn diagrams of the core-pan genes showed the existence of 3292 core genes in the five Botrysphaerace species, in addition to 4202 genes unique to the LW-Hubei isolate, 5294 for M. phaseolina, 3405 for N. parvum, 2384 for D. corticola, and 1412 for D. seriata (Additional file 8: Figure S8). Heatmaps of the dispensable genes in each strain showed the similarities and clusters of genes of the five species (Additional file 9: Figure S9).
Phylogenetic trees were constructed based on the similarities and differences in genotypes between species to reveal the relationships between fungal species in terms of genetic evolution. In addition to the reported sequences of four botryosphaeriaceous isolates, four Dothideomycetes spp. genomes were also included in comparative analyses (Additional file 24: Table S8). The sequences of the core-pan genes and sing-copy orthologs revealed that the LW-Hubei isolate is evolutionarily close to M. phaseolina, a plant pathogen that can cause stem bark sclerotia disease in jute plants (Fig. 5b). In addition, LW-Hubei is also evolutionarily close to the other three Botrysphaerace strains. Homologous proteins of LW-Hubei showed an average identity of 86.7, 83.3, 79.7, and 82.4% with those of M. phaseolina, N. parvum, D. corticola and D. seriata, respectively.
Gene families can be used to analyze the evolutionary history and differentiation paths of genes. Comparisons of the numbers of orthologs showed that LW-Hubei, M. phaseolina, N. parvum, D. corticola and D. seriata shared 4301 gene families, including 3139 single-copy orthologs, which were common to all five species. In total, 229 Botryosphaeriaceae spp.-specific groups were identified, consisting of 95 families unique to the LW-Hubei strain, 94 for M. phaseolina, 23 for N. parvum, 13 for D. corticola and four for D. seriata (Fig. 6). When compared with all other genomes, LW-Hubei had the highest orthologous (single-copy and multiple-copy orthologs) similarity to M. phaseolina (93.15%, 7334 of 7873), while sharing 90.98% (7163 of 7873), 85.84% (6758 of 7873) and 82.1% (6464 of 7873) gene orthologue similarities with N. parvum, D. corticola and D. seriata, respectively. The phylogenetic tree analysis showed that LW-Hubei is closely related to M. phaseolina based on gene family analysis (Additional file 10: Figure S10).
Genome annotation of B. dothidea LW-Hubei isolate
A total of 14,091 protein-coding genes were predicted to exist in the genome of the B. dothidea LW-Hubei strain, covering approximately 43.9% of the genome sequence. Of these, 13,135 genes were annotated by comparisons with multiple databases. Functional annotation showed that 7558 (53.6%) genes were assigned with GO terms, most of which were involved in metabolic process, catalytic activity, binding, cell process and response to stress or stimulus. A further 4523 (32.1%) genes were mapped to the KEGG pathway database (Additional file 25: Table S9 and Additional file 11: Figure S11). The major functional groups included transcription factors (573 genes), pathogen-host interaction (PHI) (1096 genes) and 351 putative CAZymes, of which 79 (22%) potential CAZymes were annotated as PHI, indicating that a large group of PHI may be involved in LW-Hubei pathogenicity (Additional file 26: Table S10). The CAZy categories contained 59 auxiliary activities (AA), including 13 AA1 with degrading lactase activity, 85 carbohydrate binding modules (CBMs), 16 carbohydrate esterases (CEs), 131 glycoside hydrolases (GHs), 44 glycosyl transferases (GTs) and 16 polysaccharide lyases (PLs). Compared with the other four botryosphaeriaceous fungi, the B. dothidea LW-Hubei isolate contained a larger group of potential GHs, followed by CBMs, including CBM1, CBM12, CBM13, CBM18, CBM20, CBM21, CBM39, CBM42, CBM47, CBM48, CBM52 and CBM63. The numbers of CAZymes possessed by B. dothidea LW-Hubei was equal to that of M. phaseolina (362) (Islam et al. 2012).
Our analyses predicted 3833 secreted proteins. In addition, 12 putative effectors were annotated using the PHI database (Additional file 27: Table S11). The corresponding numbers of secreted proteins found in M. phaseolina (1863), D. seriata (910) and N. parvum (1097) were lower than that of the B. dothidea LW-Hubei isolate.
The B. dothidea genome encodes 552 transporter genes comprising 18 families. The majority of the transporter genes (133 of 552) were similar to those cataloged in the PHI-base. In the B. dothidea LW-Hubei genome, the largest proportion of transporters belonged to the major facilitator superfamily (MFS) family (83 genes), followed by the MC family (29), the NDH family (29), the ABC superfamily (27) and the amino acid-polyamine-organocation (APC) family (24) (Additional file 28: Table S12).
In addition, the B. dothidea genome encodes 128 kinase genes, consisting of AGC (13), atypical (14), CAMK (19), CK1 (2), CMGC (27), STE (18) and others (24) (Additional file 29: Table S13), which had orthologs in the PHI-base (82 of 128). This indicates that transporters and kinases in B. dothidea play a functional role in the pathogen-host interaction (Saier Jr et al. 2014; Cheng et al. 2015).
The identification and expression levels of gene silencing components in B. dothidea strains during periods of infection with mycovirus
According to the functional annotation of the full genome in combination with de novo sequencing analysis reports (Wang et al. 2018b), several RNAi components were identified in the B. dothidea LW-Hubei genome, including two Dicers (DCL1 for GME1231_g and DCL2 for GME553_g), four argonautes (GME11353_g, GME6306_g, GME9732_g and GME8235_g) and three RNA-directed RNA polymerases (RdRp1 for GME8162_g, RdRp2 for GME11473_g, RdRp3 for GME9357_g).
The BdDicer1 ORF encodes a protein of 1895 amino acids, and the size of the predicted protein encoded by BdDicer2 is 1519 amino acids. Both BdDicers possess a DEAD box, RNA helicase domain (helC), Dicer-dimerization domain and two RNaseIII catalytic domains. Amino acid sequence analysis revealed that three BdRdRps contained an RNA-dependent RNA polymerase domain (RdRp), which are highly conserved in fungi, whereas BdRdRp1 and BdRdRp2 possess an additional RNA recognition motif (RRM) in the N-terminal region. Three BdAgo proteins (GME11353_g, GME6306_g, GME9732_g) containing three conserved domains (Argo L1, PAZ and Piwi), and one BdAgo (GME8235_g) with two conserved domains (PAZ and Piwi) were identified among other argonaute orthologs (Fig. 7a).
To explore the molecular genetic evolution of the RNA-silencing pathway in Ascomycota, phylogenetic tree analysis of RdRp-, Dicer- and Ago-like proteins from B. dothidea LW-Hubei was performed at the amino acid level using the sequences of corresponding proteins in other fungal species from the Ascomycota clade as references (Additional file 30: Table S14). The BdDicers, BdRdRps and BdAgos clustered clearly with known members from other fungi with functions in meiotic silencing via the unpaired DNA (MSUD) or quelling pathways (Romano and Macino 1992; Fulci and Macino 2007; Li et al. 2010; Dang et al. 2011). The two BdDicer proteins were grouped into two clusters, whereas the four BdAgo proteins and three BdRdRp proteins were clustered into three groups. BdDicer2, BdRdRp3, and BdAgo1,3 clustered in a group with their counterparts from the quelling pathway, which were homologous to Neurospora crassa dcl2 (27.86% identity), qde1 (14.13%), and Qde2 (30.2 and 29.29%), respectively. BdDicer1, BdRdRp1, and BdAgo2 clustered with proteins from the MSUD pathway, which were homologous to Neurospora crassa sms3 (34.76% identity), sad1 (30.19%), and sms2 (32.71%), respectively. BdAgo4 was identified within a single group (Fig. 7b). BdRdRp2 was identified in a separate group with N. crassa RdRp3 (29.64% identity), the function of which is unknown. The expression level of BdRdRp2 (GME11473_g) in mycelia was elevated in response to BdPV1 infection of LW-P and LW-CP strains for 10 d (Fig. 7c).
The expression levels of RNAi components from B. dothidea LW-Hubei strains in response to BdCV1 and BdPV1 were determined by RT-qPCR. The results showed that the expression levels of BdDicer, BdAgo and BdRdRp proteins in mycelia were upregulated to some extent in response to BdCV1 and BdPV1 infection for 5 d and 10 d, respectively (Fig. 7c, Additional file 17: Table S1). Taken together, these results show that BdCV1 and BdPV1 infection of B. dothidea strains elicits the accumulation of BdDicers, BdAgos, and BdRdRps expression at the mRNA level, indicating that these genes are involved in the antiviral silencing defense pathway in B. dothidea following infection with mycovirus.
Analysis of small RNAs of B. dothidea strains in response to dsRNA mycovirus infection
The expression of RNAi genes in B. dothidea led to the hypothesis that it possesses functional small RNAs (sRNAs) that trigger gene silencing. Four sRNA libraries were constructed from B. dothidea strains infected with mycovirus and a Mock (LW-Hubei strain) infected control and subjected to Illumina deep sequencing to investigate the existence and expression of sRNAs in B. dothidea strains; the raw sRNA sequencing data have been deposited in the NCBI sequence read archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) under the accession numbers SRP174595 (BioProject Acc. no. PRJNA51169). Fungal-specific reads were identified by mapping to the obtained B. dothidea LW-Hubei genome. A total of 11,221,049, 11,837,018, 10,949,651 and 10,880,651 clean sRNAs reads were obtained from LW-CP, LW-C, LW-P and Mock, respectively. Annotation of the sequenced clean sRNAs, including intron, exon, repeats, rRNA, tRNA, snRNA, snoRNA, miRNA sequences and other unannotated reads are shown in Table 3. Comparison of sRNA clean reads from the three LW-CP/Mock, LW-C/Mock and LW-P/Mock comparisons libraries revealed great differences in reads among the LW-CP, LW-C, LW-P and Mock samples, although the expression of the common sRNA sequences is centralized as previously reported in plant and fungal hosts (Additional file 31: Table S15). Furthermore, these comparisons revealed the high level of consistency in the high-throughput sequencing data obtained for the four B. dothidea strains.
B. dothidea miRNA-like RNAs (designated as Bd-milRNAs) accounted for 264,274 (2.36%), 222,946 (1.94%), 267,681 (2.44%) and 278,093 (2.56%) of the sequences in the LW-CP, LW-C, LW-P and Mock libraries, respectively (Table 3). Analysis of the length distribution revealed that sRNAs consisting of 21 nt were the most abundant, with 2,569,298 (191,333 unique reads), 3,582,770 (228,004 unique reads) and 3,816,137 (239,297 unique reads) in the LW-CP, LW-P and Mock libraries, respectively, while sRNAs consisting of 25 nt were the most abundant in LW-C, with 2,343,056 reads (88,608 unique reads). The next most abundant classes were the 20-nt and 22-nt sRNAs from the LW-CP, LW-P and Mock libraries, respectively, while the next most abundant class was 24-nt in LW-C (Fig. 8a). Among the unique sequences, the 21-nt sRNA sequences were most abundant, followed by the 20-nt sequences in the LW-P and Mock libraries, and the 22-nt sRNAs in the LW-CP library. The 22-nt sRNAs were the most abundant in the LW-C library (Fig. 8b).
Among the known Bd-milRNAs, the most abundant were 21-nt in length, with a total of 83,676, 73,641, 130,666 and 101,313 reads in the LW-CP, LW-C, LW-P and Mock libraries, respectively (Fig. 8c). Among the novel Bd-milRNAs, the most abundant were 21-nt in length, with a total of 209, 311, 1192 and 3623 reads in the LW-CP, LW-C, LW-P and Mock libraries, respectively. These data revealed inhibition of the Bd-milRNAs in the LW-CP and LW-C libraries (Fig. 8c). There were marked differences in the overall size distribution patterns of sRNAs between the libraries generated from the hypovirulent and virulent B. dothidea strains infected with mycoviruses, suggesting that BdCV1 and BdPV1 differentially affect sRNA and Bd-milRNAs accumulation in B. dothidea. In addition, these findings indicate that the proportion of reads from libraries of LW-CP, LW-C and LW-P, do not match the obtained target genome sequences at high proportion, mainly because of mycoviral infection of B. dothidea strains (Additional file 22: Table S6).
Prediction novel Bd-milRNAs candidate in responsive to mycovirus
After removing the Bd-milRNAs with expression levels < 5 reads for further analysis of differential expression in B. dothidea strains, a total of 68 novel Bd-milRNAs candidate were predicted. Of these, 19, 21, 16 and 19 Bd-milRNAs were identified from the LW-CP, LW-C, LW-P and Mock libraries, respectively (Fig. 9b, Additional file 32: Table S16). The precursor sequences of these novel Bd-milRNA candidates were folded into a hairpin-like structure with precursor lengths ranging from 67 to 356 nt and minimal folding energy (MFE) ranging from − 21.4 to − 136.4 kal/mol (Additional file 32: Table S16 and Additional file 12: Figure S12). The base bias in the first position among the predicted novel Bd-milRNA candidates showed a preference for uridine (U) (Additional file 13: Figure S13). The sequencing and bioinformatics results showed lower expression levels of the novel Bd-milRNAs.
Differential expression analysis of Bd-milRNAs from B. dothidea strains in response to mycoviral infection
The role of Bd-milRNAs in fungus-mycovirus interactions was investigated. A total of 167 known Bd-milRNAs were predicted, of which 64, 48, 83 and 86 were identified from the LW-CP, LW-C, LW-P and Mock libraries, respectively. Furthermore, 110, 120 and 88 known Bd-milRNAs were found to be differentially expressed in analysis of the three LW-CP/Mock, LW-C/Mock and LW-P/Mock comparison libraries, respectively (Fig. 9a, c). As expected, more changes in Bd-milRNA expression were observed in the analysis of the LW-CP/Mock and LW-C/Mock comparison libraries than those in the LW-P/Mock, revealing that BdCV1 has a greater effect than the BdPV1 on the expression of Bd-milRNAs. Evaluation of the LW-CP/Mock, LW-C/Mock and LW-P/Mock comparison libraries revealed common differential expression of 57 known Bd-milRNAs, while 12 known Bd-milRNAs were specifically expressed in the LW-CP/Mock library, 29 in the LW-C/Mock library and 20 in the LW-P/Mock library. In addition, 161 known Bd-milRNAs were differentially expressed, with the majority showing up-regulated expression in response to mycoviruses (Additional file 33: Table S17). Four randomly selected known Bd-milRNAs, Bd-milR8635, Bd-milR172c-3p, Bd-milR5636 and Bd-milR5568f-3p, were verified by stem-loop RT-PCR and sequencing analyses (Additional file 18: Table S2 and Additional file 34: Table S18).
Further, differential expression analysis of predicted novel Bd-milRNAs was performed. Twenty novel Bd-milRNAs candidates were expressed differentially, of which 12, 14 and 11 were significantly differentially expressed in the LW-CP/Mock, LW-C/Mock and LW-P/Mock libraries, respectively. Six novel Bd-milRNAs were commonly expressed differentially in the LW-CP/Mock, LW-C/Mock and LW-P/Mock libraries, while three novel Bd-milRNAs were specifically expressed in the LW-CP/Mock library, three in the LW-C/Mock library and three in the LW-P/Mock library (Fig. 9b, d). This revealed that specific and common alterations in Bd-milRNAs accumulation levels were induced in response to infection of B. dothidea strains with mycoviruses exhibiting mild and severe hypovirulence traits. The expression of two novel (Bd-milR6, Bd-milR12) and one known Bd-milRNA (Bd-milR1147.2) selected at random was verified by poly (A) tailing-mediated RT-PCR analysis, meanwhile u6 was as control (Additional file 18: Table S2 and Additional file 14: Figure S14).
Temporal expression patterns of potential milRNA-mediated target genes in response to mycovirus infection of B. dothidea strains
Based on the B. dothidea LW-Hubei genome sequence, 690 targets were predicted as target genes for 81 known Bd-milRNAs from the LW-C/Mock library by both the psRobot and TargetFinder software, while 604 target genes for were predicted for 73 known Bd-milRNAs from the LW-CP/Mock library and 568 target genes for 62 known Bd-milRNA from the LW-P/Mock library. It was observed that many Bd-milRNAs regulated multiple mRNA targets from B. dothidea (Additional file 35: Table S19). In addition, using the psRobot and/or TargetFinder software, 35 target genes were predicted for 11 novel milRNAs from the LW-C/Mock library, 51 target genes for nine novel milRNAs from the LW-CP/Mock library and 25 target genes for nine novel milRNAs from the LW-P/Mock library (Additional file 36: Table S20). The targeted genes were annotated as TFs, transporters, kinases and effectors involved in the PHI, response to stimulus and signal transduction pathways (Additional file 35: Table S19 and Additional file 36: Table S20).
To confirm the ability of mycovirus-responsive Bd-milRNAs to regulate expression of their target genes during the period of infection, the expression levels of randomly selected Bd-milRNAs and putative target mRNAs were analyzed in B. dothidea strains cultured on PDA for 5 d. RT-qPCR was performed to analyze the temporal expression patterns of the following putative target mRNAs mediated by Bd-milRNAs: putative calcium-transporting p-type ATPase (GME3643_g) and hypothetical protein MPH_03910 (GME1913_g) as predicted targets of Bd-milR12; peptidase M43 pregnancy-associated plasma-A (GME4906_g) as a predicted target of Bd-milR33; putative polyketide synthase protein (GME2439_g) and hypothetical protein MPH_01630 (GME5002_g) as predicted targets of Bd-milR60; putative ABC transporter protein ATP-dependent bile acid permease (GME4334_g) and multicopper oxidase type 1 (GME2445_g) as predicted targets of Bd-milR1147.2; zinc finger C2H2-type protein (GME4116_g) as a predicted target of Bd-milR2094-5p; glucose-6-phosphate dehydrogenase (GME4518_g) as a predicted target of Bd-milR46; mannose-6-phosphate receptor binding protein (GME2361_g) as a predicted target of Bd-milR6; sugar/inositol transporter mfs sugar transporter (GME11656_g) as a predicted target of Bd-milR65; putative wd repeat domain 5b (GME3697_g) as a predicted target of Bd-milR172c-3p; urea carboxylase (GME1751_g) and putative ABC multidrug transporter mdr1 protein (GME4636_g) as predicted targets of Bd-milR8635; The cyclin-dependent protein kinase complex component PHO85 cyclin-6 (GME10476_g), hypothetical protein MPH_03537 (GME12854_g) and BdCV1 RdRp (GenBank Acc. no, KF688736) as predicted targets of Bd-milR5636; and BdCV1 CP (GenBank Acc. no, KF688737) as a predicted target of Bd-milR5568. The correlations between the expression levels of these Bd-milRNAs and their target mRNAs during mycelial growth for 5 d following BdCV1 and BdPV1 infection are shown in Table 4 The results revealed that the expression levels of partial target mRNAs exhibited the opposite expression patterns compared with those of the corresponding Bd-milRNAs in B. dothidea strains in response to mycovirus infection. In contrast, the expression patterns of partial putative target genes showed a positive correlation with those of the corresponding Bd-milRNAs (Table 4, Additional file 18: Table S2 and Additional file 15: Figure S15).
Functional enrichment of putative milRNA-mediated target genes
The functional roles of the predicted milRNA target genes in the three comparison libraries were predicted based on analysis of the functional enrichment of all identified Bd-milRNA targets using WEGO (Web Gene Ontology Annotation Plot) software (Fig. 10). The predicted milRNA-regulated targets showed enrichment of GO terms in the biological process, cell process and molecular function categories, which is accordance with the predictions of de novo transcription analysis described previously (Wang et al. 2018b). In the biological process category, the predicted targets were involved in different biological processes, including cellular metabolic process, biosynthetic process, localization, catabolic process, cellular response to stimulus, regulation of cellular process, regulation of biological process, and signal transduction. In the cellular component category, the enriched GO terms included cell part, membrane-bounded organelle, membrane part, hydrolase activity and transferase activity. In the molecular function category, the enriched GO terms were involved in binding, cofactor binding, carbohydrate binding and response to stress. In addition, pathway enrichment analyses revealed spliceosome (ko03040), including 57 target genes with pathway annotation, was enriched to a greater level in B. dothidea strains in response to BdCV1 infection (Additional file 16: Figure S16).
Among the economically important Botryosphaeriae species of fungi, B. dothidea causes severe warts and stem canker, as well as ring rot diseases of pear worldwide (Slippers and Wingfield 2007; Zhai et al. 2014). In this study, we investigated the potential pathogenic mechanism of B. dothidea LW-Hubei isolates on pear and the interactions of B. dothidea and mycovirus. To address these issues, we performed high-throughput genome sequencing in combination comparative genomics analysis of sRNAs, including Bd-milRNA/target mRNAs from B. dothidea strains infected with dsRNA mycoviruses.
Our investigations revealed that the complete genome sequence size of the B. dothidea LW-Hubei isolate is 46.34 Mb (Table 1). This is larger than the genome sizes of sequenced but not annotated Botryosphaeriaceae taxa B. dothidea and N. parvum isolates with similar genome sizes (42.6 Mb to 45.26 Mb), and low proportions of repeat sequences determined using the NGS method (Bodega and Orlando 2014; van der Nest et al. 2014; Liu et al. 2016; Luo et al. 2018; Yan et al. 2018). In particular, the assembled B. dothidea LW-Hubei isolate genome contains 4.31 Mb repetitive sequences, representing a relatively high proportion (9.41%) of the complete genome in Botryosphaeriaceae (Table 2, Fig. 4). The repeat sequences were randomly distributed across the full genome of B. dothidea LW-Hubei isolate. It was recently reported that DNA methylation may be involved in host cell processes, such as regulation gene transcription in fungal species (Luo et al. 2018; Plissonneau et al. 2018; Yan et al. 2018). The m6A and m4C in B. dothidea were predicted for the first time using SMRT, which revealed larger numbers of m4C in the genome, occurring at a frequency of 8082/Mb, while m6A was detected at a frequency of 434/Mb. However, the mechanism of m4C and m6A regulation in B. dothidea remains to be elucidated.
The comparative genome analysis revealed gene families and core-pan genes that are specific to B. dothidea LW-Hubei as well as those that are common to Botryosphaeriaceae spp., and are more closely related to M. phaseolina than to D. corticola and D. seriata (Figs. 5, and 6). Based on the genome assembly and annotation, we analyzed the specific CAZymes, effectors protein, kinase, and transporter proteins involved in B. dothidea infection of pear, and showed that there are slightly fewer than those found in other species of the Botryosphaeriaceae taxa (Islam et al. 2012; Fernandes et al. 2014; van der Nest et al. 2014; Liu et al. 2016; Wang et al. 2018b; Yan et al. 2018), although B. dothidea was found to possess a significant number of genes involved in PHI. Annotation of these genes revealed that the B. dothidea genome possesses a large number genes involved in the pathogenesis in pear hosts that are both specific and common to those of other pathogenic fungi. These findings significantly expand our understanding of the determinants of B. dothidea pathogenicity and host-pathogen interactions.
It is worth noting the existence of RNAi in the B. dothidea LW-Hubei strain. Genes encoding two homologs of Dicer (DCL1 and DCL2), four homologs of Argonaute, and three homologs of RdRp, were identified in the genome (Fig. 7). The phylogenetic trees indicate that BdDCL2, BdRdRp3, and BdAgo1,3 play roles in the RNAi pathway, whereas the combination of BdDCL1, BdRdRp1, and BdAgo2 might be required for the MSUD pathway in B. dothidea (Fig. 7b). In particular, RT-qPCR analysis showed that the expression of RNAi components was unregulated in response to BdCV1 and BdPV1 to some extent (Fig. 7c), which is also supported by the results of de novo transcriptome sequencing of B. dothidea strains (Wang et al. 2018b). It demonstrated that knock-out BdDicer2 gene led to reduced mycelia growth, virulence and BdPV1 cp and RdRp mRNA expression in transformants in comparison with the wild strain LW-P (data not published). It revealed that RNAi was related to regulation of the pathogenesis of B. dothidea and interaction with mycovirus. In addition, the horizontal transmission between the LW-C/BdCV1 and LW-P/BdPV1 strains of B. dothidea revealed that BdPV1 can not be transmitted from LW-P to LW-C, indicating that BdCV1 might induce gene silencing that leads to the resistance of LW-C to the transmission of BdPV1 (Fig. 3). Based on these results, we hypothesize that B. dothidea possesses a RNA-silencing system preferentially involving RNAi as part of the defense mechanism against mycovirus infection.
RNA-silencing involves the function of sRNAs as small regulatory molecules that further control endogenous and exogenous RNAs in mycoviruses and regulate biological processes (Zhang et al. 2008; Li et al. 2010; Dang et al. 2011; Huntzinger and Izaurralde 2011). The RNA-silencing components encoded by the Dicer, Ago, RdRp genes have been reported to influence sRNA biosynthesis. RNA-silencing involving ChDCL1 and ChAGO1 in C. higginsianum is proposed to function as an antiviral mechanism, and FgAGO-1 and FgDICER-2 is responsible for hairpin RNA-triggered RNA-silencing and related small interfering RNA accumulation in F. graminearum (Sun et al. 2009; Huntzinger and Izaurralde 2011; Nuss 2011; Chen et al. 2015; Campo et al. 2016). It can be speculated that a sophisticated RNA-silencing system that regulates sRNA production exists in B. dothidea strains.
It has been reported that milRNA expression profiles change during mycelium growth, conidiogenesis, and pathogenicity of fungi (Campo et al. 2016; Wang et al. 2016a, 2016b; Wang et al. 2017a, 2017c). sRNAs, especially miRNAs and siRNAs, are perturbed during mycovirus infection of fungi and have important influences on biological processes (Zhang et al. 2014; Campo et al. 2016; Wang et al. 2016a, 2016b). In this study, we revealed that BdCV1 has a greater effect on the biological characteristics of B. dothidea than BdPV1 (Figs. 1, 2 and 3, Additional file 1: Figure S1). Based on previous reports of milRNA-mediated gene silencing in fungi, we investigated the existence of milRNAs in B. dothidea by constructing and sequencing four sRNA libraries. Additionally, the genome data obtained provided a convenient resource for genome-scale investigation of Bd-milRNAs where context sequences around sRNA loci were needed for secondary structure predication. The results of the size distribution, position-specific nucleotide preference, and accumulation of specific sequences all suggest that B. dothidea possesses an endogenous sRNA biogenesis pathway. This is the first report of alteration in milRNAs in hypovirulent and virulent B. dothidea strains infected with mycovirus.
These bioinformatics and RT-qPCR analyses of changes in the expression of the B. dothidea transcriptome and milRNAs revealed that the expression of Bd-milRNAs was significantly altered by mycovirus infection, especially BdCV1, which has a great effect on the expression of B. dothidea sRNAs, including Bd-milRNA. These findings indicate that such dynamic sRNA (milRNA) profiles may have important functions in suppressing mycoviral infection (Table 3, Figs. 8 and 9). Furthermore, our results indicate the hypovirulent and non-hypovirulent phenotypes and the pathogenity of B. dothidea strains caused by mycovirus infection are related to the numbers of differentially expressed milRNAs and milRNA-regulated mRNAs (Wang et al. 2018b). Previous studies have demonstrated that miRNAs regulate key genes in disease resistance pathways to affect viral infections (Wang et al. 2016a, 2016b; Wang et al. 2017a, 2017c; Fu et al. 2017). In this study, we identified a number of predicted milRNA-target genes from B. dothidea that are possibly involved in disease resistance and defense. For instance, downregulation of the expression of an ABC transporter protein that was identified as a predicted target of Bd-milR1147 in LW-C infected with BdCV1 may inhibit the B. dothidea infection. In additions, expression of the MFS sugar transporter (GME11656_g) targeted by Bd-milR65 was up-regulated in B. dothidea in response to mycovirus infection. These transporters are important in the pathogenesis of fungi (Saier Jr et al. 2014); thus, our findings indicate the involvement of such transporters in the pathogenicity of B. dothidea-infecting mycoviruses. Glucose-6-phosphate dehydrogenase (GME4518_g) was upregulated in response to BdCV1 in the B. dothidea strains LW-C and LW-CP. Zinc finger C2H2-type protein (GME4116_g) as pathogenic factor was up-regulated, targeted by down-regulated expression Bd-milR2094-5p in response to BdCV1 and/or BdPV1 infection of B. dothidea strains, which may involve in the interaction of mycovirus and B. dothidea (Table 4).
In addition, BdCV1 RdRp and CP were predicted as targets of Bd-milRNA5636 and Bd-milRNA5568, respectively. It has been reported that gga-miR-130b-3p plays a crucial role in host defense against IBDV (Infectious Bursal Disease Virus) infection by targeting the viral itself genome (Cheng et al. 2015; Fu et al. 2017); thus, in accordance with this, it can be speculated that Bd-milRNAs may be directly involved in the interaction of mycovirus with the B. dothidea host. WEGO analysis revealed that the predicted Bd-milRNA-targeted genes were enriched for kinases, small secreted proteins and transporters (Fig. 10). The majority of these targets have roles in stress-related gene regulation or PHI pathways, suggesting that milRNA-mediated regulation of the transcriptome is an essential process in the mycovirus-stress response in the B. dothidea host. The mechanistic details of this process in resistance of B. dothidea to BdPV1 and BdCV1 remain to be elucidated.
In this study, we sequenced the 46.34 M genome of the B. dothidea LW-Hubei strain and annotated 13,135 genes to reveal the pathogenic mechanism. We identified RNAi components and mycovirus-altered milRNA expression patterns in the B. dothidea host and showed that BdCV1 has more effect on milRNA expression than BdPV1. This information not only contributes to a more comprehensive understanding of the complex processes involved in the regulation of the biological characteristics of B. dothidea, but also provides insights into the molecular mechanisms of hypovirulence of B. dothidea in pear following mycovirus infection. Our analysis revealed that the interaction of B. dothidea and mycovirus involves the combined action of the antiviral gene silencing pathway, and milRNA-mediated regulation of target gene mRNA expression in B. dothidea.
- B. dothidea miRNA-like RNAs:
Botryosphaeria dothidea chrysovirus 1
Botryosphaeria dothidea partitivirus 1
Cluster of Orthologous Groups
Cetyltrimethyl ammonium bromide
RNA-dependent RNA polymerase
Single Molecule Real-Time
transmission electron microscopy
Web Gene Ontology Annotation Plot
Adhikari S, Curtis PD (2016) DNA methyltransferases and epigenetic regulation in bacteria. FEMS Microbiology Reviews 40:575–591
Anagnostakis SL (1982) Biological control of chestnut blight. Science 215:466–471
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H et al (2000) Gene ontology: tool for the unification of biology. Nature Genetics 25:25–29
Audic S, Claverie JM (1997) The significance of digital gene expression profiles. Genome Research 7:986–995
Badouin H, Hood ME, Gouzy J, Aguileta G, Siguenza S et al (2015) Chaos of rearrangements in the mating-type chromosomes of the anther-smut fungus Microbotryum lychnidis-dioicae. Genetics 200:1275–1284
Benson G (1999) Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Research 27:573–580
Blanco-Ulate B, Rolshausen P, Cantu D (2013) Draft genome sequence of Neofusicoccum parvum isolate UCR-NP2, a fungal vascular pathogen associated with grapevine cankers. Genome Announcements 1:e00339–e00313
Bodega B, Orlando V (2014) Repetitive elements dynamics in cell identity programming, maintenance and disease. Current Opinion in Cell Biology 31:67–73
Bologna NG, Voinnet O (2014) The diversity, biogenesis, and activities of endogenous silencing small RNAs in Arabidopsis. Annual Review of Plant Biology 65:473–503
Cai Q, Qiao L, Wang M, He B, Lin FM et al (2018) Plants send small RNAs in extracellular vesicles to fungal pathogen to silence virulence genes. Science 360:1126–1129
Campo S, Gilbert KB, Carington JC (2016) Small RNA-based antiviral defense in the phytopathogenic fungus Colletotrichum higginsianum. PLoS Pathogens 12:e1005640
Chen CF, Ridzon DA, Broomer AJ, Zhou ZH, Lee DH et al (2005) Real-time quantitative of microRNAs by stem-loop RT-PCR. Nucleic Acids Research 33:e179
Chen L, Zheng D, Liu B, Yang J, Jin Q (2016) VFDB 2016: hierarchical and refined dataset for big data analysis-10 years on. Nucleic Acids Research 44:D694–D697
Chen R, Jiang N, Jiang Q, Sun X, Wang Y et al (2014) Exploring microRNA-like small RNAs in the filamentous fungus Fusarium oxysporum. PLoS One 9:e104956
Chen Y, Gao Q, Huang M, Liu Y, Liu Z et al (2015) Characterization of RNA silencing components in the plant pathogenic fungus Fusarium graminearum. Scientific Reports 5:12500
Cheng Y, Wang X, Yao J, Voegele RT, Zhang Y et al (2015) Characterization of protein kinase PsSRPKL, a novel pathogenicity factor in the wheat stripe rust fungus. Environmental Microbiology 17:2601–2617
Chiba S, Lin YH, Kondo H, Kanematsu S, Suzuki N (2013) Effects of defective interfering RNA on symptom induction by, and replication of, a novel Partitivirus from a Phytopathogenic fungus, Rosellinia necatrix. Journal of Virology 87:2330–2341
Chiba S, Salaipeth L, Lin YH, Sasaki A, Kanematsu S et al (2009) A novel bipartite double-stranded RNA mycovirus from the white root rot fungus Rosellinia necatrix: molecular and biological characterization, taxonomic considerations, and potential for biological control. Journal of Virology 83:12801–12812
Chin CS, Alexander DH, Marks P, Klammer AA, Drake J et al (2013) Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data. Nature Methods 10:563–569
Cohen NR, Ross CA, Jain S, Shapiro RS, Gutierrez A et al (2016) A role for the bacterial GATC methylome in antibiotic stress survival. Nature Genetics 48:581–586
Dahlmann TA, Kück U (2015) Dicer-dependent biogenesis of small RNAs and evidence for MicroRNA-like RNAs in the penicillin producing fungus Penicillium chryspgenum. PLoS One 10:e0125989
Dai X, Zhao PX (2011) psRNATarget: a plant small RNA target analysis server. Nucleic Acids Research 39:W155–W159
Dang Y, Yang Q, Xue Z, Liu Y (2011) RNA interference in fungi: pathways, functions, and applications. Eukaryotic Cell 10:1148–1155
Davis BM, Chao MC, Waldor MK (2013) Entering the era of bacterial epigenomics with single molecule real time DNA sequencing. Current Opinion in Microbiology 16:192–198
Dean RA, Talbot NJ, Ebbole DJ, Farman ML, Mitchell TK et al (2005) The genome sequence of the rice blast fungus Magnaporthe grisea. Nature 434:980–986
Dong Y, Li Y, Zhao M, Jing M, Liu X et al (2015) Global genome and transcriptome analyses of Magnaporthe oryzae epidemic isolate 98-06 uncover novel effectors and pathogenicity related genes, revealing gene gain and lose dynamics in genome evolution. PLoS Pathogens 11:e1004801
Edgar RC (2004a) MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics 5:113
Edgar RC (2004b) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research 32:1792–1797
Elbourne LD, Tetu SG, Hassan KA, Paulsen IT (2017) TransportDB 2.0: a database for exploring membrane transporters in sequenced genomes from all domains of life. Nucleic Acids Research 45:D320–D324
Fahlgren N, Carrington JC (2010) miRNA target prediction in plants. Methods in Molecular Biology 592:51–57
Faino L, Seidl MF, Datema E, van den Berg GC, Janssen A et al (2015) Single-molecule real-time sequencing combined with optical mapping yields completely finished fungal genome. MBio 6:e00936–e00915
Fang G, Munera D, Friedman DI, Mandlik A, Chao MC et al (2012) Genome-wide mapping of methylated adenine residues in pathogenic Escherichia coli using single-molecule real-time sequencing. Nature Biotechnology 30:1232–1239
Fernandes I, Alves A, Correia A, Devreese B, Esteves AC (2014) Secretome analysis identifies potential virulence factors of Diplodia corticola, a fungal pathogen involved in cork oak (Quercus suber) decline. Fungal Biology 118:516–523
Fischer M, Knoll M, Sirim D, Wagner F, Funke S et al (2007) The cytochrome P450 engineering database: a navigation and prediction tool for the cytochrome P450 protein family. Bioinformatics. 23:2015–2017
Friedlander MR, Chen W, Adamidi C, Maaskola J, Einspanier R et al (2008) Discovering microRNAs from deep sequencing data using miRDeep. Nature Biotechnology 26:407–415
Fu L, Niu B, Zhu Z, Wu S, Li W (2012) CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28:3150–3152
Fu M, Wang B, Chen X, He Z, Wang Y, Li X et al (2017) MicroRNA gga-miR-130b suppresses infectious bursal disease virus replication via targeting of the viral genome and cellular suppressors of cytokine signaling 5. Journal of Virology 92:e01646–e01617
Fulci V, Macino G (2007) Quelling: post-transcriptional gene silencing guided by small RNAs in Neurospora crassa. Current Opinion in Microbiology 10:199–203
Galperin MY, Makarova KS, Wolf YI, Koonin EV (2015) Expanded microbial genome coverage and improved protein family annotation in the COG database. Nucleic Acids Research 43:D261–D269
Garibaldi A, Bertetti D, Poli A, Gullino ML (2012) First report of fruit rot in pear caused by Botryosphaeria dothidea in Italy. Plant Disease 96:910
Ghabrial SA, Suzuki N (2009) Viruses of plant pathogenic Fungi. Annual Review of Phytopathology 47:353–384
Guo L, Li J, Li B, Zhang X, Zhou Z et al (2009) Investigation on the occurrence and chemical control of Botryosphaeria canker of apple in China. Plant Protection 35:120–123 in Chinese
Huntzinger E, Izaurralde E (2011) Gene silencing by microRNAs: contributions of translational repression and mRNA decay. Nature Reviews Genetics 12:99–110
Islam MS, Haque MS, Islam MM, Emdad ME, Halim A et al (2012) Tools to kill: genome of one of the most destructive plant pathogenic fungi Macrophomina phaseolina. BMC Genomics 13:493.
Janbon G, Maeng S, Yang DH, Ko YJ, Jung KW et al (2010) Characterizing the role of RNA silencing components in Cryptococcus neoformans. Fungal Genetics and Biology 47:1070–1080
Jiang N, Yang Y, Janbon G, Pan J, Zhu X (2012) Identification and functional demonstration of miRNAs in the fungus Cryptococcus neoformans. PLoS One 7:e52734
Johnson AD, Handsaker RE, Pulit SL, Nizzari MM, O'Donnell CJ et al (2008) SNAP: a web-based tool for identification and annotation of proxy SNPs using HapMap. Bioinformatics 24:2938–2939
Jones P, Binns D, Chang HY, Fraser M, Li W et al (2014) InterProScan 5: genome-scale protein function classification. Bioinformatics 30:1236–1240
Kadotani N, Nakayashiki H, Tosa Y, Mayama S (2004) One of the two dicer-like proteins in the filamentous Fungi Magnaporthe oryzae genome is responsible for hairpin RNA-triggered RNA silencing and related small interfering RNA accumulation. Journal of Biological Chemistry 279:44467–44474
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M et al (2008) KEGG for linking genomes to life and the environment. Nucleic Acids Research 36:D480–D484
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M (2016) KEGG as a reference resource for gene and protein annotation. Nucleic Acids Research 44:D457–D462
Kang K, Zhong J, Jiang L, Liu G, Gou CY et al (2013) Identification of microRNA-like RNAs in the filamentous fungus Trichoderma reesei by solexa sequencing. PLoS One 8:e76288
Kim KE, Peluso P, Babayan P, Yeadon PJ, Yu C et al (2014) Long-read, whole-genome shotgun sequence data for five model organisms. Scientific data 1:140045
Kim VN, Han J, Siomi MC (2009) Biogenesis of small RNAs in animals. Nature Reviews Molecular Cell Biology 10:126–139
Krishnan P, Meile L, Plissonneau C, Ma X, Hartmann FE et al (2018) Transposable element insertions shape gene regulation and melanin production in a fungal pathogen of wheat. BMC Biology 16:78
Kumar S, Stecher G, Tamura K (2016) MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Molecular Biology and Evolution 33:1870–1874
Langmead B, Trapnell C, Pop M, Salzberg SL (2009) Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology 10:R25
Lau SK, Chow WN, Wong AY, Yeung JM, Bao J et al (2013) Identification of microRNA-like RNAs in mycelial and yeast phases of the thermal dimorphic fungus Penicillium marneffei. PLoS Neglected Tropical Diseases 7:e2398
Lee HC, Li L, Gu W, Xue Z, Crosthwaite SK et al (2010) Diverse pathways generate microRNA-like RNAs and dicer-independent small interfering RNAs in fungi. Molecular Cell 38:803–814
Levasseur A, Drula E, Lombard V, Coutinho PM, Henrissat B (2013) Expansion of the enzymatic repertoire of the CAZy database to integrate auxiliary redox enzymes. Biotechnology for Biofuels 6:41
Li L, Chang SS, Liu Y (2010) RNA interference pathways in filamentous fungi. Cellular and Molecular Life Science 67:3849–3863
Lin R, He L, He J, Qin P, Wang Y et al (2016) Comprehensive analysis of microRNA-seq and target mRNAs of rice sheath blight pathogen provides new insights into pathogenic regulatory mechanisms. DNA Research 23:415–425
Liu Z, Lian S, Li B, Lu H, Dong X et al (2016) Draft genome sequence of Botryosphaeria dothidea, the pathogen of apple ring rot. Genome Announcements. 4:e01142–e01116
Livak KJ, Schmittgen TD (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods 25:402–408
Luo X, Cao J, Huang J, Wang Z, Guo Z et al (2018) Genome sequencing and comparative genomics reveal the potential pathogenic mechanism of Cercospora sojina Hara on soybean. DNA Research 25:25–37
Makarova KS, Wolf YI, Koonin EV (2015) Archaeal clusters of orthologous genes (arCOGs): an update and application for analysis of shared features between thermococcales, methanococcales, and methanobacteriales. Life 5:818–540
Manso AS, Chai MH, Atack JM, Furi L, De Ste Croix M et al (2014) A random six-phase switch regulates pneumococcal virulence via global epigenetic changes. Nature Communications 5:5055
Mueth NA, Ramachandran SR, Hubert SH (2015) Small RNAs from the wheat stripe rust fungus (Puccinia striiformis f.sp. tritici). BMC Genomics 16:718
Murray MG, Thompson WF (1980) Rapid isolation of high molecular weight plant DNA. Nucleic Acids Research 8:4321–4325
Nandi T, Ong C, Singh AP, Boddey J, Atkins T et al (2010) A genomic survey of positive selection in Burkholderia pseudomallei provides insights into the evolution of accidental virulence. PLoS Pathogen 6:e1000845
Nuss DL (2005) Hypovirulence: Mycoviruses at the fungal-plant interface. Nature Reviewers Microbiology 3:632–642
Nuss DL (2011) Mycoviruses, RNA silencing, and viral RNA recombination. Advances in Virus Research 80:25–48
Pearson MN, Beever RE, Boine B, Arthur K (2009) Mycoviruses of filamentous fungi and their relevance to plant pathology. Molecular Plant Pathology 10:115–128
Peter M, Kohler A, Ohm RA, Kuo A, Krützmann J et al (2016) Ectomycorrhizal ecology is imprinted in the genome of the dominant symbiotic fungus Cenococcum geophilum. Nature Communications 7:12662
Plissonneau C, Hartmann FE, Croll D (2018) Pangenome analyses of the wheat pathogen Zymoseptoria tritici reveal the structural basis of a highly plastic eukaryotic genome. BMC Biology 16:5.
Romano N, Macino G (1992) Quelling: transient inactivation of gene expression in Neurospora crassa by transformation with homologous sequences. Molecular Microbiology 6:3343–3353
Saier MH Jr, Reddy VS, Tamang DG, Västermark A (2014) The transporter classification database. Nucleic Acids Research 42:D251–D258
Segers GC, Zhang X, Deng F, Sun Q, Nuss DL (2007) Evidence that RNA silencing functions as antiviral defense mechanism in fungi. Proceedings of the National Academy of Sciences of the United States of America 104:12902–12906
Sit CS, Ruzzini AC, Van Arnam EB, Ramadhar TR, Currie CR et al (2015) Variable genetic architectures produce virtually identical molecules in bacterial symbionts of fungus-growing ants. Proceedings of the National Academy of Sciences of the United States of America 112:13150–13154
Slippers B, Crous PW, Denman S, Coutinho TA, Wingfield BD et al (2004) Combined multiple gene genealogies and phenotypic characters differentiate several species previously identified as Botryosphaeria dothidea. Mycologia 96:83–101
Slippers B, Wingfield MJ (2007) Botryosphaeriaceae as endophytes and latent pathogens of woody plants: diversity, ecology and impact. Fungal Biology Reviews 21:90–106
Smit AFA, Hubley R, Green P (2014) RepeatMasker Open-4.0. 2013–2015. URL http://www.repeatmasker.org
Stanke M, Diekhans M, Baertsch R, Haussler D (2008) Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics 24:637–644
Sun Q, Choi GH, Nuss DL (2009) A single Argonaute gene is required for induction of RNA silencing antiviral defense and promotes viral RNA recombination. Proceedings of the National Academy of Sciences of the United States of America 106:17927–17932
Tang W, Ding Z, Zhou ZQ, Wang YZ, Guo LY (2012) Phylogenetic and pathogenic analyses show that the causal agent of apple ring rot in China is Botryosphaeria dothidea. Plant Disease 96:486–496
Tebeest DO, Shilling CW, Riley LH, Weidemann GJ (1989) The number of nuclei in spores of three species of Colletotrichum. Mycologia 81:147–149
Ter-Hovhannisyan V, Lomsadze A, Chernoff YO, Borodovsky M (2008) Gene prediction in novel fungal genomes using an ab initio algorithm with unsupervised training. Genome Research 18:1979–1990
Torto-Alalibo T, Collmer CW, Gwinn-Giglio M (2009) The plant-associated microbe gene ontology (PAMGO) consortium: community development of new gene ontology terms describing biological processes involved in microbe-host interactions. BMC Microbiology 9(Suppl 1):S1
Tsuji M, Kudoh S, Hoshino T (2015) Draft genome sequence of cryophilic basidiomycetous yeast Mrakia blollopis SK-4, isolated from an algal mat of Naga-ike Lake in the Skarvsnes ice-free area, East Antarctica. Genome Announcements 3:e01454–e01414
van der Nest MA, Bihon W, De Vos L, Naidoo K, Roodt D et al (2014) Draft genome sequences of Diplodia sapinea, Ceratocystis manginecans, and Ceratocystis moniliformis. IMA Fungus 5:135–140
Vargas WA, Martín JM, Rech GE, Rivera LP, Benito EP et al (2012) Plant defense mechanisms are activated during biotrophic and necrotrophic development of Colletotricum graminicola in maize. Plant Physiology 158:1342–1358
Wang B, Liang X, Gleason ML, Zhang R, Sun G (2018a) Comparative genomics of Botryosphaeria dothidea and B. kuwatsukai, causal agents of apple ring rot reveals both species expansion of pathogenicity-related genes and variations in virulence gene content during speciation. IMA Fungus 9:243–257
Wang B, Sun Y, Song N, Zhao M, Liu R et al (2017a) Puccinia striiformis f.sp. tritici microRNA-like RNA 1 (Pst-milR1), an important pathogenicity factor of Pst, impairs wheat resistance to Pst by suppressing the wheat pathogenesis-related 2 gene. New Phytologist 215:338–350
Wang L, Jiang J, Wang Y, Hong N, Zhang F et al (2014) Hypovirulence of the phytopathogenic fungus Botryosphaeria dothidea: association with a coinfecting chrysovirus and a partitivirus. Journal of Virology 88:7517–7527
Wang L, Luo H, Hu W, Yang Y, Hong N et al (2018b) De novo transcriptomic assembly and mRNA expression patterns of Botryosphaeria dothidea infection with mycovirus chrysovirus 1 (BdCV1) and partitivirus 1 (BdPV1). Virology Journal 15:126
Wang L, Luo H, Wang G, Wang L (2017b) Effect of Botryosphaeria dothidea Chrysovirus 1 isolate belonging to the Chrysoviridae family on growth and pathogenity of the B. dothidea strain infection in pears. Journal of Fruit Science 34:1330–1339 in Chinese
Wang M, Weiberg A, Dellota E Jr, Yamane D, Jin H (2017c) Botrytis small RNA Bc-siR37s suppresses plant defense genes by cross-kingdom RNAi. RNA Biology 14:421–428
Wang M, Weiberg A, Lin FM, Thomma BP, Huang HD et al (2016a) Bidirectional cross-kingdom RNAi and fungal uptake of external RNAs confer plant protection. Nature Plants 2:16151
Wang S, Li P, Zhang J, Qiu D, Guo L (2016b) Generation of a high resolution map of sRNAs from Fusarium graminearum and analysis of responses to viral infection. Scientific Reports 6:26151
Weiberg A, Wang M, Bellinger M, Jin H (2014) Small RNAs: a new paradigm in plant-microbe interactions. Annual Review of Phytopathology 52:495–516
Weiberg A, Wang M, Lin FM, Zhao H, Zhang Z et al (2013) Gungal small RNAs suppress plant immunity by hijacking host RNA interference pathways. Science 342:118–123
Wu J, Wang Y, Xu J, Korban SS, Fei Z et al (2018) Diversification and independent domestication of Asian and European pears. Genome Biology 19:77
Xie J, Jiang D (2014) New insights into Mycoviruses and exploration for the biological control of crop fungal diseases. Annual Review of Phytopathology 52:45–68
Yan J, Xie Y, Zhang W, Wang Y, Liu J et al (2013) Species of Botryosphaeriaceae involved in grapevine dieback in China. Fungal Diversity. 61:221–236
Yan JY, Zhao WS, Chen Z, Xing QK, Zhang W et al (2018) Comparative genome and transcriptome analyses reveal adaptations to opportunistic infections in woody plant degrading pathogens of Botryosphaeriaceae. DNA Research 25:87–102
Yang F, Wang GP, Xu WX, Hong N (2017a) A rapid silica spin column-based method of RNA extraction from fruit trees for RT-PCR detection of viruses. Journal of Virological Methods 247:61–67
Yang T, Groenewald JZ, Cheewangkoon R, Jami F, Abdollahzadeh J et al (2017b) Families, genera, and species of Botryosphaeriales. Fungal Biology 121:322–346
Ye J, Fang L, Zheng H, Zhang Y, Chen J et al (2006) WEGO: a web tool for plotting GO annotations. Nucleic Acids Research 34:W293–W297
Ye Z, Pan Y, Zhang Y, Cui H, Jin G et al (2017) Comparative whole-genome analysis reveals articial selection effects on Ustilago esculenta genome. DNA Research 24:635–648
Yu X, Li B, Fu Y, Xie J, Cheng J et al (2013) Extracellular transmission of a DNA mycovirus and its use as a natural fungicide. Proceedings of the National Academy of Sciences of the United States of America 110:1452–1457
Zhai L, Zhang M, Lv G, Chen X, Jia N et al (2014) Biological and molecular characterization of four Botryosphaeria species isolated from pear plants showing stem wart and stem canker in China. Plant Disease 98:716–726
Zhang DX, Spiering MJ, Nuss DL (2014) Characterizing the roles of Cryphonectria parasitica RNA-dependent RNA polymerase-like genes in antiviral defense, viral recombination and transposon transcript accumulation. PLoS One 9:e108653
Zhang X, Segers GC, Sun Q, Deng F, Nuss DL (2008) Characterization of Hypovirus-derived small RNAs generated in the chestnut blight fungus by an inducible DCL-2-dependent pathway. Journal of Virology 82:v2613–v2619
Zhao Z, Cao K (2012) Epidemiological studies and prevention and control of apple ring rot. Northern Horticulture 1:127–129 in Chinese
Zheng A, Lin R, Zhang D, Qin P, Xu L et al (2013a) The evolution and pathogenic mechanisms of the rice sheath blight pathogen. Nature Communications 4:1424
Zheng W, Huang L, Huang J, Wang X, Chen X et al (2013b) High genome heterozygosity and endemic genetic recombination in the wheat stripe rust fungus. Nature Communications 4:2673
Zhou J, Fu Y, Xie J, Li B, Jiang D et al (2012) Identification of microRNA-like RNAs in a plant pathogenic fungus Sclerotinia sclerotiorum by high-throughput sequencing. Molecular Genetics and Genomics 28:275–282
We would like to thank the native English speaking scientists of Elixigen Company (Huntington Beach, California) for editing our manuscript.
This research was financially supported by grants from the National Natural Science Foundation of China (No. 31471862), the key research and development plans of the Ministry of Science and Technology [grant no. 2017YFD0201103 and 2018YFD0201406], and the earmarked fund for Pear Modern Agro-industry Technology Research System (CARS-29-15).
Availability of data and materials
This whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession SACU00000000. The version described here is version SACU00000000. Small RNA sequencing raw reads have been deposited in the NCBI database, with the Sequence Read Archive accession number SRP174595, BioProject PRJNA511629. The data generated or analyzed during this study are included in the supplementary files.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Evaluation of the pathogenicity of Botryosphaeria dothidea strains. (a) Lesions on ‘Hohsui’ pear fruit and branches, and ‘Fushi’ apple fruit induced by Botryosphaeria dothidea strains; (b) Lesion length on ‘Hohsui’ pear fruit at 9 d (b-I), ‘Fuji’ apple fruit at 4 d (b-II) and ‘Hohsui’ pear branches at 20 d (b-III) after inoculation with LW-1(LW-CP), LW-C, LW-P, Mock, HL-1 and HBWH-1 B. dothidea strains. (DOCX 653 kb)
Figure S2. The colony morphology of mycovirus-infected Botryosphaeria dothidea strains in MS culture at 25 °C darkness for 3 d (I) and developing conidiomata and conidial angle under black light with 365 nm wavelength for 5 d with the naked eye (II) and 11 d observed under stereo microscope (III). Scale bars = 0.5 mm. (DOCX 2450 kb)
Figure S3. Conidiogenous cells with developing conidia and mature conidia from B. dothidea strains cultured on PDA culture. Scale bars = 20 μm. (DOC 1051 kb)
Figure S4. The mycelial morphology of Botryosphaeria dothidea strains cultured on PDA for 3 d at 25 °C in continuous darkness observed under a light microscope. Scale bar = 10 μm. a, b, c and d represent LW-CP, LW-C, LW-P and Mock samples, respectively. (DOCX 614 kb)
Figure S5. The detection of dsRNA patterns from the offspring condia of (a) LW-CP, (b) LW-C and (c) LW-P by 1.2% agarose gel electrophoresis. (DOCX 826 kb)
Figure S6. Statistical analysis of high quality subreads with estimated genome coverage 80 X were obtained from LW-Hubei genome sequencing data. (DOCX 48 kb)
Figure S7. Whole genome synteny analysis of LW-Hubei strain by comparison with Macrophomina phascolina MS6, Neofusicoccum parvum UCRNP2, Diplodia corticola CBS 112549 and Dothistroma septosporum NZE10 at the nucleotide (nt) level. (a) The 49 scaffolds of LW-Hubei and Macrophomina phascolina MS6 were compared (Mb scale). (b) The 49 scaffolds of LW-Hubei and Neofusicoccum parvum UCRNP2 were compared (Mb scale). (c) The 49 scaffolds of LW-Hubei and Diplodia corticola CBS 112549 were compared (Mb scale). (d) The 49 scaffolds of LW-Hubei and Diplodia seriata were compared (Mb scale). (DOCX 3136 kb)
Figure S8. Venn diagram of the common and unique Core-pan genes among five Botryosphaeriaceae strains: Botryosphaeria dothidea LW-Hubei, M. phaseolina MS6, N. parvum UCRNP2, D. corticola CBS 112549 and D. seriata. (DOCX 114 kb)
Figure S9. Heatmap of the dispensable gene in five Botryosphaeriaceae strains. (DOCX 69 kb)
Figure S10. Phylogenetic tree analysis based on gene families from Botryosphaeria dothidea LW-Hubei and the eight reference fungi by the neighbor-joining method. The scale number represents branch lengths. (DOCX 65 kb)
Figure S11. Annotation of the LW-Hubei genome by Nr and KOG databases. (a) Annotation of 12,273 proteins by the Nr database. Large percentage of genes in LW-Hubei genome with homologs in Macrophomina phascolina. (b) Annotation of 2536 proteins by the KOG database for delineation of 25 Clusters of Orthologous Groups of proteins (COG). The majority of proteins (n = 292) were classified into General function prediction only. (DOCX 141 kb)
Figure S12. The precursors of five novel milRNAs and their hairpin structures in Botryosphaeria dothidea strains. The mature Bd-milRNAs are shown in yellow and milRNA*s are underlined in green. The numbers show the base locations. (PDF 203 kb)
Figure S13. First nucleotide bias in novel Bd-milRNA candidates isolated from Botryosphaeria dothidea strains in (a) Mock, (b) LW-C, (c) LW-P and (d) LW-CP libraries. (DOCX 215 kb)
Figure S14. Putative Bd-milRNAs (two novel and one known) detected by poly (A) RT-PCR. (DOCX 233 kb)
Figure S15. Expression levels of Bd-milRNA target mRNAs detected by RT-qPCR. (PDF 282 kb)
Figure S16. KEGG pathway classifications of functional enrichment for differentially expressed known (a) and novel (b) Bd-milRNA target mRNAs in mycovirus-infected Botryosphaeria dothidea strains of LW-CP, LW-C and LW-P. (PDF 125 kb)
Table S1. Sequences of qPCR primers used for analysis of the expression level of RNAi components in Botryosphaeria dothidea strains. (DOCX 14 kb)
Table S2. The sequences of primers used for analysis of the expression of putative Bd-milRNAs and target gene mRNAs in Botryosphaeria dothidea strains. (DOCX 16 kb)
Table S3. The sizes of conidia from Botryosphaeria dothidea strains. (DOCX 13 kb)
Table S4. Summary of original sequencing data for Botryosphaeria dothidea LW-Hubei isolate generated using the Illumina platform. (DOCX 13 kb)
Table S5. Summary of original sequencing data for Botryosphaeria dothidea LW-Hubei isolate generated using the PacBio platform. (DOCX 14 kb)
Table S6. Statistical analysis of basic RNA-seq and small RNA sequencing information corresponding to the LW-Hubei strain genome used in this study. (DOCX 15 kb)
Table S7. Statistical analysis of specific motifs for transmethylases. (DOCX 19 kb)
Table S8. The complete genome features of the eight fungi used as references in the comparative genome analyses. (DOCX 14 kb)
Table S9. Different methods used for prediction of functional annotation of Botryosphaeria dothidea LW-Hubei isolate. (DOCX 13 kb)
Table S10. Annotation of the CAZy categories, PHI and TF in LW-Hubei isolate. (XLSX 149 kb)
Table S11. Predicted effectors (n = 12) from Botryosphaeria dothidea LW-Hubei can be annotated by PHI. (DOCX 15 kb)
Table S12. Annotated transporters (n = 552) of Botryosphaeria dothidea LW-Hubei isolate. (DOCX 13 kb)
Table S13. Predicted kinases (n = 128) from Botryosphaeria dothidea LW-Hubei isolate. (DOCX 14 kb)
Table S14. Summary of accession numbers for RNAi components from the fungal strains analyzed in this study. (DOCX 15 kb)
Table S15. Summary of common and specific sRNA sequences and mean frequencies in the (a) LW-CP and Mock libraries, (b) LW-P and Mock libraries, and (c) LW-C and Mock libraries constructed from Botryosphaeria dothidea strains. (DOCX 15 kb)
Table S16. The 62 candidate novel Bd-milRNAs and six novel Bd-milRNAs among the milRNA*s detected in Botryosphaeria dothidea strains. (DOCX 32 kb)
Table S17. The differentially expressed known Bd-milRNAs in the LW-CP/Mock, LW-C/Mock and LW-P/Mock libraries constructed from Botryosphaeria dothidea strains by sRNA sequencing. (DOCX 29 kb)
Table S18. Three known Bd-milRNAs sequences were cloned and sequenced by stem-loop RT-PCR (DOCX 50 kb)
Table S19. The target gene of known Bd-milRNAs predicted by psRobot and/or TargetFinder software. (XLSX 575 kb)
Table S20. The target gene of novel Bd-milRNAs predicted by psRobot and/or TargetFinder software. (XLSX 20 kb)
About this article
- Pear ring rot
- Botryosphaeria dothidea
- Botryosphaeria dothidea chrysovirus 1 (BdCV1)
- Botryosphaeria dothidea partitivirus 1 (BdPV1)