Diapause induction, termination, sequencing, and gene identification
After the induction process was completed, the parasitized host eggs were dissected to verify whether T. dendrolimi entered diapause successfully or not. Diapausing parasitoids should stay at the prepupal stage. If diapause induction failed, they would continue to develop into pupae or 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 L:D). If diapause was disrupted, the parasitoid would convert from prepupa to pupa within several days; conversely, if diapause was undisrupted, the T. dendrolimi would remain in the prepupa stage. In this study, 99% of the parasitoids entered diapause successfully, and about 95% broke out of the diapause state during the 70-day termination treatment.
RNA samples obtained from distinct stages of this species were prepared and sequenced using the Illumina Hiseq2000 sequencing platform. Four cDNA libraries were constructed from the samples described above (Dpre, Dp, NDpre, and NDp). After filtering raw reads (reads containing adaptor, reads containing N larger than 10%, and low-quality reads (Qphred < 20) were removed), clean reads were retained (Table 3). The clean data were assembled and hierarchically clustered by Trinity and Corset with 87,022 unigenes produced, and an average length of 1,604 nt and an N50 of 3,148. Of the unigenes, 35,231 (40.5%) were longer than 1000 bp (Table 4).
In order to gain information of gene function, transcripts were annotated using BLASTX searches against the non-redundant (NR) sequence database; 39,969 (45.92%) of them displayed homology to known proteins (E < 1e-5; Fig. 1A). Nearly 25,000 annotated unigenes, over 65% of the annotated unigenes, were homologous to T. pretiosum, probably because the genome of T. pretiosum was the only available genome in Trichogramma. Fewer unigenes 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). Less than 40 unigenes were matched to those from Microplitis demolitor. Based on these results, we can infer that there are significant genetic differences between Trichogramma and other parasitoid species. Among all annotated unigenes, 73.0% had significant homology with an E-value of <10-30 (Fig. 1C), and 52.3% had a similarity higher than 80% (Fig. 1D). After filtering and removing redundant sequences, we kept the unigenes with significant BLAST hits and constructed a non-redundant dataset containing 7,593 sequences with an average length of 3351 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 GAPDH was selected as the reference gene after measuring its stable expression level for diapause and non-diapause. The tendencies of these genes’ expressions were similar according to RNA-Seq and qPCR (Fig. 2). Among these 10 selected genes, all except tre were up-regulated during diapause.
To get more acquainted with diapause-specific transcriptional changes in T. dendrolimi induced by low temperature, we carried out pairwise comparisons between different libraries to identify the DEGs. The result of DEG analysis revealed that more unigenes showed striking changes in mRNA levels during diapause stage. DESeq2 identified 1,204 DEGs exclusively changed in Dpre compared to NDpre, Dp, and NDp. We used the expressions of these genes at four time points to construct a gene expression profile composed of two groups: 820 genes were up-regulated and thus belonged to group 1, and the others were down-regulated significantly and thus belonged to group 2 (Fig. 3). All these genes might be involved in diapause maintenance of T. dendrolimi.
To figure out the potential function of identified DEGs, GO enrichment was performed. Compared with non-diapause individuals, the number of up-regulated genes in the diapause development process was significantly higher than down-regulated genes. In the gene repertoire of Dpre, more DEGs were assigned to protein binding, oxidation-reduction processes, and metal ion binding. Among them, metal ion binding and ribosome biogenesis were concerned in the subsequent analysis (Fig. 4).
Comparative analysis of genes involved in diapause
Diapause is a dynamic process accompanied by a series of physiological transitions. Several studies have focused on the general gene expression pattern of insect diapause without clear elucidation of the diapause mechanism. This is due to the complexity of the diapause process as well as the variations among insect species.
We identified 79 genes potentially involved in T. dendrolimi diapause based on their homology to genes with known functions from other insect species. Several new genes with the capability of being included in the diapause genetic toolkit were also analyzed. In the following section, they are organized into different categories according to their function.
According to the results of GO enrichment, we focused on the genes enriched in the metal ion binding term. There were eight genes coded zinc finger protein were identified in the transcriptome. Zinc finger protein gene 271 had an SFP domain (Fig. 5). Genes containing this domain were considered to be a putative transcriptional repressor during G2/M transition. One of the most noticeable characteristics of diapause is the blockage of ontogeny, and this blockage always occurs with the cell cycle cessation [38, 39]. In Nasonia vitripennis, the S phase of the cell cycle disappeared in the beginning stage of diapause due to the cells being arrested in the G0/G1 and G2 phases . Similarly, in drosophilid fly, Chymomyza costata, the cell cycle of CNS cells was arrested in the G0/G1 (86.6%) and G2 (12.8%) division phases during diapause. The wee1 gene, coding for a kind of inhibitory kinase, was up-regulated during diapause in this species . We obtained similar results in T. dendrolimi (Fig. 5). Therefore, the wee1 gene might be a molecular marker of diapause.
It was reported that cyclin-dependent kinase 4 (CDK4) and cyclin D regulated the G1/S transition during cell cycle . The expressed level of CDK4 and cyclin D played an essential role in cell cycle arrest during diapause, while in T. dendrolimi, the expression of these two genes were up-regulated during diapause. We speculate that the cell cycle might be stuck at other phases, like G2 or G0. In addition, cyclin A1/G/G2/K was up-regulated during diapause. This suggested that these cyclin-related genes might play a part in T. dendrolimi diapause maintenance.
p53 induces growth arrest or apoptosis and has a negative regulation on cell division by controlling an array of genes in this process . In addition, p53 and DNA damage-regulated gene 1 (pdrg1) are involved in multiple cellular, such as apoptosis, DNA damage repair, cell cycle. In Artemia sinica, PDRG1 takes an essential role in diapause termination and regulation of cell cycle during early embryonic development . By blocking the expression of pdrg1 in human colon cancer cells, cell growth reduced significantly . Transcriptome based analysis indicated that the expression of p53 was significantly increased during diapause, and the contrary result was observed for pdrg1 in T. dendrolimi. Moreover, this result also revealed that the apoptosis activity was enhanced although the diapause individual remained in a dormancy state. It is speculated that, during diapause, the wasps had accumulated harmful substances to some degree to maintain diapause.
Ribosome biogenesis has an vital role in the control of cell growth and division in eukaryotes . In the present study, we identified 31 DEGs involved in ribosome biogenesis, of which 28 were up-regulated during diapause in the prepupae stage (Fig. 6). It was reported that the rate of ribosome synthesis during diapause was lower than that of non-diapause eggs in Bombyx mori, so the upregulation of ribosomal proteins played an important role in blocking diapause [47, 48]. In mosquitos, Culex pipiens, the expression of rpS3a was dramatically reduced for a short time in the diapause stage. After the injection of rpS3a dsRNA into non-diapaused females, the development of follicle was arrested, similarly to in the diapause state . Conversely, in T. dendrolimi, diapause prepupae had a higher expression of ribosomal proteins compared to non-diapause wasps. It is noticeable that B. mori and C. pipiens enter diapause as adults, whereas T. dendrolimi enter diapause as prepupae. Although diapause decreased metabolism level, they may still need more energy to maintain this diapause condition.
UDP-glucuronosyltransferase (UDPGT) is of central importance in the union and next elimination of toxic xenobiotics and endogenous compounds. Glutathione-S-transferase (GST) belongs to a multifunctional protein family mainly located in the cytoplasm. Both are related to cellular detoxification. The transcriptional expression of GST and UDPGT were up-regulated during T. dendrolimi prepupal diapause. This was different from other species, such as solitary bee Tetrapedia diversipes  and Tetranychus urticae in which transcripts are down-regulated during diapause. In other species, the downregulation of GST and UDPGT might be correlated with non-feeding conditions, so the amount of exogenous substances decreases accordingly. However, T. dendrolimi is a kind of endoparasitoid wasp that spends its whole life cycle within the host egg, meaning starvation is a rare situation until the prepupae stage. Therefore, we speculated that the reason for the observed up-regulation of these two genes in diapause T. dendrolimi might be due to the increased resistance under unfavorable environmental conditions.
Cytochrome P450s (CYP450s) are hemoproteins involved in several physiological processes such as biosynthesis of hormones and degradation of xenobiotics . Previous studies showed that CYP450 of Schistosoma mansoni was essential for worm survival and egg development . Besides, CYP4G1 was proved to be related to cuticular hydrocarbon biosynthesis in Drosophila . In T. denrolimi transcriptome, we screened three genes (CYP4G15, CYP6A2, CYP6K1) belonging to two CYP families, CYP4 and CYP6. These genes were up-regulated in diapause individuals, suggesting that when T. denrolimi entered diapause, the environmental conditions were not suitable for survival. In this process, many harmful substances may be produced. So, a possible function for these genes is to balance out harmful substances, maintaining cellular homeostasis.
Autophagy is a conserved degradative pathway to remove unnecessary or dysfunctional cellular components [54, 55]. It is an indispensable catabolic and evolutionarily conserved process. Cytoplasmic components, including macromolecules and organelles are the targets in this pathway, which for degradation in response to stress conditions, such as developmental tissue remodeling, prolonged starvation periods, and nutritional fluctuations in the environment[56, 57]. To date, at least 30 autophagy-related genes (ATG) have been identified [58-60]. In T. dendrolimi, there were ATG1, ATG2, ATG9, and ATG13, which were significantly up-regulated in diapause.
In yeast, Saccharomyces cerevisiae, autophagy-related protein 13 (ATG13) is required for glycogen storage during the stationary phase . This result agreed with our results that the mRNA abundance of ATG13 was higher in diapause T. dendrolimi than non-diapause T. dendrolimi. In addition, for all we know, the present study is the first to report that ATG9 might be engaged in diapause. It has been reported that ATG1 and ATG13 might be regulated by the mTOR signaling pathway . In this study, the expressing profiles of genes in the mTOR signaling pathway were distinct among different diapause stages. Most of them were up-regulated, including ATG1 and ATG13. These findings were largely different from observations of other insect species (Fig. 7). The mTOR signaling pathway and ATGs may have an unknown regulatory relationship in terms of controlling diapause. A phylogenetic tree was constructed using ATG13 from different insect species (Fig. 8). TdATG13 did not cluster with species in the same genus but instead with species in Lepidoptera. This might be due to an evolutionary adaption between endoparasitoids and hosts.
Besides, ATG2 is the most downstream ATG protein in the preautophagosomal structure organization process. Valverde (2019) implied that lipid homeostasis mainly regulated by the activity of this autophagy-dependent protein , which was in line with Velikkakath’s previous finding (2012) . ATG2 also plays a significant role in life span extension . Certainly, extended longevity is one of the typical characteristics of diapause. Through phylogenetic analysis of ATG2, T. dendrolimi and N. vitripennis have high similarity, suggesting that ATG2 is possibly conservative to parasitoid wasp species (Fig. 9).
Lipid metabolism is essential for energy homeostasis of any organism. It has been shown that some diapausing insects use lipids as predominant energy stores [65, 66]. During diapause, almost all selected lipid metabolism-related genes were up-regulated, coinciding with the mobilization of TAG reserves (Fig. 10). These genes have important influence in the formation of triacylglycerol, which is the main caloric reserve during diapause. In addition, several genes’ expressions significantly increased after diapause termination. This might be due to post-diapause development. In our previous study, diapause T. dendrolimi had higher numbers of parasitized hosts than non-diapause T. dendrolimi . The increasing lipid storage can provide more energy for maintaining activities.
Trehalose, the main hemolymph sugar in many insects, functions in energy metabolism and protection in adverse ambient stresses, such as dehydration, heat, freezing, desiccation, and oxidation [67, 68]. During diapause, trehalose is accumulated to increase survival . Trehalase (tre) is an essential enzyme in trehalose catabolism. The trehalose transporter (tret) plays a decisive role in regulating trehalose concentration. In our results, tre was down-regulated during diapause of T. dendrolimi. However, the opposite was verified for tret, which was significantly up-regulated at the diapause stage (Fig. 11). The downregulation of tre and upregulation of tret increased the amount of trehalose during diapause.
Heat shock protein 68 (hsp68) was expressed constitutively in T. dendrolimi, although the expression level during the diapause stage was higher than that during the other three stages. HSPs are known as stress proteins, and Hsp68 belongs to the Hsp70 super family . As in T. dendrolimi, the hsp68 of Drosophila melanogaster are cold inducible. A study on D. melanogaster found that after inducing diapause by low temperature, hsp68 mRNA accumulated to higher amounts than in normal developmental conditions . According to this result, hsp68 might take effect in the process of T. dendrolimi diapause.