Characterization of the Rosellinia necatrix Transcriptome and Genes Related to Pathogenesis by Single-Molecule mRNA Sequencing.

White root rot disease, caused by the pathogen Rosellinia necatrix, is one of the world's most devastating plant fungal diseases and affects several commercially important species of fruit trees and crops. Recent global outbreaks of R. necatrix and advances in molecular techniques have both increased interest in this pathogen. However, the lack of information regarding the genomic structure and transcriptome of R. necatrix has been a barrier to the progress of functional genomic research and the control of this harmful pathogen. Here, we identified 10,616 novel full-length transcripts from the filamentous hyphal tissue of R. necatrix (KACC 40445 strain) using PacBio single-molecule sequencing technology. After annotation of the unigene sets, we selected 14 cell cycle-related genes, which are likely either positively or negatively involved in hyphal growth by cell cycle control. The expression of the selected genes was further compared between two strains that displayed different growth rates on nutritional media. Furthermore, we predicted pathogen-related effector genes and cell wall-degrading enzymes from the annotated gene sets. These results provide the most comprehensive transcriptomal resources for R. necatrix, and could facilitate functional genomics and further analyses of this important phytopathogen.

Rosellinia necatrix, one of the most prevalent soil-borne fungal phytopathogens, infects both the surfaces of roots and the tissue under the bark in a broad range of commercially important crop plants such as cotton, nuts, apples, cherries, and pears (Eguchi et al., 2009). Many of these species are of considerable economic importance, resulting in the listing of 343 fungus-host combinations by the United States Department of Agriculture (http://nt.arsgrin.gov/fungaldatabases/index.cfm). Owing to the physically devastating symptoms of R. necatrix-mediated disease in fruit trees, it has been recently considered to be of the most significant agricultural problems in the world (ten Hoopen and Krauss, 2006). The white root rot disease caused by R. necatrix results in varying degrees of canopy decline, followed by leaf drop, wilting, and eventually the death of infected plants (Kondo et al., 2013). Recently, a detailed process of the infection of avocado roots by R. necatrix has been reported using a GFP-harboring strain (Pliego et al., 2009). Upon infection, the mycelia of the fungus dispersed throughout the soil initially aggregated over the root surface. Penetration occurs by either directly forming a sclerotium or via natural opening sites such as lenticels and wounds (Pliego et al., 2012). Although the means of physical penetration by R. necatrix have been well characterized, the genetic mechanisms underlying virulence and fungal growth factors are largely unknown.

The mechanisms of general virulence factors and infection strategies used by plant pathogens to gain entry into the interior of their hosts have been well studied (Macheleidt et al., 2016). To overcome physical barriers, such as rigid plant cell walls and waxy cuticle layers covering the epidermis, the phytopathogens—particularly filamentous fungi—have evolved a diverse array of mechanisms. An example of such strategies are specialized infection structures, which secrete plant cell wall degrading enzymes (CWDEs) and various plant immune suppressors (i.e., effectors). It has been well reported that the various enzymes secreted by fungal pathogens to depolymerize the components of the plant cell wall include cellulose, hemicellulose, and pectin (Kubicek et al., 2014). In addition, many types of pathogenic effectors have been well characterized and their biological roles in disabling the defense/immune systems of host plants have been extensively studied (Dou and Zhou, 2012). A combination of advanced genome sequencing technologies and systems biology has begun to unravel the functional diversity and genetic structures of pathogenic effectors and CWDEs (Chatnaparat et al., 2016). However, the lack of more basic genetic information regarding R. necatrix, most notably that of the transcriptome, is currently limiting our understanding of the mechanisms underlying the pathogenhost interactions of this harmful fungal plant pathogen and economically important crop species.

Recent advances in next-generation sequencing (NGS) technologies have led to the discovery of an enormous variety of genomic information. The development of high-throughput RNA-Seq technology has enabled the identification of de novo assembled unigene sets and the quantification of gene expression levels (Goodwin et al., 2016). However, the reduced sensitivity of short-read assemblies from RNA-Seq data could be a barrier to improving the accuracy of gene model predictions in the absence of reference genome data. The Pacific Biosciences (PacBio) single-molecule sequencing platform could alternatively provide direct and reliable sequence information as full-length transcripts with a mean read length of > 7 kb (Tilgner et al., 2014). Thus, in the present study, we used PacBio long-read mRNA sequencing technology to obtain reliable and comprehensive transcriptome sequences from R. necatrix. As a result, here we present 10,616 full-length transcriptomes and gene sets related to cell cycle control and potential pathogenicity. These results of the present study could improve genome annotation as well as our understanding of the pathogenicity of this harmful plant pathogen.

Materials and Methods

Preparation of fungi samples and RNA isolation

Two strains (KACC 40445 and 40168) of R. necatrix were kindly distributed by the Korean Agricultural Culture Collection (KACC) (http://genebank.rda.go.kr/microbeMain. do). To test the growth rate of R. necatrix, two strains were incubated in the dark for indicated times at 25°C on potato dextrose agar (PDA). To extract total RNA, strain KACC 40445 was cultured in potato dextrose broth (PDB) liquid medium for 5 days at 25°C. Total RNA was extracted using a RNA extraction kit (iNtRon Biotech, Seoul, Korea).

Library construction and PacBio single-molecule long-read sequencing

High-quality RNA was used for cDNA synthesis and library construction. cDNA was synthesized using the Clontech SMARTer PCR cDNA Synthesis Kit (Clontech Laboratories, Mountainview, CA, USA). Cycle optimization was performed to determine the optimal number of cycles for large-scale PCR. After performing large-scale PCR, 3 cDNA size fractions (1–2 kb, 2–3 kb, and 3–6 kb) were prepared using the BluePippin size selection system for NGS (Sage Science, Beverly, MA, USA). Each sample was used as input for library preparation. The SMRTbell library was constructed by using the SMRTbell™ Template Prep Kit 1.0 (PN 100-259-100; Pacific Biosciences, Menlo Park, CA, USA). After a sequencing primer is annealed to the SMRTbell template, DNA polymerase is bound to the complex (DNA/Polymerase Binding Kit P6; Pacific Biosciences). Following the polymerase-binding reaction, the MagBead Kit is used to bind the library complex with MagBeads prior to sequencing. MagBead-bound complexes provide for more reads per SMRT cell. This polymerase-SMRTbell-adaptor complex is then loaded into zero-mode waveguides (ZMWs). The SMRTbell library was sequenced using a total of 8 SMRT cells (1–2 kb, 2 cells; 2–3 kb, 2 cells; 3–6 kb, 4 cells; Pacific Biosciences) using C4 chemistry (DNA sequencing Reagent 4.0), and 1 × 240 min movies were captured for each SMRT cell using the PacBio RS sequencing platform (Pacific Biosciences). Long reads were identified using the SMRT Analysis v2.2 RS_IsoSeq.1 classification protocol. Subsequently, reads were clustered into full-length reads derived from the same isoform and polished consensus sequences using TOFU pipeline (Gordon et al., 2015).

Gene predictions and annotations

Automated gene prediction was performed using the automated annotation pipeline MAKER (Cantarel et al., 2008; Holt and Yandell, 2011). Gene annotations were made using all sequences from the fungal order Xylariomycetidae from the Uniprot/Swiss-Prot protein database. Ab initio gene predictions were created by MAKER using the programs SNAP (Korf, 2004) and Augustus (Stanke et al., 2006). Gene models were further improved by providing MAKER with PacBio RNAseq data. Three iterative runs of MAKER were used to produce the final gene set (Supplementary Fig. 1; Cantarel et al., 2008). Post annotation was performed to add putative gene functions protein domains by BLASTing against National Center for Biotechnology Information (NCBI) NR, UniProt/Swiss-Prot, and InteProScan databases (Jones et al., 2014). Based on the annotated gene, gene ontology (GO) and KEGG analysis were used with the default setting “Blast2GO” to confirm overall function (Conesa and Götz, 2008).

Real-time RT-PCR

Primers were designed using Primer 3.0 (http://primer3plus.com/). PCR reactions were carried out using the Applied Biosystems™ SYBR Green Master Mix on an Applied Biosystems Quant Studio3 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). PCR conditions were as follows: 95°C for 30 s; followed by 40 cycles of 95°C for 5 s, 58°C for 30 s, and a dissociation stage of 95°C for 15 s; 60°C for 60 s, and 95°C for 15 s. The expression of selected genes was normalized against GAPDH of R. necatrix, which was used as an internal reference gene. The relative gene expression was calculated using the 2−ΔΔCt method. These experiments were repeated using two biological replicates. Primer sequences are listed in Supplementary Table 1.

Florescence microscopy

For florescence microscopy, mycelial samples from R. necatrix were stained with Congo Red (0.5% w/v in water) for 15 min, and then briefly washed in distilled water. Specimens were observed using the Texas red channel of a Nikon florescence microscope (Eclipse Ti; Nikon, Tokyo, Japan) equipped with a 40×.

Accession codes

The PacBio sequence data generated for this work are accessible via the NCBI Sequence Read Archive under accession number (SRA: SRP096541 and BioProject: PRJNA358652).

Results and Discussion

In order to obtain as many expressed transcript sequences as possible from R. necatrix mycelial tissues, we extracted high-quality total RNA from R. necatrix (KACC 40445 strain) (Supplementary Fig. 2A). Then, we prepared three size-fractionated cDNA libraries (1–2 kb, 2–3 kb, and 3–6 kb) and subjected them to PacBio single-molecule, long-read sequencing (Supplementary Fig. 2B). A total of 605,153 reads were obtained using a total of 8 SMRT Cells (average size of 2,584 bp), and the sequences were processed using ToFu (Fig. 1A). After filtering low-quality and chimeric regions, 109,251 high-quality isoforms were selected. The high-quality isoforms contained 73,819 redundant isoforms, thereby 35,432 isoforms were finally subjected to further analyses for gene annotation and predicting gene functions. As a result, we successfully filtered out 28,953 repeated sequences, and annotated 10,616 uniquely expressed transcripts of R. necatrix mycelial tissues (Supplementary Table 2).

Fig. 1

PacBio Iso-Seq library construction and comparison of full-length transcripts and isoforms between the Rosellinia necatrix draft genome and PacBio data. (A) Summary of sequence from R. necatrix PacBio SMRT cells. FL, full length. (B) Comparison of full-length transcripts and isoforms between R. necatrix draft genome and PacBio data. (C) Comparison of isoforms produced by PacBio and R. necatrix draft genome.

Next, we compared the 10,616 unique isoforms to the R. necatrix draft genome (NCBI GCA_001445595.1), and classified the isoforms into eight groups (Fig. 1B). Of the isoforms, approximately 5,455 (43%) transcripts were predicted to be novel isoforms where at least one splicing junction was shared with a reference transcript. Other isoforms were categorized as: isoforms completely matching a reference intron (3,907; 31%), isoforms with exonic overlap with a reference transcript (1,591; 13%), isoforms overlapping both exons and introns of a reference transcript (506; 4%), isoforms contained in a reference gene (419; 3%), isoforms exonic overlap with reference on opposite strand (393; 3%), intergenic isoforms (318; 3%), polymerase run-on isoforms (26; 0.2%), and isoforms overlapping a reference intron on the opposite strand (16; 0.1%). Using PacBio Iso-Seq data, we identified a large quantity of novel gene sets with two more isoforms. Subsequently, to assess the quality and novelty of our sequencing data, we compared the transcript lengths between PacBio data and the R. necatrix draft genome. As shown in Fig. 1C, our selected sequenced transcripts recovered transcripts longer than those analyzed by draft genome sequence of R. necatrix deposited in NCBI. Total 10,255 predicted nucleotide sequences of R. necatrix are available in NCBI GenBank, but only 4,674 transcript sequences were completely overlapped with 10,616 PacBio sequences. In addition, our analyzed data displayed that 91.3% of transcripts (9,692) were matched over 50% identity with analyzed transcript sequences. These results indicate that our newly obtained transcriptome data of R. necatrix by a PacBio single-molecule long-read sequencing platform could be a more reliable resource than previously reported data for genome information of this important plant pathogenic fungi.

GO could define and classify the gene functions of annotated sets (Ashburner et al., 2000). To predict the functions of the annotated R. necatrix genes, the GO database and Blast2GO were used. In GO annotation, 7,012 (66.05%) genes were matched to the database. Of the 31 resulting functional groups, the majority of sequences where grouped with genes involved in metabolic processes and cellular processes. In the metabolic pathway, we compared sequences to genes involved in nitrogen-related metabolism and carbohydrate metabolism, which enhance host-pathogen interactions (Fig. 2). Cellular processes represent the work that takes place on a cellular level, and thus we identified genes in the categories cell cycle and single-organism cellular process (Fig. 2).

Fig. 2

Frequency distribution of gene ontological (GO) classification in Rosellinia necatrix. The results are categorized into biological processes, molecular functions, and cellular components using Blast2GO (https://www.blast2go.com/). The upper y-axis indicates the number of genes in category.

A KEGG analysis facilitates the understanding of a high-level of functional study performed using high-throughput technologies (Kanehisa et al., 2016). It can also allow for the identification of enzymatic function via BLAST. A total of 1,613 annotated genes belonged to the pathway in 126 KEGG databases (Supplementary Table 2). The majority of genes in this analysis were assigned to the metabolism-related pathway, similar to the results of GO annotation. The largest quantity of matched genes were identified as hydrolases, some of which are known to cause specific interactions in host plant-pathogen systems (Yarullina et al., 2016).

The cell cycle has been shown to be one of the most important factors in the establishment of plant disease via appressorium morphogenesis and formation in the rice blast fungus, Magnaporthe oryzae (Saunders et al., 2010).

Interestingly, in the present study, we observed different growth rates of the two distributed strains, where the mycelia of KACC 40445 strain grew more vigorously in PDA media in comparison to the 40168 strain (Fig. 3A). We also confirmed that the 40445 R. necatrix strain grew approximately 2–3 times faster (Fig. 3B), indicating a role for the activation of cell cycle control in this strain. To test whether the cell cycle was activated in the 40445 strain, we first compared hyphae size between the two strains by staining cell walls with Congo red. The morphology and size of hyphae between the two strains were similar, suggesting that different mechanisms govern the regulation of cell cycle and division in these two strains (Fig. 3C). Next, we selected putative genes related to cell cycle checkpoints and cell division in R. necatrix (Supplementary Table 3, 4, Supplementary File 1), and assessed their expression levels in actively dividing mycelial tissues. As expected, two representative positive regulators of cell cycle control (CDC123 and PAS1) were upregulated, and most negative regulators (i.e., CDC5, CDC14, WHI5, CDT2, and RES2) were downregulated in the 40445 strain (Fig. 3D). These results indicate that the control of the cell division of hyphae is one of the most critical factors for mycelia growth in R. necatrix. Recent research has shown that independent S-phase cell cycle checkpoints could be critical factors in the infection of host plants by M. oryzae (Osés-Ruiz et al., 2017). Our results indicate that the increased expression of the G1/S specific cyclin PAS1 in the KACC 40445 strain might not only enhance the activity of cell cycles, but could also affect its pathogenicity. More detailed studies of pathogenesis in different strains of R. necatrix should be carried out to further increase our understanding of pathogenicity in this species.

Fig. 3

Identification and functional analysis of genes related to cell cycle and division in the Rosellinia necatrix genome. (A) Growth rate of mycelial tissue of two R. necatrix strains on potato dextrose agar (PDA). Five-day-old mycelial tissues of KACC 40445 and 40168 are displayed. (B) The length of mycelial tissues of the two different strains were measured and plotted for five days. (C) The cell phenotype of mycelium of the two different R. necatrix strains. Congo Red-stained hyphae of KACC 40445 and 40168 were observed using florescence microscopy. Yellow arrows indicate septa of hyphae. Scale bars = 50 μm. (D) The expression levels of genes related to cell cycle or division were determined using real time RT-PCR (n = 3 biological replicates). Error bars indicate standard errors.

Fungal phytopathogens are known to be responsible for huge losses of economically important crop plants. The secretion systems of these fungi deliver special secretory proteins called effectors that trigger injection and translocation, along with signal sequences, and secrete specific effectors into the host via different systems (Tseng et al., 2009). Thus, while plants have evolved an innate immune system that recognizes specific patterns of pathogens and rapidly induces diverse immune responses, phytopathogens have developed effectors to avoid the innate immunity of host plants or inactivate toxic metabolites (Lo Presti et al., 2015). Recently, the development of high-throughput technologies has yielded the genomes and transcriptomes of a variety of fungal species. These analyses have revealed several critical proteins involved in pathogen-host interactions, particularly CWDEs and effector proteins (Sonah et al., 2016). Thus, we next aimed to define gene sets related to fungal pathogenicity, including effectors and enzymes involved in the degradation of plant cell walls. To predict effector genes, we used previously published secretion system effector genes (Costa et al., 2015) and the PHI (Pathogen Host Interactions) database (http://www.phi-base.org/), where PHI-base is a database of virulence and effector genes that have been experimentally proven via pathogen-host interactions (Urban et al., 2017). We used BLASTp to match PHI-base with an e-value cutoff of 1E-10 and 30% identity to more accurately predict effector genes. As a result, 54 genes found in the PacBio gene sets were annotated to 26 effector genes. Each of the matched gene functions and secretion system types are presented in Supplementary Table 5 and Supplemenatry File 2. The presence of 26 types, with at least 54 gene sets encoding pathogenic effector genes, is expected to have a wide range of effects on the host plant through a variety of functions, including transporters for effector secretion, the direct release of effector proteins into the target cells, defects in organelle trafficking, etc. Additionally, we used an EffectorP protocol (Sperschneider et al., 2016) and predicted putative 1,672 effector genes (Supplementary File 2). We also predicted 74 putative novel CWDEs among the 10,616 unique annotated genes by BLAST search against to NCBI NR and UniProt database (Supplementary Table 1 and Supplementary File 3). The CWDEs in the transcriptome of R. necatrix were commonly categorized as various celluloses, hemicelluloses, and pectinases. Owing to the critical roles of CWDEs in overriding the barrier of the host plant cell wall during infection, the numerous kinds of CWDEs present in R. necatrix likely result in compatibility with a broad spectrum of plant hosts for pathogenesis. In addition, during the late stages of pathogenic invasion, CWDEs are important because of requirement of all phytopathogenic fungi (Quoc and Chau, 2016).

In the present study, we analyzed the full-length transcriptome of R. necatrix using SMRT sequencing techniques and compared it to draft genome. In comparison to de novo-assembled draft genome information, more than 70% of the isoforms identified in this study are novel, and over 5,000 genes have been identified (Fig. 1). These results indicate that the PacBio Iso-Seq platform could accurately provide more reliable genomic resources from organisms for which genomic information is still incomplete. In addition to providing a more complete genomic picture, Iso-Seq techniques could provide precise information regarding transcription start sites and alternatively spliced exons. To define the functions of the resulting genes and their importance at the genomic level, unidentified genes are currently under investigation using a combination of SMRT sequencing techniques and biochemical/phenotypical analyses. Further work in the functional genomics of pathogen-related genes and mechanisms governing the control of growth in R. necatrix will enable the elucidation of novel and specific targets to protect crop plants from this devastating phytopathogen.