Transcriptomic analysis of Bursaphelenchus xylophilus treated by a potential phytonematicide, punicalagin

Punicalagin (≥98%, HPLC, pomegranate) was purchased from Sigma-Aldrich (St. Louis, MO, USA). Our strain of B. xylophilus was isolated from chips of infected pine wood collected in Nanjing, Jiangsu Province, China using the Baermann funnel technique (Viglierchio and Schmit, 1983). Nematodes were reared on a lawn of Botrytis cinerea cultured on potato dextrose agar (PDA) medium in the dark at 26°C for seven days, at which point the ratio of adult male to adult female to juvenile was approximately 1:1:2. An aqueous suspension of nematodes (about 25 nematodes per μl) was prepared and the nematode suspension mixed with punicalagin (500 μM) as the treatment group (P). Treated (P) and normal (CK) nematode samples were cultured in the dark at 26°C for 36 h and then were collected by centrifugation at 4°C for RNA extraction. All samples were assayed in triplicate.

Total RNA was extracted with TRIzol (Invitrogen, Carlsbad, CA, USA) according to manufacturer’s instructions and was then treated with DNase I (Takara, Dalian, China). RNA quantity and quality were assessed by UV spectroscopy using a NanoDrop spectrophotometer and a Qubit fluorimeter (Thermo Fisher Scientific, Pittsburg, PA, USA), respectively. RNA-Seq transcriptome libraries were constructed using 5 μg total RNA using a TruSeq RNA sample preparation Kit (Illumina, San Diego, CA, USA). mRNAs were isolated using oligo-dT magnetic beads. Fragmentation, cDNA synthesis, end repair, adenylation of 3´ ends, and ligation of the Illumina-indexed adapters were performed according to the manufacturer’s protocol. Libraries were enriched in a 15-cycle PCR reaction and then size-selected using 2% certified Low Range Ultra Agarose (Biorad, Hercules, CA, USA) and quantified using the dye TBS 380 Picogreen (Invitrogen). Clusters were generated by bridge PCR amplification on a cBot System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, libraries were sequenced on the Illumina HiSeq 4000 platform (2 × 151 bp read length). To ensure higher accuracy of the successive bioinformatics analysis, clean reads were obtained (Tables A1, A2) by removing reads containing adapters, reads in which more than 10% of the bases were unknown and low quality reads from the raw data using SeqPrep (https://github.com/jstjohn/SeqPrep; -q 20-L 20) and Sickle (https://github.com/ajoshi/sickle; default parameters).

Transcriptome assembly was carried out using the short-read assembly program Trinity (http://trinityrnaseq.sourceforge.net/, v2013-02-25) using the parameters “--max_memory 50 G --min_kmer_cov 3 --min_contig_length 350 --bfly_opts --V 10” (Haas et al., 2013; Grabherr et al., 2011). Unigenes were compared against those in the NCBI non-redundant protein sequence database and annotated in detail using the BLAST alignment algorithm. Unigenes were annotated by gene ontology (GO, http://www.geneontology.org/) terms describing biological processes, molecular functions and cellular components. BLAST2GO (http://www.blast2go.com/b2ghome) was used to get GO annotations (Conesa and Götz, 2008). The sequences were aligned using the KOG (Eukaryotic Orthologous Groups) database to predict and classify their functions (Grabherr et al., 2011).

Candidate gene pathways were identified and annotated using Kyoto encyclopedia of genes and genomes (KEGG, http://www.genome.jp/kegg/pathway.html) pathway software, KOBAS v. 2.0 (Xie et al., 2011).

The expression level for each transcript was calculated using the fragments per kilobase of exon per million mapped reads method to identify DEGs between the two different samples. The program RSEM (http://deweylab.biostat.wisc.edu/rsem/) was used for differential expression analysis (Li and Dewey, 2011). The genes were selected using a false discovery rate (FDR) ≤0.05 and log2(fold-change) ≥1. In addition, functional-enrichment analysis using GO and KEGG was performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at Bonferroni-corrected from p < 0.05 when compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out using Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do).

The punicalagin-treated (P) and normal (CK) nematode samples were obtained using the same method as those for RNA-Seq. The first-step solution for cDNA synthesis was prepared using 1 μl random hexamers (50 µM), 1 μl dNTPs (at 10 mM each), RNA template (<5 ug), and RNase-free ddH₂O, topped up to a final solution volume of 10 μl, and incubated at 65°C for 5 min, then cooled in an ice bath. The second-step solution was subsequently prepared using 10 μl the first-step reaction solution, 4 μl 5x PrimeScript II Buffer, 0.5 μl RNase Inhibitor (40 U/µl), 1 μl PrimeScript II RTase (200 U/µl) and RNase-free ddH₂O topping up to a final solution volume of 20 μl. The final solution were incubated successively at 30°C for 10 min, 42°C for 30 min, and 95°C for 5 min, and it was then cooled in an ice bath.

For validation, we selected six representative unigenes related to basic physiological and metabolic functions of B. xylophilus, including: DN8346_c0_g1, encoding dynein heavy cytoplasmic; DN18432_c0_g1, encoding ATP synthase F0 subunit 6 (mitochondrion); DN6713_c0_g1, encoding twitchin; DN14208_c0_g1, encoding small HSP21-like protein; DN8075_c0_g1, encoding a nematode cuticular collagen (collagen triple helix repeat domain containing protein); DN14349_c1_g1, encoding glucosidase 2 subunit beta. A PrimeScript II 1st Strand cDNA Synthesis Kit (Takara, 6210B) was used for first-strand cDNA synthesis experiment. The PCR primers were designed using Primer Premier 5.0 (Table 1). The actin-4 gene of B. xylophilus served as an internal control. Quantitative real-time PCR was carried out using SYBR Premix Ex Taq II (Takara). Relative gene expression values was calculated using the 2^−ΔΔC _T method (Livak and Schmittgen, 2001; Revel et al., 2002), and the Pearson correlation coefficient was computed in R.

We identified 21,750 unigenes with a mean length of 1,492 bp and an N50 value of 2,654 bp through de novo transcriptome assembly. The longest unigene was 19,867 bp and the smallest was 351 bp. Unigenes with lengths between 400 and 5,000 bp accounted for approximately 90.29% of the total (Fig. 1). We functionally annotated 12,608 (57.97%) of the unigenes (Fig. 2). The NCBI non-redundant (NR) database contained 12,587 (57.87%) of the genes and 3,752 (17.25%) of these were assigned with the GO terms including biological process, cellular component and molecular function (Fig. 3a). We also classified 2,739 (12.60%) of the genes into 26 functional families (Fig. 3b) and 7,064 (32.48%) were annotated to 321 pathways in five KEGG categories (Fig. 3c).

Figure 1:

Figure 2:

Figure 3:

We found 2,575 DEGs comprising 1,428 up-regulated and 1,147 down-regulated genes. These genes were related to the basic physiologic and metabolic functions of B. xylophilus (Table 2). Our GO functional analysis organized 381 of the DEGs into three groups containing 271, 48 and 62 members that were assigned to 23 categories of biological process, 18 categories of cellular component, and 12 categories of molecular function, respectively. The most highly enriched genes fell into the following categories: those belonging to biological process were single-organism process, cellular process, metabolic process, developmental process, biological regulation, regulation of biological process, localization, locomotion, reproduction and response to stimulus; those belonging to cellular component were cell, cell part, organelle, membrane, organelle part, macromolecular complex and membrane part; those belonging to molecular function were binding and catalytic activity. GO annotation results indicated that punicalagin mainly changed the expression of the genes related to metabolism, development, biological regulation, localization, locomotion and stimulus response that could contribute to its nematotoxic activity (Figs 4a,b).

Figure 4:

KOG analysis identified 499 DEGs that were annotated to 26 categories (Fig. 4c). The top 10 KOG categories were: general function (61 DEGs), signal transduction (58 DEGs), energy production and conversion (47 DEGs), translation, ribosomal structure and biogenesis (45 DEGs), posttranslational modification, protein turnover, chaperones (40 DEGs), intracellular trafficking, secretion and vesicular transport (37 DEGs), transcription (32 DEGs), cytoskeleton (31 DEGs), RNA processing and modification (26 DEGs), and amino acid transport and metabolism (23 DEGs). Genes related to signal transduction, energy metabolism, translation, ribosomal structure and biogenesis, posttranslational modification, protein turnover and chaperones might be the main nematotoxic targets of punicalagin.

The KEGG annotation aligned the DEGs to 279 pathways. The top 10 pathways that were related to life activities of B. xylophilus included oxidative phosphorylation (38 DEGs), endocytosis (35 DEGs), spliceosome (31 DEGs), ribosome (29 DEGs), RNA transport (28 DEGs), carbon metabolism (27 DEGs), ubiquitin mediated proteolysis (24 DEGs), MAPK signaling pathway (24 DEGs), peroxisome (24 DEGs) and biosynthesis of amino acid (22 DEGs).

We selected six genes for qRT-PCR analysis to validate the expression profiles obtained by RNA-Seq. The correlation between these two methods was 100% (Fig. 5).

Figure 5: