Beauveria bassiana ERL836 and JEF-007 with similar virulence show different gene expression when infecting western flower thrips, Frankniella occidentalis

doi:10.21203/rs.3.rs-17095/v2

Download PDF

Research article

Beauveria bassiana ERL836 and JEF-007 with similar virulence show different gene expression when infecting western flower thrips, Frankniella occidentalis

https://doi.org/10.21203/rs.3.rs-17095/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 27 Nov, 2020

Read the published version in BMC Genomics →

You are reading this older preprint version

Read the latest preprint version →

Background: Insect-killing fungal species, Beauveria bassiana, is as an environment-friendly pest management tool, and many isolates are on the track of industrialization. However, some of B. bassiana isolates show similar morphology and virulence against insect pests, and so it is hard to differentiate them. Herein we used two patented isolates, ERL836 and JEF-007, and investigated their virulence against western flower thrips, Frankliniella occidentalis, and further analyzed genome structures and transcriptional responses when infecting the thrips to see possible differences.

Results: The two isolates showed no significant differences in fungal growth, conidial production, and virulence against thrips, and they were structurally similar in genome. But, in transcription level, ERL836 appeared to infect thrips easily, while JEF-007 appeared to have more difficulty. In the GO analysis of ERL836 DEGs (differentially expressed genes), the number of up-regulated genes was much larger than that of down-regulated genes, when compared to JEF-007 DEGs (more genes down-regulated). Interestingly, in the enrichment analysis using shared DEGs between two infecting isolates, plasma membrane-mediated transporter activity and fatty acid degradation pathway including cytochrome P450 were more active in infecting ERL836.

Conclusion: The two B. bassiana isolates had similar morphology and virulence as well as genome structure, but in transcription level they differently infected western flower thrips. This comparative approach using shared DEG analysis could be easily applied to characterize the difference of the two B. bassiana isolates, JEF-007 and ERL836.

Epigenetics & Genomics

Beauveria bassiana

western flower thrips

genome

transcription

cytochrome P450

Beauveria bassiana has received great interest from both academia and industry because of its great potential as an alternative environment-friendly pest management tool to substitute synthetic pesticides that can damage the environment and cause insect resistance. Popularly available B. bassiana products are BotaniGard^® (isolate name: GHA), Naturalis-L^® (ATCC74040), Chongchaesak^® (ERL836), Broadband^® (PPRI 5339), and BioCeres^® (ANT-03) and they are used to control agricultural insect pests such as moths, beetles, whitefly, aphid, mite, and thrips among others [1]. Particularly B. bassiana ERL836 GR has been developed by LG-Chemical-affiliated FarmHannong and successfully launched and sold out in a local market. Insect resistance to pesticides may not be an issue when using these fungal pathogens for pest control because this fungal group uses complicated mechanical hyphal penetration and enzymatic degradation to kill the host. This mode of fungal action is completely different from that of chemicals targeting synapses or energy metabolism. Although other entomopathogenic fungi such as Metarhizium, Cordyceps (previously Isaria or Paecilomyces), Akanthomyces (previously Lecanicillium and Verticillium) are also being investigated as pest management tools, B. bassiana has several advantages over these other fungi in terms of field application, mass production, and long-term storage.

Numerous B. bassiana isolates have been collected form natural environments and have been transferred to fungal libraries around the world. USDA-ARS manages a huge entomopathogenic fungal library and support academic research worldwide. Japan and China have big fungal libraries and Brazil and European countries also have entomopathogenic fungal libraries. Industrially important isolates from these libraries are placed in the research and development track for development as commercial or government-supported products. Very few isolates successfully pass all steps of industrialization to be launched in the biopesticide market. Nevertheless, the number of B. bassiana isolates being commercially developed is increasing. An important intellectual property (IP) issue is how to scientifically determine whether a particular B. bassiana isolate is different from another isolate. Attempts have been made to use phylogenetic analysis based on one or a couple of housekeeping genes, but this approach is not always adequate. As a prerequisite for intellectual property in South Korea, such as a patent for a particular isolate, clear experimental data that describes the uniqueness of an isolate needs to be submitted; this type of data could be used to resolve possible isolate-mediated conflicts in the biopesticide industry.

As introduced above, phylogenetic approaches have been mainly used to determine the evolutionary origins of isolates and their genetic relationships. Rehner and Buckley investigated the correspondence between Beauveria and Cordyceps using EF1-α and ITS phylogenies [2]. Glare et al. analyzed rDNA to determine phylogenetic relationships among 26 Beauveria isolates; more specifically, they used sequence information from the 3' end of 16s rDNA across ITS 1, 5.8s rDNA, and ITS 2 to the 5' end of 28s rDNA [3]. In phylogenies, three nuclear genes encoding elongation factor 1-α (TEF1), RNA polymerase II largest subunit (RPB1), and RNA polymerase II second largest subunit (RPB2) were used by Rehner and colleagues [4]. Furthermore, a new species, B. lii, was found by four-locus based phylogenetic analysis [5]. However, these previous phylogenetic analyses were based on one or a couple of house-keeping genes, and could not be easily used to suggest that one isolate is different from another isolate in B. bassiana under IP issues.

Although fungal phylogenetic analysis can be used to suggest the uniqueness of one isolate in patent submission and acquisition, but the question is how many genes should be analyzed for this purpose. Whole genome sequencing (WGS) can provide the sequences of the entire genome of an organism and to compare the roles of genes or the diversity of genes in the isolates of same species of fungi. WGS of 10 B. bassiana isolates have been performed to date. ARSEF2860 (33.70 Mb, ADAH00000000.1) was the first WGS of B. bassiana and was obtained using Roche 454 system and Illumina paired-end sequencing in 2012 [6]. Genomes of other B. bassiana isolates, such as D1-5 (36.69 Mb, ANFO00000000) in 2014, ARSEF1520 (36.97 Mb, JTCW01000000), ARSEF 2597 (38.83 Mb, JTCX00000000), ARSEF8028 (35.02 Mb, JRHA00000000), ARSEF5078 (34.45 Mb, JTCZ00000000), ARSEF4305 (34.77 Mb, JTCY00000000) in 2016, BCC2660 (34.56 Mb, MWYT00000000) in 2017, and JEF-007 (36.54 Mb, MRVG00000000 from our previous work) and Bv062 (34.84 Mb, MAQY00000000) in 2018 have been generated. Whole genome sequences of D1-05, ARSEF1520, ARSEF 2597, ARSEF5078, and ARSEF4305 were obtained by Illumina sequencing.

To more deeply understand fungal mode of action, transcriptome analyses have been conducted, and this approach could be used to detect differences in patterns of gene expression among isolates of B. bassiana. De novo sequencing of non-model organisms is the most widely used strategy for transcriptomic profiling [7, 8, 9]. Some studies have investigated the interaction between B. bassiana and host insects at the molecular level. Transcriptome analysis of the initial phases of B. bassiana infection of the coffee berry borer, Hypothenemus hampei, was conducted [10]. Gene expression in silverleaf whitefly, Bemisia tabaci, in response to the infection by B. bassiana has also been analyzed [11]. Immunity-related gene expression in the Asian corn borer, Ostrinia furnacalis, when infected by B. bassiana was studied using an RNA-seq approach [12, 13]. Infection of B. bassiana to bean bug, Riptortus pedestris has also been investigated [14], as has the resistance and susceptibility of two silkworm species to B. bassiana infection [15]. Many virulence-related genes were found to be up-regulated in B. bassiana when infecting Anopheles stephensi mosquito [16]. Chen et al. studied the molecular mechanisms of B. bassiana infection to wax moth, Galleria mellonella [17]. Recently, in our laboratory Lee et al. performed transcriptome analysis of the bean bug in response to infection by B. bassiana [18]. In 2017, Wang et al. conducted a study to evaluate transcriptomic differences in two B. bassiana isolates when infecting silkworms [19].

Of the insects to be managed, thrips receives many attentions because of its high resistance to chemicals and cryptic behaviors. Transcriptome analyses of thrips infected with various entomopathogens have been conducted, but most studies have focused on gene expression of thrips, rather than fungal genes. Differentially expressed genes in western flower thrips, Frankliniella occidentalis, in response to tomato spotted wilt virus infection were analyzed by RNA-seq [20-23]. Zhang and colleagues identified a total of 36,339 thrips unigenes, and among them, 278 genes were involved in insecticide resistance [20]. Similarly, molecular responses of melon thrips, Thrips palmi, to capsicum chlorosis virus infection were analyzed by RNA-seq [24].

In this study, two patented B. bassiana isolates, JEF-007 (Patent No: 10-1666968, South Korea) and ERL836 (Patent No:10-1974265, South Korea; commercialized in a local market) were used and the following biological features were compared: hyphal growth, conidial production and virulence against thrips. Secondly, the genomes of the two isolates were compared. The whole genome of ERL836 (GenBank accession: PPTI00000000) was sequenced using PacBio technology and the genome sequence data for JEF-007 (GenBank accession: MRVG00000000) was obtained from our previous work [1]. Lastly the gene expression patterns of ERL836 and JEF-007 when infecting western flower thrips, F. occidentalis, were analyzed using an RNA-seq technology.

Fungal morphological growth

The two isolates showed similar hyphal growth and conidial production (Figure 1). The two isolates grew similarly on the 1/4SDA medium with white hyphal mass and conidia (Figure 1A). On the millet-based solid culture, JEF-007 and ERL836 produced 3.98 ± 0.71 × 10⁹ conidia g^-1 and 5.12 ± 1.06 × 10⁹ conidia g^-1, respectively and no significant difference of conidial production was observed (F_1,12 = 1.7, p = 0.783) (Figure 1B).

Virulence against thrips

In the bioassay against western flower thrips, JEF-007 and ERL836 showed similarly high virulence against the adults under laboratory conditions (Figure 2). When the adults were sprayed with the two isolates, the treated thrips showed fungus dose-dependent response and no significant mortality between the two isolate treatments. In 3 days after treatment, LC₅₀ of JEF-007 was 7.9×10⁶ (1.8×10⁵ ~ 3.4×10⁸) conidia ml^-1 (R² = 0.881) and LC₅₀ of ERL836 was 2.2×10⁷(3.6×10⁶ ~ 1.4×10⁸) conidia ml^-1 (R² = 0.916) (Figure 2A). The lethal time 50 (LT₅₀) of the two isolates were similar in each conidial concentration. At 1×10⁷conidia ml^-1, the LT₅₀ of JEF-007 was 3.04 (2.41~3.91) days and that of ERL836 was 2.87 (1.92~4.30). Similar hyphal growth and sporulation on the insect cadavers (mycosis) were observed in the two fungal isolates-treated adults in 5 days after the treatment (Figure 2B).

Genome features of two B. bassiana isolates

From the de novo sequencing of the two isolates (Supplementary Table S1), the genome sizes of B. bassiana JEF-007 and ERL836 were similar to each other; 36.5 Mb of JEF-007 and 35.5 Mb of ERL836 (Table 1). These genome sizes were similar to the previously analyzed B. bassiana ARSEF2860 (33.7 Mb). Most genomic features of the three B. bassiana isolates were not significantly different, except for different scaffold numbers due to different sequencing platforms (ARSEF2860: 242, ERL836 and JEF-007: 15 and 39, respectively). The G+C content of the three isolates was approximately 48~52%, and the numbers of protein-coding genes were 10,631 (ERL836), 10,857 (JEF-007) and 10,366 (ARSEF2860).

Comparison of genome structure of two B. bassiana isolates

No big significant differences in genome structure were observed between the two B. bassiana genomes (Figure 3A). When ERL836 genome was aligned with JEF-007 genome, although there were some minor differences, but mostly genome structures (synteny) were similar to each other. When the genome of ERL836 was aligned to that of ARSEF2860 as a reference genome, two genome structures were similar, although some minor parts of the middle of the ARSEF2860 genome did not match the ERL836 genome. Additionally, when the three B. bassiana isolates were subjected to KEEG pathway analysis, no significant differences were detected (Figure 3B). Metabolic pathway (9.7%), biosynthesis of secondary metabolites (3.9%) and biosynthesis of antibiotics (2.8%) were commonly major pathways in the three isolates. When orthologs of the two B. bassiana isolates were analyzed, ERL836 shared 8,847 genes with JEF-007 and the two isolates shared 8,409 genes with ARSEF2860 as a reference (Figure 3C). In the three B. bassiana isolates, the percentages of the unique gene of each isolate were 0.19~0.35% of the total shared genes and the numbers are quite small.

Comparative DEGs of two infecting B. bassiana isolates

B. bassiana ERL836 or JEF-007 was cultured on the 1/4SDA medium and adults of western flower thrips were exposed to the fungal mass for three days. Collected infecting fungus was subjected to RNA extraction and differentially expressed genes were compared. In the DEG analysis using the JEF-007 sequencing reads (Supplementary Table S2), a total of 643 genes were up-regulated, and 982 genes were down-regulated when infecting western flower thrips (Figure 4A). Approximately 1.5 times more contigs belonged to down-regulated genes than up-regulated genes. Most fold changes were in the range of -2 to 2 and accounted for 80.9% of total DEGs (1,625 DEGs). Among up-regulated genes (643 genes), 51% were assigned to the biological process category, 32% to the molecular function category, and 17% to the cellular component category. Among down-regulated genes (982 genes), 44% were assigned to the molecular function category, 41% to the biological process category, and 15% to the cellular component category. Both up- and down-regulated genes in the biological process category were involved in metabolic processes. Most up-regulated genes in the molecular function category had catalytic activity, whereas most down-regulated genes in this category were involved in binding. A similar number of contigs in most GO categories of the cellular component category were found among up-regulated genes, while genes with a membrane function were some of the most down-regulated genes.

In the DEG analysis using the ERL836 sequencing reads (Supplementary Table S2), a total of 1,197 genes were up-regulated, and 360 genes were down-regulated (Figure 4B). ERL836 treatment resulted in approximately twice more up-regulated contigs than down-regulated contigs. Most fold changes were in the range of–2 to 2 and accounted for 71.96% of total DEGs (2,471 DEGs). Among up-regulated genes (1,197 genes), 38.4% were in the biological process category, 43.1% in the molecular function category, and 18.5% in the cellular component category. Similarly, among down-regulated genes, 39.8% were in the biological process category, 40.0% in the molecular function category, and 20.1% in the cellular component category. Metabolic process-related genes were the most abundant up-regulated genes. Among genes with a molecular function, most contigs were annotated as encoding proteins with catalytic activity followed by binding and heterocyclic compound binding. Among cellular component genes, similar numbers of contigs were found in most of the GO categories.

Randomly selected up-regulated genes from the ERL836 and JEF-007 contigs were subjected to qRT-PCR for validation. The fold change levels of the selected genes as analyzed by qRT-PCR were consistent with the Illumina sequencing analysis results (Supplementary Figure S1).

GO enrichment of shared DEGs between two isolates

A list of shared DEGs between infecting B. bassiana ERL836 and JEF-007 was derived from RNA-seq data for enrichment analysis. A total 6,355 DEGs were shared by the two infecting B. bassiana isolates, but 3,576 genes without gene IDs were excluded from the analysis (Figure 5A). To investigate differences, 81 genes were trimmed based on the fold change value between infecting ERL836 and JEF-007 (|FC(ERL836-JEF-007)|>3). From the GO enrichment analysis (Supplementary Table S3), 30 GO terms were significantly affected by shared DEG groups (FDR<0.05) (Figure 5B). Transporter activity (GO:0005215), plasma membrane part (GO: 0044459), and intrinsic component of plasma membrane (GO:0031226) functions were profoundly affected. According to pathway (KEGG) from the GO enrichment, the significantly affected pathway was the fatty acid degradation pathway (KEGG:00071) (-log(p)>1.4), which is related to cytochrome P450 (Figure 5C). Some cytochrome P450 genes (CYP gene) in the shared gene groups showed different gene expression levels; highly up-regulated in ERL836 but down-regulated in JEF-007. The fold change value of CYP539B1 in ERL836 was 2.34, while in JEF-007 it was -0.73. The CYP655C1 had fold change value of 3.87 in ERL836 but -0.1 in JEF-007. Additionally, CYP5099A1 was upregulated at 2.66-fold in ERL836, but downregulated at -3.21-fold in JEF-007.

In this work, the two B. bassiana isolates, ERL836 and JEF-007 showed high similarities in biological features, such as morphology, conidial productivity and virulence against western flower thrips. Not only these two isolates, but also many other B. bassiana isolates look similar and are regarded to be similar active ingredient when being developed as biopesticides. The highly similar morphological characteristics of isolates can lead to potential conflicts when they are submitted to patent application and product development. So, in this work to further characterize two patented B. bassiana isolate ERL836 and JEF-007, their genome structures and gene expression patterns when infecting western flower thrips were analyzed and compared. The genome structures of the two isolates were very similar, but gene expression patterns were quite different. Different transcriptional responses could therefore be used as a characteristic to differentiate the two B. bassiana isolates.

The two B. bassiana isolates, JEF-007 & ERL836 had similar genome structures, but their scaffold numbers were different from B. bassiana ARSEF2860, which might be due to the different sequencing technologies used to obtain the genomes of the isolates. JEF-007 and ERL836 isolates were sequenced using the recent PacBio sequencing technology which reads ca. 10 Kb at a time (maximum 20 Kb) followed by error correction, whereas the genome of ARSEF2860 was generated using Illumina pair-end sequencing technology, which reads short fragments of 100 bp. Ten and nine contigs of more than 1 Mb (considered the minimum size of chromosome) were generated for JEF-007 and ERL836, respectively. It is therefore possible that the two isolates have ca. 10 chromosomes, whereas ARSEF2860 has been reported to have six chromosomes [25]. The ERL836, JEF-007 and ARSEF2860 as a reference shared 8,409 genes, but each of isolate had less than 30 unique ortholog genes, very small. Unique genes of each isolate were predicted to be involved in fungal life cycles including growth (FluG domain-containing protein, flocculin, methyltransferase, and protein kinase domain-containing protein), stress response (indoleamine-dioxygenase subfamily, tlh3, metallophosphoesterase, and transposase-like protein) and nutrient uptake (alpha/beta hydrolase fold-3 and ribose-phosphate pyrophosphokinase). The rest half of them is hypothetical protein gene.

In GO analysis of ERL-836 DEGs, there were significantly more up-regulated genes than down-regulated genes in various GO categories, when compared with JEF-007. Particularly highly up-regulated ERL836 GO terms were detected in the biological process and molecular function of ERL836. This result indicates that when infecting thrips, possibly ERL836 much easily overcomes the defense mechanisms of thrips compared to JEF-007, but JEF-007 looks having difficulty in infecting the thrips. In ERL836, many fungal genes were highly up-regulated in the infection process. For example, the eukaryotic aspartyl protease showed an 11.2-fold increase in ERL836 when infecting western flower thrips. This aspartyl protease has a potential to degrade host insect tissues and disable anti-microbial peptides [26].

In the GO enrichment analysis, S. cerevisiae was used as a reference rather than B. bassiana because Saccharomyces GO IDs have been established much better than B. bassiana GO IDs and provide a good eukaryotic model system for genetic and biochemical mechanisms. However, no gene IDs were found for 3,576 genes which were shared by both of ERL836 and JEF-007. This is a limitation of this GO enrichment analysis. In spite of this limitation, GO categories of membrane transporters were significantly different between infecting JEF-007 and ERL836. ABC transporters are membrane transporters that play an essential role in B. bassiana infection as well as in infection by the phytopathogenic fungi ABC1 and BcatrB [27]. Transporter capability depends on activation of antioxidative enzymes [28]. B. bassiana harbors 21 ABC transporter genes, but not all ABC transporters contribute to fungal pathogenicity [29].

In the enrichment analysis, fatty acid degradation pathway was significantly different between infecting JEF-007 and ERL836. The Cytochrome P450 genes have been found to be involved in fatty acid degradation. Interestingly, most cytochrome P450 genes in infecting ERL836 were up-regulated, but in JEF-007 they were down-regulated. This suggests that these cytochrome P450 genes of ERL836 could play an important role in thrips infection in contrast to in JEF-007. In eukaryotic microorganisms, cytochrome P450 is involved in alkane oxidation. At least 83 genes encoding cytochrome p450s have been reported in B. bassiana, and a study of B. bassiana infection of silkworms revealed that cytochrome P450s are involved in infection and proliferation of B. bassiana [30]. In the previous study, CYP5099A1, CYP617A1, CYP617A2, CYP52G8, CYP539B1, and CYP655C1 were found to be involved in fatty acid degradation (KO00071), which is related to insect hydrocarbon degradation [31]. Cytochrome P450 enzymes can catalyze cascade formation of mono-oxidation products followed by diterminal oxidation and finally produce α-ω acids [32]. CYP505D4 is involved in sulfur metabolism (KO00920). CYP540B16, CYP61A1, CYP586B1, and CYP682H1 are involved in steroid biosynthesis (KO00100). The two B. bassiana isolates, JEF-007 and ERL836, shared the expression of CYP.

In summary, the two B. bassiana isolates ERL836 and JEF-007 showed similar biological characteristics, such as morphology and virulence against western flower thrips as well as similar genome structure, but gene expression patterns were quite different when infecting western flower thrips. B. bassiana ERL836 appeared to infect thrips easily, while JEF-007 appeared to have more difficulty based on transcriptional analyses. Interestingly, in the enrichment analysis using shared DEGs between two infecting isolates, plasma membrane-mediated transporter activity and fatty acid degradation pathway including cytochrome P450 were more active in infecting ERL836. This comparative approach using shared DEG analysis could be applied to easily characterize the difference of the two patented B. bassiana isolates, ERL836 and JEF-007.

Fungal isolates

B. bassiana ERL836 was obtained from the Entomology Research Laboratory (ERL) Worldwide Collection of Entomopathogenic Fungi of the University of Vermont, USA, and was originally collected from soil in California, USA. B. bassiana JEF-007 was obtained from the Insect Microbiology and Biotechnology Laboratory (IMBL), Chonbuk National University, Republic of Korea. Fungal isolates were grown on quarter strength Sabouraud dextrose agar (1/4SDA; Difco, USA) in the dark at 25 ± 1°C and stored in 20% (v/v) glycerol at -80°C.

Hyphal growth and conidial production

Fungal mycelia masses from 7-day-old B. bassiana isolates were collected into 1.5 ml Eppendorf tubes containing 0.03% (v/v) siloxane solution. After 30 seconds of shaking on a vortex mixer (Vortex-Genie 2TM; VWR Scientific, NY, USA), 2 ml of conidial suspension (1×10⁷conidia ml^-1) was dropped on the middle of a ¼ SDA plate (90 mm diameter). Plates were kept in a 25°C incubator and observed for 9 days under a stereoscopic microscope. Fungal spore productivity was evaluated using millet-based solid cultured granules. Fungal granules were produced by modification of the protocol of Kim et al [33]. Briefly, 200 g of millet was mixed with 100 ml of distilled water with 50% citric acid (160 ml) in a polyvinyl bag. Then the millet was autoclaved at 121°C for 15 minutes and cooled down at room temperature. One milliliter of fungal conidial suspension (1×10⁷conidia ml^-1) was inoculated and incubated at 25 ± 2°C for 10 days. After drying the fungal granules till they had a moisture content of less than 10%, the 0.1 g millet-based medium was collected into 1.5 ml microtubes with 0.03% (v/v) siloxane solution (Silwet, FarmHannong, Seoul, Republic of Korea). The fungal suspension was then vortexed vigorously for 3 min and the number of conidia was counted using a hemocytometer (Cat No. 2960408, Marienfeld, Bad Mergentheim, Germany) under a microscope (400×).

Bioassay against western flower thrips

To compare the virulence of the two B. bassiana JEF-007 and ERL836, a colony of western flower thrips was received from the National Institute of Horticultural and Herbal Science, Republic of Korea. Insects were placed on filter paper moistened with 2 ml of distilled water in a breeding dish (90 mm diameter, 50 mm height). Each stage of thrips was reared in a different breeding dish and provided with fresh bean sprouts every day. Old sprouts with eggs were transferred to new breeding dish every second day. Thrips were kept at 25 ± 2°C, 40 ± 10% relative humidity, and a 14:10 (L:D) photoperiod. For fungal bioassay against the thrips, conidial suspensions were prepared and adjusted to 1×10⁵, 1×10⁶, 1×10⁷, and 1×10⁸conidia ml^-1 using 0.03% (v/v) siloxane solution. A cucumber disc (diameter 60 mm) was placed on a filter paper in a Petri-dish (diameter 60 mm) and 1 ml of conidial suspension was sprayed on the cucumber disc and dried for 10 min at room temperature. Distilled water containing 0.03% Silwet (v/v) served as a non-treated control. To maintain high humidity, 100 µl distilled water was added to the dishes. Eleven two-day-old thrips adults were transferred to the fungus-treated dish (11 adults dish^-1) and coved with lids. Treated dishes were kept at 24 ± 2°C and living adults were counted every day. Each treatment was replicated three times (3 plates per treatment). Lethal concentration 50 (LC₅₀) values of the two isolates were calculated using a Probit analysis (SPSS Inc., 2018).

Whole genome sequencing

For whole genome sequencing of B. bassiana ERL836 at Macrogen (www.macrogen.com; Macrogen Inc., Seoul, Korea), genomic DNA and RNA were extracted from 7-day-old fungal mycelia. DNA quantity was assessed using Pico-green staining (Invitrogen, Cat No. P7589) and Victor 3 fluorometry. To assess DNA quality, gel electrophoresis was performed. The concentration of genomic DNA (81.79 ng/µl) was measured using a Nano Drop spectrophotometer (Thermo Scientific) and a Qubit fluorometer (Life Technology). For PacBio RS sequencing, 8 g of input genomic DNA was used for 20 kb library preparation. For gDNA with a size range less than 17 Kb, a Bioanalyzer 2100 (Agilent) was used to determine the actual size distribution. Genomic DNA was sheared with g-TUBE (Covaris Inc., Woburn, MA, USA) and purified using AMPure PB magnetic beads (Beckman Coulter Inc., Brea, CA, USA) if the apparent size was greater than 40 kb. The gDNA concentration was measured using both a Nano Drop spectrophotometer and a Qubit fluorometer, and approximately 200 ng ml^-1 of gDNA was run on a field-inversion gel. Total of 10 ml of library was prepared using the PacBio DNA Template Prep Kit 1.0 (for 3~10 Kb). SMRT bell templates were annealed using PacBio DNA/Polymerase Binding Kit P6. PacBio DNA Sequencing Kit 4.0 and eight SMRT cells were used for sequencing. Subsequent steps were based on the PacBio Sample Net-Shared Protocol, which is available at http://pacificbiosciences.com/.

Genome analysis

All the genome-level comparisons were outsourced to Macrogen except genome alignment and KEEG analysis. When B. bassiana JEF-007 and ERL836 genomes were compared using a Repeat Masker program (v4.0.5), another B. bassiana ARSEF2860 was used as a reference. To obtain fungal secreted protein proportions, TMHMM (v2.0) was used and the protein sequences of the B. bassiana isolates were subjected to ortholog analysis using OrthoMCL (v.2.0.3). Sequences encoding peptides shorter than 10 amino acids and those with more than 20% stop codons were removed prior to blastp analysis (v2.2.25+; E-value 1E-5). To compare the sequences of two B. bassiana isolates, the Nucmer (NUCleotide MUMmer) module of the MUMmer package was used and data were analyzed using the following three steps: maximal extract matching, match clustering, and alignment extension. Alignments of ERL836 with JEF-007, ERL836 with ARSEF2860, and JEF-007 with ARSEF2860 were analyzed by Assemblytics (http://assemblytics.com/) and dot blot figures were drawn. Lastly KEGG orthology (KO) was analyzed by the web-based server KAAS (http://www.genome.jp/tools/kaas/). For gene ortholog analysis and pathway mapping, query sequences of JEF-007, ERL836 and ARSEF2860 were inputted, respectively. Identification of orthologous genes was conducted by bi-directional best hit (BBH). A total of 348,165 ascomycete sequences (taken from the KEGG databases) were used as reference sequences.

RNA extraction from thrips-infecting B. bassiana isolates

B. bassiana ERL836 or JEF-007 which was infecting adults of western flower thrips was subjected to RNA extraction. Conidial suspension (70 ml, 1×10⁷ conidia ml^-1) of each isolate was spread on ¼SDA plate and cultured at 27°C for 3 days. Then ca. 600 adults of western flower thrips were transferred to the cultured plate and the plate was incubated at 27°C for 3 days. Infected thrips were collected into 1.5 ml microtubes with 0.03% (v/v) siloxane solution, and vortexed for 3 minutes. To exclude thrips bodies, the suspensions were filtered by sterilized iron mesh (ca. 80 mm²pore size) and a fungal suspension was collected. As a control, fungal mycelia of each isolate were harvested from 6-day-old cultures on ¼SDA. The two different treatments were replicated three times (3 samples from the thrips-treated B. bassiana and another 3 samples from the B. bassiana only). Total RNAs were extracted by TRIzol reagent (Invitrogen Life Technologies, CA, USA) following the manufacturer’s instructions. RNA purity and integrity were quantified by ASP-2680 spectrophotometer (ACTGene, Piscataway, NJ, USA) and an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies, Palo Alto, Cambridge, USA).

RNA sequencing and analysis

Libraries of non-infecting and infecting ERL836 or JEF007 were made using the Truseq RNA kit (Illumina, San Diego, USA) following the manufacturer’s protocol at Macrogen. Multiple indexing adapters were ligated to the ends of the double-stranded cDNA and then enriched by PCR to create DNA library templates. In each isolate, infecting and non-infecting samples were sequenced in parallel using an Illumina HiSeq 2000 sequencer. Before analyzing the sequences, quality scores were checked by Fast QC (ver 0.11.7). After receiving sequencing raw data from Macrogen, the following analyses were conducted in our laboratory. For efficient and robust de novo reconstruction of transcriptomes, Trinity (ver 2.8.3) was used (https://github.com/trinityrnaseq/). To identify candidate coding regions within transcript sequences, Trans-Decoder (ver 5.5.0) was used. In the first step, long ORFs were extracted, and then likely coding regions were predicted. To cluster similar nucleotide sequences into clusters meeting a user-defined similarity threshold (0.9 in this analysis), CD-HIT-EST (version 4.7) was used. To quantify transcript abundances, Kallisto (ver 0.45.0) was used to build an index form from the fasta form of target sequences and non-infecting and infecting libraries were compared. Transcripts per million (TPM) of non-infecting and infecting samples was calculated. Raw signals were normalized using a log₂-based transformation. Fold-change statistical tests were performed and log₂|FC|≧2 was defined as statistically significant differential expression. Genes were blasted using the Blast2Go program. Analysis was conducted by local blast with B. bassiana ARSEF2860 as a reference sequence. Statistical significance threshold was 1.0E-10 and the number of blast hits was set to one. Gene ontology (GO) analysis of up- and down-regulated genes was performed using InterPro (online) in the Blast2Go program. The public EMBL-EBI database was used to scan sequences against InterPro's signatures. Up- and down-regulated genes were annotated at GO level 2.

Validation of RNA-sequencing

Non-infecting and infecting RNA samples were subjected to reverse transcription (RT) using AccuPower® RT PreMix (Bioneer, Daejeon, Republic of Korea) with the oligo (dT) 15 primer (Promega, MI, USA). A set of primers for RT-PCR (Supplementary Table S4) were designed at Snap Dragon (http://www.flyrnai.org/snapdragon). RT-PCR was performed as follows: an initial denaturation at 94°C of 5 minutes followed by 34 cycles of 30 seconds at 94°C, 30 seconds at (Tm)°C, and 30 seconds at 74°C, followed by a final extension for 10 minutes at 74°C (C-1000, Bio Rad, Hercules, CA, USA). To validate up-regulated genes, five genes of JEF-007 and three genes of ERL836 were randomly selected and subjected to quantitative RT-PCR (qRT-PCR). To remove gDNA contamination, RNA extracted by the Trizol method (described above) was treated with 1 mg of DNaseI (Invitrogen Life Technologies, CA, USA). cDNAs were generated using the AccuPowder^® Rocketscript RT PreMix kit (BioNeer) and oligo(dT) primer (Promega), and they were used as the template for PCR amplification. Then qRT-PCR was performed using Thunderbird^® Syber^® qPCR mix (QPS-201, TOYOBO, Japan) on a 96-well Bio-Rad CFX96 Real-Time PCR System (Bio-Rad, USA). Cycling parameters for qRT-PCR were as follows: denaturation for 1 minute at 95°C, and then 40 cycles of 15 seconds at 95°C, 1 minute at 60°C followed by melting with an increase in temperature of 0.5°C per 5 seconds starting from 65°C to 95°C. Primers for B. bassiana actin (=γ-actin, GenBank Accession No: HQ232398) were used as an internal control to obtain relative expression levels [34]. △Ct (threshold cycle) was calculated as (Ct value of up-regulated genes) - (Ct value of Bb-actin) and subjected to the calculation of fold change value (2^-^△△^Ct).

GO enrichment analysis

Shared DEG contigs between infecting JEF-007 and infecting ERL836 which had more than 3-fold change difference were subjected to GO enrichment analysis to investigate the different mechanisms of the two fungal isolates in pathogenesis. In this analysis, blastn was performed with an E-value of 1.0E-10. Non-blasted and non-GO ID contigs were removed in this analysis. Functional enrichment was performed using the g:Profiler web server (http://bitt.cs.ut.ee/gprofiler/) [35]. Input query data list was matched with Saccharomyces cerevisiae (https://flybase.org/). The B. bassiana does not provide enough gene IDs for enrichment analysis, so alternatively S. cerevisiae with much larger analyzed gene IDs was used as a reference. The related p-value was corrected for multiple testing using the Benjamini-Hochberg False Discovery Rate (FDR) procedure with a threshold of 0.05.

Data analysis

Data on the numbers of conidia, percentage of live thrips and expression level in qRT-PCR were arc-sine transformed and analyzed using an ANOVA or generalized linear model (GLM) followed by Tukey's honestly significant difference (HSD) for multiple comparisons. All the analyses were conducted using SPSS (SPSS Inc., 2018) at the 0.05 (α) level of significance.

1/4SDA: quarter strength Sabouraud dextrose agar; ANOVA: analysis of variance; DEG: differentially expressed gene; FC: fold change; FDR: false discovery rate; GLM: generalized linear model; GO: gene ontology; HSD: honestly significant difference; KEEG: Kyoto encyclopedia of genes and genomes; LC₅₀: lethal concentration 50%; and LT₅₀: lethal time 50%

Acknowledgements

We appreciate Dr. Woo Jin Kim for comments and advice for the analyses of genome and RNA sequencing data, and National Institute of Horticultural and Herbal Science, Republic of Korea for providing western flower thrips colonies.

Authors’ contributions

JSK designed the whole experiments and wrote the manuscript. SK conducted bioassay against thrips and analyzed the NGS data and wrote manuscript. JCK conducted analyzed the WGS and RNA-seq data using raw sequencing materials. SJL designed the analyses of WGS and RNA-seq data. MRL, SEP, DL and SB reared the thrips and supported the bioassays and statistical data analysis. TYS supported the design of experiments and helped troubleshooting. All authors reviewed this manuscript.

Funding

This work was supported by Korea Institute of Planning and Evaluation for Technology in Food, Agriculture and Forestry (IPET) through Agriculture and Livestock Machinery / Equipment industry Technology Development Program, funded by Ministry of Agriculture, Food and Rural Affairs (MAFRA) (Grant No: 118103-03).

Availability of data and materials

The genome of B. bassiana ERL836 has been deposited in the GenBank database under the project accession no. PPTI00000000.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

¹Department of Agricultural Biology, Chonbuk National University, Jeonju 54596, Korea. ² Department of Microbiology and Cell Science, University of Florida, Gainesville FL 32611-0700, USA.

Lee SJ, Lee MR, Kim S, Kim JC, Park SE, Li D, Shin TY, Nai YS, Kim JS, Genomic Analysis of the Insect-Killing Fungus Beauveria bassiana JEF-007 as a Biopesticide. Sci Rep. 2018;8:12388.
Rehner SA, Buckley EA, Beauveria phylogeny inferred from nuclear ITS and EF1-α sequences: evidence for cryptic diversification and links to Cordyceps teleomorphs. Mycologia 2005;97:84-98.
Glare TR, Inwood AJ, Morphological and genetic characterisation of Beauveria spp. from New Zealand. Mycol Res. 1998;102:250-256.
Rehner SA, Minnis AM, Sung GH, Luangsa-ard JJ, Devotto L, Humber RA, Phylogeny and systematics of the anamorphic, entomopathogenic genus Beauveria. Mycologia 2011;103:1055-1073.
Zhang S, Xia Y, Keyhani NO, Contribution of the gas1 gene of the entomopathogenic fungus Beauveria bassiana, encoding a putative glycosylphosphatidylinositol-anchored β-1, 3-glucanosyltransferase, to conidial thermotolerance and virulence. Appl Environ Microbiol. 2011;77:2676-2684.
Xiao G, Ying SH, Zheng P, Wang ZL, Zhang S, Xie XQ, Shang Y, St Leger RJ, Zhao GP, Wang C, Genomic perspectives on the evolution of fungal entomopathogenicity in Beauveria bassiana. Sci Rep. 2012;2:483.
Meyer E, Aglyamova GV, Wang S, Buchanan-Carter J, Abrego D, Colbourne JK, Willis BL, Matz MV, Sequencing and de novo analysis of a coral larval transcriptome using 454 GSFlx. BMC Genomics 200;10:219.
Shu S, Chen B, Zhou M, Zhao X, Xia H, Wang M, De novo sequencing and transcriptome analysis of Wolfiporia cocos to reveal genes related to biosynthesis of triterpenoids. PLoS One 2013;8:e71350.
Yang P, Zhu JY, Gong ZJ, Xu DL, Chen XM, Liu WW, Lin XD, Li YF, Transcriptome analysis of the Chinese white wax scale Ericerus pela with focus on genes involved in wax biosynthesis. PLoS One 2012;7: e35719.
Mantilla JG, Galeano NF, Gaitan AL, Cristancho MA, Keyhani NO, Gongora CE, Transcriptome analysis of the entomopathogenic fungus Beauveria bassiana grown on cuticular extracts of the coffee berry borer (Hypothenemus hampei). Microbiol. 2012;158:1826-1842.
Xia J, Zhang CR, Zhang S, Li FF, Feng MG, Wang XW, Liu SS, Analysis of whitefly transcriptional responses to Beauveria bassiana infection reveals new insights into insect-fungus interactions. PLoS One 2013;8:e68185.
Liu Y, Shen D, Zhou F, Wang G, An C, Identification of immunity-related genes in Ostrinia furnacalis against entomopathogenic fungi by RNA-seq analysis. PLoS One 2014;9:e86436.
Chu Y, Liu Y, Shen D, Hong F, Wang G, An C, Serine proteases SP1 and SP13 mediate the melanization response of Asian corn borer, Ostrinia furnacalis, against entomopathogenic fungus Beauveria bassiana. J Invertebr Pathol. 2015;128:64-72.
Yang YT, Lee SJ, Nai YS, Kim S, Kim JS, Up-regulation of carbon metabolism-related glyoxylate cycle and toxin production in Beauveria bassiana JEF-007 during infection of bean bug, Riptortus pedestris (Hemiptera: Alydidae). Fungal Biol. 2016;120:1236-1248.
Xing D, Yang Q, Jiang L, Li Q, Xiao Y, Ye M, Xia Q, RNA-seq analyses for two silkworm strains reveals insight into their susceptibility and resistance to Beauveria bassiana infection. Int J Mol Sci. 2017;18:234.
Lai Y, Chen H, Wei G, Wang G, Li F, Wang S, In vivo gene expression profiling of the entomopathogenic fungus Beauveria bassiana elucidates its infection stratagems in Anopheles mosquito. China Life Sci. 2017;60:839-851.
Chen A, Wnag Y, Shao Y, Zhou Q, Chen S, Wu Y, Chen H, Liu E, Genes involved in Beauveria bassiana infection to Galleria mellonella. Arch Microbiol. 2018;200:541-552.
Lee SJ, Yang YT, Kim S, Lee MR, Kim JC, Park SE, Hossain MS, Shin TY, Nai YS, Kim JS, Transcriptional response of bean bug (Riptortus pedestris) upon infection with entomopathogenic fungus, Beauveria bassiana JEF‐007. Pest Manage Sci. 2019;75:333-345.

19 Wang JJ, Bai WW, Zhou W, Liu J, Chen J, Liu XY, Xiang TT, Liu RH, Wang WH, Zhang BL, Transcriptomic analysis of two Beauveria bassiana strains grown on cuticle extracts of the silkworm uncovers their different metabolic response at early infection stage. J Invertebr Pathol. 2017;145:45-54.

20 Zhang Z, Zhang P, Li W, Zhang J, Huang F, Yang J, Bei Y, Lu Y, De novo transcriptome sequencing in Frankliniella occidentalis to identify genes involved in plant virus transmission and insecticide resistance. Genomics 2013;101:296-305.

21 Schneweis DJ, Whitfield AE, Rotenberg D, Thrips developmental stage-specific transcriptome response to tomato spotted wilt virus during the virus infection cycle in Frankliniella occidentalis, the primary vector. Virol. 2017;500:226-237.

22 Shrestha A. Champagne DE, Culbreath AK, Rotenberg D, Whitfield AE, Srinivasan R, Transcriptome changes associated with Tomato spotted wilt virus infection in various life stages of its thrips vector, Frankliniella fusca (Hinds). J Gen Virol. 2017;98:2156-2170.

23 Gupta R, Kwon SY, Kim ST, An insight into the tomato spotted wilt virus (TSWV), tomato and thrips interaction. Plant Biotechnol Rep. 2018;12:157-163.

24 Gamage SMW, Rotenberg D, Schneweis DJ, Tsai CW, Dietzgen RG, Transcriptome-wide responses of adult melon thrips (Thrips palmi) associated with capsicum chlorosis virus infection. PLoS One 2018;13:e0208538.

25 Padmavathi J, Devi KU, Rao CUM, Reddy NNR, Telomere fingerprinting for assessing chromosome number, isolate typing and recombination in the entomopathogen Beauveria bassiana. Mycol Res. 2003;107:572-580.

26 Freimoser FM, Screen S, Hu G, Leger RS, EST analysis of genes expressed by the zygomycete pathogen Conidiobolus coronatus during growth on insect cuticle. Microbiol. 2003;149:1893-1900.

27 Stefanato FL, Abou-Mansour E, Buchala A, Kretschmer M, Mosbach A, Hahn M, Bochet CG, Métraux JP, Schoonbeek HJ, The ABC transporter BcatrB from Botrytis cinerea exports camalexin and is a virulence factor on Arabidopsis thaliana. Plant J. 2009;58:499-510.

28 Wang ZL, Zhang LB, Ying SH, Feng MG, Catalases play differentiated roles in the adaptation of a fungal entomopathogen to environmental stresses. Environ Microbiol. 2013;15:409-418.

29 Song TT, Zhao J, Ying SH, Feng MG, Differential contributions of five ABC transporters to mutidrug resistance, antioxidion and virulence of Beauveria bassiana, an entomopathogenic fungus. PLoS One 2013;8:e62179.

30 Hou C, Qin G, Liu T, Geng T, Gao K, Pan Z, Qian H, Guo X, Transcriptome analysis of silkworm, Bombyx mori, during early response to Beauveria bassiana challenges. PLoS One 2014;9:e91189.

31 Pedrini N, Ortiz-Urquiza A, Zhang S, Keyhani NO, Targeting of insect epicuticular lipids by the entomopathogenic fungus Beauveria bassiana: hydrocarbon oxidation within the context of a host-pathogen interaction. Front Microbiol. 2013;4:24.

32 Scheller U, Zimmer T, Becher D, Schauer F, Schunck WH, Oxygenation cascade in conversion of n-alkanes to α, ω-dioic acids catalyzed by cytochrome P450 52A3. J Biol Chem. 1998;273:32528-32534.

33 Kim JS, Kassa A, Skinner M, Hata T, Parker BL, Production of thermotolerant entomopathogenic fungal conidia on millet grain. J Ind Microbiol Biotechnol. 2011;38:697-704.

34 Zhou YH, Zhang YJ, Luo ZB, Fan YH, Tang GR, Liu LJ, Pei Y, Selection of optimal reference genes for expression analysis in the entomopathogenic fungus Beauveria bassiana during development, under changing nutrient conditions, and after exposure to abiotic stresses. Appl Microbiol Biotechnol. 2012;93:679-685.

35 Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, Vilo J, g:Profiler - a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res. 2016;44:W83-W89.

Table S1. de novo assembly of B. bassiana JEF-007 (A) and ERL846 (B) after the sequencing of whole genomes using Pac Bio RSII technology with error correction.

Table S2. de novo assembly of B. bassiana RNA-sequencing raw data.

Table S3. GO enrichment analysis of differentially expressed and shared genes between infecting B. bassiana JEF-007 and ERL836.

Table S4. Primers used in qRT-PCR for validation of B. bassiana RNA-sequencing.

Figure S1. Validation of RNA-sequencing of B. bassiana ERL836 and JEF-007 using qRT-PCR.

SupplementaryInformationRevised.pdf

Download PDF

Journal Publication

published 27 Nov, 2020

Read the published version in BMC Genomics →

Editorial decision: Major revision
02 Nov, 2020
Review #1 received at journal
09 Aug, 2020
Reviewer #1 agreed at journal
20 Jul, 2020
Editor assigned by journal
16 Jul, 2020
Reviewers invited by journal
16 Jul, 2020
Submission checks completed at journal
15 Jul, 2020
Editor invited by journal
15 Jul, 2020

You are reading this older preprint version

Read the latest preprint version →

Beauveria bassiana ERL836 and JEF-007 with similar virulence show different gene expression when infecting western flower thrips, Frankniella occidentalis

Status:

Journal Publication

Version 2

Abstract

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Supplementary Information

Supplementary Files

Status:

Journal Publication

Version 2