Diapause induction, termination, sequencing, and gene identification
After the completion of the diapause induction process, the parasitized host eggs were dissected to verify whether T. dendrolimi entered diapause successfully. Diapausing parasitoids remain at the prepupal stage (Dpre). If diapause induction fails, the prepupae would die or continue to develop into pupae and adults. Following the diapause termination process, the parasitized host eggs were transferred to normal development conditions (26°C ± 1°C, 60% ± 5% RH, 16:8 h L:D). If diapause was disrupted, the parasitoids would transform from prepupa to pupa within several days, noted as Dp. However, if diapause was not disrupted, the T. dendrolimi would remain in the prepupa stage. In this study, 99% of the parasitoids entered diapause successfully, and about 95% resumed development after the 70 d termination treatment. The prepupae and pupae of T. dendrolimi that developed under normal conditions were obtained as reference groups noted as NDpre and NDp, respectively.
RNA samples obtained from distinct stages of T. dendrolimi were prepared and sequenced using the Illumina Hiseq2000 sequencing platform. Four cDNA libraries were constructed from the samples of Dpre, Dp, NDpre, and NDp. After filtering raw reads (reads containing adaptors, reads containing N larger than 10%, and low-quality reads (Qphred < 20) were removed), clean reads were retained (Table 1). The clean data were assembled by Trinity and Corset with 87,022 transcripts, and an average length of 1,604 nt and an N50 of 3,148. Of the transcripts, 35,231 (40.5%) were longer than 1000 bp (Table 2).
To study gene function, transcripts were annotated using BLASTX searches against the non-redundant (NR) sequence database; 39,969 (45.92%) displayed homology to known proteins (E < 1e-5; Fig. 1A). Nearly 25,000 annotated transcripts, over 65% of the annotated transcripts, were homologous to T. pretiosum, probably because the genome of T. pretiosum was the only available Trichogramma. Fewer transcripts were homologous to Nasonia vitripennis (1,217, 3.1%), Apis florea (26, 0.03%), A. dorsata (12, 0.01%), or A. cerana (15, 0.01%) (Fig. 1B). Fewer than 40 transcripts matched those from Microplitis demolitor. Among all annotated transcripts, 73.0% had significant homology with an E-value of <10-30 (Fig. 1C), and 52.3% had a similarity greater than 80.0% (Fig. 1D). After filtering and removing redundant sequences, we retained those with significant BLAST hits and constructed a non-redundant dataset containing 7,593 unigenes with an average length of 3,351 nt. Based on the annotation, such as gene length, ID, and speculative function, the diapause-related genes and potential genes involved in diapause were sorted out for further analysis.
Identification of DEGs and functional classification
Ten genes were selected for validation with qPCR, and glyceraldehyde phosphate dehydrogenase (GAPDH) was selected as the reference gene after measuring its stable expression level in diapause and non-diapause groups. The tendencies of the expression profiles of these genes were similar according to RNA-Seq and qPCR (Fig. 2). Among these 10 selected genes, all except trehalase (tre) were up-regulated during diapause.
To study diapause-specific transcriptional changes in T. dendrolimi induced by low temperature, we made pairwise comparisons between different libraries to identify the DEGs. A total of 5,702 DEGs were identified among four groups. Among these DEGs, there were 3,182 DEGs changed in Dpre compared to NDpre. DESeq identified 3,251 and 3,442 DEGs exclusively changed in Dp vs NDp and Dpre vs Dp, respectively. In addition, the DEGs changed in NDpre vs NDp were 1,511. This group of DEGs may be the genes related to normal development, namely from prepupa to pupa, not to diapause development. According to the Venn diagram, there were 463 genes changed throughout the diapause development process, while in the normal development process, the expression of these genes did not change (Fig. 3).
To determine the potential function of identified DEGs, GO enrichment was performed. In all combinations, except for Dpre vs Dp, more genes were up-regulated. However, when we compared Dpre to Dp, there was little difference in the number of up- and down-regulated genes. Furthermore, more DEGs were assigned to the same category among different groups. Regulation of transcription, DNA-templated process, oxidation-reduction process, and signal transduction process were the top three in these four groups. The number of DEGs involved in ribosome biogenesis was much higher during diapause development than during normal development (Fig. 4). The subsequent analyses are based directly upon these results.
Comparative analysis of genes involved in diapause
Based on the results of GO enrichment, we focused on the genes enriched in the oxidation-reduction process, regulation of transcription, DNA-templated process and signal transduction process, which were processes enriched in most DEGs. In addition, we also examined ribosome biogenesis.
A total of 342 genes were identified in the oxidation-reduction process, and 16 of these belong to cytochrome P450s (CYP450s). In the T. dendrolimi transcriptome, 22 CYP450s were identified, and 16 were differentially expressed. These 16 genes belonged to 4 clans. In diapause stages (Dpre), 10 genes were up-regulated. In the pupae after diapause (Dp), five genes were highly expressed. Only one gene (CYP9E2) was highly up-regulated in normal pupae (NDp) (Fig. 5). These results show that the number of up-regulated genes during diapause was significantly higher than that in other stages.
There were 36 transcription factors differentially expressed during diapause development, and it appears that three kinds of transcription factors might be associated with diapause in T. dendrolimi.
The first kind is zinc finger protein. Three genes encoded zinc finger protein. Zinc finger protein 271, zinc finger 184, and zinc finger 544, were identified in the transcriptome. They were all up-regulated such that the expressions of these three genes in Dpre were higher than expression in Dp. Zinc finger protein gene 271 had an SFP domain. Genes containing this domain are putative transcriptional repressors during the G2/M (second gap period to mitotic period) transition. The wee1 gene, encoding an inhibitory kinase, was up-regulated during diapause in N. vitripennis [36]. We obtained similar results in T. dendrolimi (Fig. 6). In addition, zinc finger protein 184 contained a GDT1 domain, which is a putative Ca2+/H+ antiporter. Ca2+/H+ antiporter, which maintains homeostasis, has been studied in plants, but there are few studies.
The second type of transcription factor is homeobox domain protein. In the transcriptome, 11 homeobox-containing genes were differentially expressed during diapause except for pit1, which was significantly up-regulated in the individuals that terminated diapause. Among these genes, homeobox protein homothorax (hth) had the greatest change in expression, followed by homeotic protein distal-less (dll) and homeobox protein six1 (six1). The gene expression of another homeobox protein six1 was similar to hth; expression increased after entering diapause stage. And this gene may be involved in the regulation of cell proliferation, apoptosis, and embryonic development.
The third group of transcription factors is forkhead box protein. FOXOs have been identified as candidates for the molecular control of embryonic diapause in some species, like Culex pipiens [22, 37]. In T. dendrolimi, three forkhead box proteins (foxo), forkead box protein E3 (foxe3), and forkhead box protein D3 (foxd3) were identified. These genes were up-regulated both in diapause prepupae and resulting pupae, and they likely play a role in diapause development.
In the T. dendrolimi transcriptome, the expression of Protein phosphatase 2A (PP2A), which belonging to signal transduction process, was down-regulated during the diapause stage. This result was consistent with that obtained in the cotton bollworm (Helicoverpa armigera). Low PP2A expression in diapause individuals contribute to the accumulation of p-Akt, and p-Akt leads to H. armigera diapause [24, 38].
In addition to these three biological processes, ribosome biogenesis is also important in the control of cell growth and division in eukaryotes [39]. In this study, ribosome biogenesis involved 31 DEGs, and 29 genes were up-regulated during prepupal diapause. Only two genes, 40S ribosomal protein S11 (rpS11) and 28S ribosomal protein S5 (rpS5) were down-regulated during prepupal diapause. All of the 60S ribosomal proteins were up-regulated (Fig. 7).
Some genes, even those not involved in the biological process, with a considerable number of genes enriched, appear to be important in the diapause development of T. dendrolimi, such as p53 and the DNA damage-regulated gene 1 (pdrg1), which gene expressions were significantly changed during diapause development. The transcriptional expressions of Glutathione-S-transferase (GST) and UDP-glucuronosyltransferase (UDPGT) were also up-regulated during T. dendrolimi prepupal diapause.
Lipid metabolism is essential for energy homeostasis. Some diapausing insects use lipids for energy storage [40, 41]. During diapause, almost all selected lipid metabolism related genes were up-regulated, coinciding with the mobilization of TAG reserves (Fig. 8).