When neonate larvae were reared on an artificial diet with different low concentrations of EO and D-camphor for 3 d, significant deaths initially occurred at 12 mg·g− 1 and 20 mg·g− 1, respectively (P < 0.05, Fig. 1a). Therefore, the LC00 values of the EO and D-camphor treatments were 10 mg·g− 1 and 17 mg·g− 1, respectively. As shown in Fig. 1b, the unique working concentrations of EO and D-camphor based on their calculated LC15 values used in subsequent bioassays were 30 mg·g− 1 and 41 mg·g− 1, respectively.
There were variations in developmental time, body weight, and larval viability between the terpenoid-exposed and control groups (Fig. 2). The developmental time of larvae feeding on the terpenoid-infused diets was longer than that of the control group when they reached the 2nd instar. However, a significant increase in the growth rate of larvae exposed to terpenoids at LC00 (EO_10mg_g & D_camphor_17mg_g) was observed afterwards (Fig. 2a). Body weight was either unaffected or slightly inhibited at LC15 (EO_30mg_g & D_camphor_41mg_g), and we even found that the highest weight gain was observed at LC00, with a gain 10% greater than when terpenoid was absent in terms of pupal weight (Fig. 2b). Terpenoid exposure significantly affected the age-specific survival rate (log-rank test, chi-square = 247.772, P < 0.001). The survivorship curves showed that the larvae reared on terpenoid-infused diets had rapidly declining survivorship from 9 to 24 days after egg hatching, and increasing doses of EO and D-camphor markedly affected larval viability (Fig. 2c). Our final step was to combine the developmental parameters described above into a single index (i.e., the RPI) that serves as a proxy of overall performance. According to the heatmap, the populations reared on terpenoid-infused diets at lower dosages had greater z-scores than the control group, whereas the lowest value occurred at higher dosages (Fig. 2d).
Transcriptome responses to EO and D-camphor exposure
RNA-seq of P. tsushimanus whole larvae fed terpenoid-infused or terpenoid-free diets yielded a total of 794,599,342 high-quality, clean reads, with a 42.94% − 43.62% GC content and over 96% and 91% Q20 and Q30 values, respectively (Table S1). De novo assembly of transcriptome data resulted in 34,486 unigenes and 44,513 transcripts. The total and average lengths of these unigenes were 41,545,312 bp and 1204.7 bp, respectively, with an N50 (length N for which 50% of all bases in the assembly were located in a unigene of length < N) of 2,245 bp and a mapped fragment percentage of 59.75% (Table S2).
A total of 34,452 unigenes were annotated by comparing them to the COG, GO, KEGG, NR, Swiss-Prot or Pfam databases (Table S3). As shown in Figure S1, NR, GO and KEGG annotation provided comprehensive information on gene function. The species distribution of the top BLAST hits in the NR database showed that 6303 unigenes (34.1%) had the highest scoring blast matches to Dendroctonus ponderosae, followed by Sitopjilus oryzae (3434 unigenes, 18.6%) (Fig. S1a). The observed GO terms for each domain (biological process, cellular component, and molecular function) and KEGG pathways for each category (metabolism, human diseases, organismal systems, cellular processes, environmental information processing, and genetic information processing) are detailed in Figure S1b and c.
To screen the genes potentially involved in P. tsushimanus resistance to host-specific terpenoid defence (i.e., EO and D-camphor), differentially expressed genes (DEGs) were detected between the terpenoid-exposed and nonexposed groups (CK_vs_EO_10mg_g; CK_vs_EO_30mg_g; CK_vs_D_camphor_17mg_g; CK_vs_D
_camphor_41mg_g), with the criteria P-value < 0.05 and |log2FC| ≥ 1. In total, 14 upregulated genes and 26 downregulated genes were observed in the EO_10mg_g treatment, whereas 104 upregulated genes and 99 downregulated genes were observed in the EO_30mg_g treatment (Fig. 3a). Additionally, 116 upregulated genes and 48 downregulated genes were observed in the D_camphor_17mg_g treatment, whereas 143 upregulated genes and 68 downregulated genes were observed in the D_camphor_41mg_g treatment (Fig. 3b). The Venn diagram showed that 3, 111, 120, and 50 specific DEGs (i.e., DEGs occurring in only a particular treatment) were found in the EO_10mg_g, EO_30mg_g, D_camphor_17mg_g, and D_camphor_41mg_g treatments, respectively (Fig. 4c). In addition, 18 common DEGs (i.e., DEGs occurring in all treatments) were significantly different between the terpenoid-exposed and nonexposed groups (Fig. 3c).
To annotate the putative functions of the DEGs identified in each treatment, GO enrichment analysis was performed (Fig. 4). In brief, terpenoid exposure induced similar molecular responses in all groups: most upregulated unigenes were significantly enriched in monooxygenase activity related to CYP450s and structural constituent of cuticle related to CPs, while downregulated unigenes significantly enriched in protein folding and de novo posttranslational protein folding were only observed in the EO_10mg_g treatment (Fig. 4a). Subsequently, we placed emphasis on two GO pathways: monooxygenase activity (GO: 0004497) and structural constituent of cuticle (GO: 0042302). CYP450 genes accounted for 87.6% of the monooxygenase activity (Fig. 5a). Among these genes, several CYP450 genes were significantly upregulated after terpenoid exposure, notably the genes called ‘TRINITY_DN986_C0_g1’, ‘TRINITY_DN32261_C0_g1’, and ‘TRINITY_DN1 2993_C0_g1’ (Fig. 5b-e). Analogous to the results for the monooxygenase activity pathway, among the genes associated with structural constituent of cuticle pathway, 88.9% belonged to the CP gene family (Fig. 6a). We also found that most of the CP genes were commonly upregulated in all treatments, while a few genes were found to be downregulated after EO treatments, including the genes called ‘TRINITY_DN3873_C0_g1’ and ‘TRINITY_DN17856_C0_g1’ (Fig. 6b-e).
WGCNA was used to further identify the target genes involved in countering host-specific terpenoid defence. As shown in the gene dendrogram, the genes with similar expression patterns were globally divided into 13 unique module eigengenes using average linkage clustering (Supplementary Fig. S3). The relevance between each module and terpenoid-exposure treatment showed that the MEblack module had a higher correlation with D-camphor treatment than other modules (MS-value = 0.441), whereas the MEtan module was the most relevant to EO treatment (MS-value = 0.787) (Fig. 7a). Next, the top 30 genes with the highest connectivity (i.e., hub genes) in the MEblack and MEtan modules were used to analyse gene expression and co-expression networks (Fig. 7b-g). We found that the hub genes in the MEtan module exhibited higher expression levels in EO treatments than in D-camphor treatments (Fig. 7b & c). With respect to hub genes in the MEblack module, there was no obvious difference in expression level between the EO and D-camphor groups (Fig. 7e & f). Using co-expression networks, five core hub genes in the MEtan and MEblack modules were identified (Fig. 7d & g). Remarkably, the core hub genes in the MEtan module, namely TRINITY_DN986_C0_g1, TRINITY_DN32261_C0_g1, TRINITY_DN12993_C0_g1, and TRINITY_DN1904_
C0_g1, belong to the CYP450 family, whereas the gene TRINITY_DN1820_C0_g1 is of unknow function (Fig. 7d). In the MEblack modules, the core hub genes, namely, TRINITY_DN39845_C0_g1, TRINITY_DN1984_C0_g1, TRINITY_DN624_C0_g1, and TRINITY_DN5193_C0_g1, were annotated to the cuticular protein family, whereas the gene TRINITY_DN129_C0_g1 was uncharacterized (Fig. 7g).