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 s of shaking on a vortex mixer (Vortex-Genie 2TM; VWR Scientific, NY, USA), 2 ml of conidial suspension (1 × 107 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 μl) in a polyvinyl bag. Then the millet was autoclaved at 121 °C for 15 min and cooled down at room temperature. One milliliter of fungal conidial suspension (1 × 107 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 × 105, 1 × 106, 1 × 107, and 1 × 108 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 (LC50) 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 μl− 1 of gDNA was run on a field-inversion gel. Total of 10 μl 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 interacting with cuticle of western flower thrips adults was subjected to RNA extraction. Conidial suspension (70 μl, 1 × 107 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 min. To exclude thrips bodies, the suspensions were filtered by sterilized iron mesh (ca. 80 μm2 pore size) and a fungal suspension was collected. As a control, fungal mass of each isolate was 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 cuticle-interacting and cuticle-non-interacting 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, cuticle-interacting and cuticle-non-interacting 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 log2-based transformation. Fold-change statistical tests were performed and log2|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
Cuticle-interacting and cuticle-non-interacting 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 min followed by 34 cycles of 30 s at 94 °C, 30 s at (Tm)°C, and 30 s at 74 °C, followed by a final extension for 10 min 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 μg 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 min at 95 °C, and then 40 cycles of 15 s at 95 °C, 1 min at 60 °C followed by melting with an increase in temperature of 0.5 °C per 5 s 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 JEF-007 and 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.