Pine wood nematode (PWN),
In recent years, we have studied the effects of phytonematicides against
Transcriptomic analysis is essential to interpret the functional elements of the genome and understand development and disease (Wang et al., 2009). Transcriptomic profiling allows the comparison of transcriptomes across a disease state compared to normal cells or of specific experimental stimuli compared to normal physiological conditions. It is thus a potential tool for interpreting the functional elements of the genome and uncovering the biological mechanisms of development and diseases (Wang et al., 2009; Han et al., 2015). In the present study, we performed transcriptomic analysis of
Punicalagin (≥98%, HPLC, pomegranate) was purchased from Sigma-Aldrich (St. Louis, MO, USA). Our strain of
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 (
Data statistics of raw reads from six samples for RNAseq.
Sample | Raw reads | Raw reads base (bp) | Q20 (%) | Q30 (%) |
---|---|---|---|---|
CK1 | 31,003,104 | 4,650,465,600 | 95.72 | 92.01 |
CK2 | 37,538,058 | 5,630,708,700 | 95.1 | 91.11 |
CK3 | 37,726,236 | 5,658,935,400 | 95.31 | 90.95 |
P1 | 35,871,190 | 5,380,678,500 | 96.92 | 93.53 |
P2 | 46,131,806 | 6,919,770,900 | 96.48 | 92.26 |
P3 | 33,488,164 | 5,023,224,600 | 97.32 | 93.81 |
Note: CK represents the normal pine wood nematodes, and P represents the punicalagin-treated pine wood nematodes.
Data statistics of clean reads from six samples after quality control.
Sample | Clean reads | Clean reads base (bp) | Q20 (%) | Q30 (%) |
---|---|---|---|---|
CK1 | 28,585,924 | 4,136,450,402 | 98.48 | 95.62 |
CK2 | 34,085,674 | 4,903,647,270 | 98.36 | 95.36 |
CK3 | 34,177,082 | 4,955,215,491 | 98.19 | 94.97 |
P1 | 33,854,004 | 4,938,319,926 | 98.6 | 95.91 |
P2 | 43,447,298 | 6,287,507,779 | 98.11 | 94.77 |
P3 | 32,165,338 | 4,669,426,049 | 98.45 | 95.57 |
Transcriptome assembly was carried out using the short-read assembly program Trinity (
Candidate gene pathways were identified and annotated using Kyoto encyclopedia of genes and genomes (KEGG,
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 (
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 ddH2O, 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 ddH2O 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
Primer sequences used for qRT-PCR validation of differentially expressed genes.
Gene ID | Primer sequence | Amplicon length (bp) |
---|---|---|
DN8346_c0_g1 | Forward: CTGCTGAATGAGTGGGTA | 197 |
Reverse: AGAAGTTTGAAAGGAGGC | ||
DN18432_c0_g1 | Forward: AAGTGTCACCTCCTTTAC | 132 |
Reverse: GGTATTCGTTTTGTCCT | ||
DN6713_c0_g1 | Forward: CGAGGTCCGTTAGAAGTG | 128 |
Reverse: TTGCCAGTCTCAGTGTCC | ||
DN14208_c0_g1 | Forward: GTAAGCCTGGAGAAAAG | 195 |
Reverse: AGTTGACGGTGTTGGTG | ||
DN8075_c0_g1 | Forward: GCAACCACCAGGAGCAAC | 88 |
Reverse: CGGAAATGATGGAGAACCC | ||
DN14349_c1_g1 | Forward: AAGTGACGAGCCAGGTA | 168 |
Reverse: TCACAAACATTCGGACA | ||
DN6594_c0_g1a | Forward: CAACCCCAAGGCTAACA | 303 |
Reverse: TCACGCACGATTTCACG |
Note: arepresents internal control.
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).
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
Annotated genes differentially expressed in response to punicalagin and related to physiological processes in
TRINITY_Gene ID | Annotation | Log2 FC | FDR | Type |
---|---|---|---|---|
DN8346_c0_g1 | Cytoplasmic dynein heavy chain | 3.0314 | 2.54E-24 | Up |
DN18432_c0_g1 | ATP synthase F0 subunit 6 (mitochondrion) | −3.5432 | 8.44E-20 | Down |
DN6713_c0_g1 | Twitchin | 2.6719 | 1.19E-18 | Up |
DN2658_c0_g1 | Cytochrome c oxidase subunit I (mitochondrion) | −2.2625 | 1.66E-12 | Down |
DN11917_c0_g2 | Cytochrome b, partial (mitochondrion) | −2.3332 | 3.83E-09 | Down |
DN9703_c0_g1 | Cytochrome c oxidase subunit 3 (mitochondrion) | -2.2922 | 3.88E-09 | Down |
DN1325_c0_g2 | Heat shock protein 20 | −2.3758 | 1.70E-05 | Down |
DN14208_c0_g1 | Small HSP21-like protein | −1.4702 | 6.73E-05 | Down |
DN14130_c0_g1 | Heat shock protein Hsp-12.2 | −1.3353 | 0.0005945 | Down |
DN10040_c0_g1 | Electron-transfer-flavoprotein | −1.1439 | 0.01728 | Down |
DN8075_c0_g1 | Nematode cuticle collagen and collagen triple helix repeat domain containing protein | −1.1863 | 0.01908 | Down |
DN14349_c1_g1 | Glucosidase 2 subunit beta | −1.0567 | 0.02580 | Down |
DN22_c0_g1 | NADH dehydrogenase subunit 1 (mitochondrion) | 2.2809 | 0.04550 | Up |
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
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).
Through pathway analysis, we identified gene regulation patterns correlated with a response to the nematoxin punicalagin. A significantly up-regulated unigene was that encoding cytoplasmic dynein heavy chain, a putative factor in phagosome maturation of phagosome pathway (Fig. 6a). This protein has been reported to be involved in insulating and decomposing exogenous substances (Rabinovitch, 1995). These clues lead us to speculate that
We found that the genes encoding proteins involved in the respiratory electron-transport chain including NADH dehydrogenase, cytochrome c oxidase, cytochrome b and flavoproteins were differently expressed under punicalagin treatment. In addition, KEGG annotation identified additional DEGs involved in oxidative phosphorylation (Fig. 6b), an essential process for basic metabolism (Hatefi, 1985; Gray and Winkler, 1996). Furthermore, energy production and conversion in the KOG database contained the third most DEGs (Fig. 4c). All of these clues indicate that punicalagin may influence the energy metabolism of
Lastly, the unigene encoding twitchin was also significantly up-regulated (shown in Table 1). Twitchin is a silk protein related to muscle function in mollusks (Funabara et al., 2007), coinciding with our observation that the PWNs treated with punicalagin under a microscope are abnormally twisted compared with controls (Fig. A3). In addition, the unigene encoding nematode cuticle collagen was down-regulated and this may be related to our previous finding that the body walls of punicalagin-treated nematodes were abnormal (Guo et al., 2017).
In conclusion, we found 2,575 DEGs in total between the two libraries of punicalagin-treated PWNs and control from our comparative transcriptomic analysis. Specifically, we obtained the main DEGs closely related to life activities of PWN, GO terms, KOG calories and KEGG pathways containing more DEGs, and speculated that PWNs could give response to the stimulus from punicalagin through phagosome, endocytosis, peroxisome and MAPK signaling pathways. In addition, punicalagin was supposed to influence PWN energy metabolism. The genes encoding twitchin and nematode cuticular collagen could be crucial regulation targets of punicalagin, which might own to its nematotoxic activity against PWN. Together, these clues may prove valuable for elucidating the nematotoxic mechanism of punicalagin against