Analysis of the transcriptomes of Galeruca daurica(Coleoptera: Chrysomelidae) adults at different summer diapause stages

LI Yan-Yan, CHEN Long, LI Ling, TAN Yao, PANG Bao-Ping,*

(1. Research Center for Grassland Entomology, Inner Mongolia Agricultural University, Hohhot 010020, China;2. Department of Agronomy, Hetao College, Bayannoer, Inner Mongolia 015000, China)

Abstract: 【Aim】 This study aims to explore the crucial genes and metabolic pathways involved in obligatory summer diapause in Galeruca daurica, a new pest with great outbreak in the Inner Mongolia grasslands, northern China. 【Methods】 Using RNA-Seq, we performed sequencing, transcriptional profiling and functional prediction for G. daurica adults at different summer diapause stages, i.e. pre-diapause (PD), diapause (D) and post-diapause (TD), and screened differentially expressed genes (DEGs) at different summer diapause stages based on the RNA-Seq data. The expression levels of ten DEGs screened based on the RNA-Seq data were verified by qPCR. 【Results】 A total of 202 770 198 clean reads from nine libraries were filtered, and 12 078 060 transcripts were assembled into 82 292 unigenes with an average length of 783.59 bp and a N50 of 1 545 bp. The 2 395 (2 119 up-regulated and 277 down-regulated) and 62 (59 up-regulated and 3 down-regulated) DEGs were identified in the D vs PD and TD vs D comparison groups, respectively. The KEGG analysis revealed that the glycolysis/gluconeogenesis and fatty acid biosynthesis pathways were significantly enriched in the D vs PD and TD vs D comparison groups, respectively. Additionally, many DEGs related to the Ca2+ signaling were differentially expressed during diapause. Finally, the expressional analysis result by qPCR for ten DEGs showed a high consistency between qPCR and RNA-Seq results. 【Conclusion】 The glycolysis/gluconeogenesis, fatty acid biosynthesis, and Ca2+signaling pathways may play an important role in diapause regulation of G. daurica. This study establishes a foundation for future studies on the molecular mechanism underlying obligatory summer diapause in G. daurica.

Key words: Galeruca daurica; summer diapause; obligatory diapause; RNA-Seq

1 INTRODUCTION

The obligatory diapause is a fixed component that occurs at the particular stage of univoltine insect’s ontogenetic program, regardless of the ambient environment. On the other hand, facultative diapause is induced by environmental factors, usually found in multivoltine species, which can be “broken” under the favorite photoperiod or temperature conditions (Denlinger, 2002). According to the species, insects can enter diapause in any life cycle stage-embryonic, larval, pupal, or adult stages (Gilletal., 2017). Intriguingly, depending on the season when diapause happens, diapause can be divided into summer diapause and winter diapause. Much progress has not only been made in the physiological mechanisms of winter diapause (Hahn and Denlinger, 2011; Denlingeretal., 2012), but also some functional genes and categories involved in winter diapause regulation have been identified (Robichetal., 2007; Baker and Russell, 2009; Rinehartetal., 2010; Bryonetal., 2013; Qietal., 2015; Meyersetal., 2016). Summer diapause is frequently found in insects living in the hot and dry area, or triggered by food shortage, which facilitates insect’s survival in adverse environmental conditions during summer. Induction of summer diapause has been reported in insects among Lepidoptera (Xueetal., 1997; Liuetal., 2006), Orthoptera (Maigaetal., 2010; Kearneyetal., 2018), and Coleoptera (Xueetal., 2002; Ohashietal., 2003). Renetal. (2018) used the high-throughput RNA-Seq to identify candidate genes between summer diapause stage and non-diapause stage of the onion maggot,Deliaantiqua, and indicated that the circadian clock might play an important role in summer diapause induction. However, up to now, nearly all previous studies on the molecular mechanisms of diapause have focused on facultative diapause, and the molecular mechanism of obligatory diapause is nearly unknown. Therefore, an improved understanding on obligatory summer diapause is essential for understanding the special response of insects to hot and dry environment, and may provide better strategies for pest control.

Galerucadauricais an important grassland pest that is mainly distributed in China, Mongolia, Russia (Siberia) and Korea, and has recently extensively outbroken in Inner Mongolia, northern China since 2009 (Zhouetal., 2016).G.dauricahas typical diapause process including summer and winter diapause at different developmental stages in its life.G.dauricaoccurs one generation each year in Inner Mongolia, and the females deposit their eggs under cow dung, stone and in dense stands of grass in autumn. Overwintered eggs begin to hatch as early as in mid-April next year. The adults start to emerge in late-May. After feeding for about one week, adults gather under cow dung, stone and in dense stands of grass, and stop feeding and moving, and then enter obligatory summer diapause. Diapause terminates about three months later (Maetal., 2019). In previous studies, we found that the contents of lipid and total protein decreased gradually inG.dauricaadults during oversummering, whereas the water, total sugar and glycogen contents showed no obvious change (Chenetal., 2018a), the expression levels of diapause-related genesGdTRE1 andGdHsp10ain diapause adults were higher than those in pre-diapause adults (Chenetal., 2018b, 2019), and the geneJHBPwas lowly expressed in diapause adults (Chenetal., 2020). Both heat and cold stress could induce the expressions ofGdHsp20.6,GdHsp10 andGdHsp60 (Tanetal., 2018; Huoetal., 2019). However, the candidate genes and processes linked to summer diapause were unknown. In this study, we used RNA-Seq to identify the differentially expressed genes (DEGs) related to summer diapause inG.dauricaadults by comparing the pre-diapause, diapause and post-diapause groups. Additionally, 10 DEGs were selected to further test by qPCR to validate the transcriptome results. This transcriptome work provides an important molecular resource for obligatory summer diapause.

2 MATERIALS AND METHODS

2.1 Test insects

As described in Maetal. (2019),G.dauricaeggs were collected from the Sunite grasslands (44°21′N, 114°00′E) of Inner Mongolia, China, on April 17, 2014. Eggs were incubated (RH: 70%±5%; temperature: 25±1℃; photoperiod: 14L∶10D) until hatching. Larvae were transferred to natural conditions, and reared onAlliummongolium. Sampling was conducted at 3, 40, and 90 d after adult emergence to obtain adult samples at the pre-diapause (PD), diapause (D), and post-diapause (TD) stages, respectively. Thirty pairs of adults for each period were collected and stored at -80℃ for RNA extraction. Each time point had three biological replicates (10 pairs/replicate).

2.2 cDNA library construction and RNA sequencing

Samples of section 2.1 were ground in liquid nitrogen, and total RNA extraction was randomly selected from prepared samples using Trizol reagent following the manufacturer’s instructions (Invitrogen, Carlsbad, CA, USA) with three replicates. The purity, concentration, and integrity of the RNA samples were measured using the NanoDrop (Thermo Scientific, Massachusetts, USA), Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA), and Agilent 2100 (Agilent Technologies, Santa Clara, CA, USA), to ensure that transcriptome sequencing was performed using qualified samples. mRNA was enriched with Oligo (dT) magnetic beads and the mRNA was interrupted by fragmentation buffer. The first cDNA strand was synthesized by random hexamers using mRNA as a template. Buffer, dNTPs, RNase H and DNA polymerase I were added to synthesize the second cDNA strand. The cDNA was purified by AMPure XP beads. The purified double-stranded cDNA for end repair was A-tailed and ligated to the sequencing adapter. AMPure XP Beads were used for fragment size selection. Finally, the cDNA library was enriched by PCR amplification. The concentration of the library and insert size was detected by Qubit2.0 and Agilent 2100. The qPCR was utilized to quantify the effective concentration of the library. High-throughput sequencing was performed using Illumina Hiseq 2500 (Illumina, San Diego, CA, USA) in Biomaker Technologies Co, LTD (Beijing).

2.3 Sequence assembly and functional annotation

Raw data were filtered from the original sequencing results through in-house Perl scripts, removing the linker sequences containing poly-N and low-quality reads to obtain high-quality clean data.Denovoassembly of RNA-Seq was performed with Trinity (Grabherretal., 2011). The annotation of assembled sequences was performed by BLASTn and BLASTx (E-values≤1e-5) based on NCBI non-redundant protein/nucleotide sequences (Nr/Nt, http:∥www.ncbi.nlm.nih.gov/). Meanwhile, we also conducted annotation through Swiss-Prot (http:∥www.ebi.ac.uk/uniprot), Pfam (http:∥pfam.xfam.org/), GO (http:∥www.geneontology.org/), KOG/COG (https:∥www.ncbi.nlm.nih.gov/COG/) and KEGG (http:∥www.genome.jp/kegg/) (Tatusovetal., 2000; Kanehisaetal., 2004; Kooninetal., 2004; Sherlock, 2009).

2.4 Screening of differentially expressed genes (DEGs) based on the RNA-Seq data

Reads obtained from the sequencing of each sample were compared with the unigene library (assembled in section 2.3) using Bowtie (Langmeadetal., 2009). According to the comparison results, the expression level was estimated in combination with RSEM (RNA-Seq by Expectation-Maximization)(Li and Dewey, 2011). FPKM (fragments per kilobase of transcript per million mapped reads) value was used to represent the expression levels of the corresponding unigenes. The DEGs were analyzed using the DESeq2 v1.6.3 program package (http:∥www.bioconductor.org/packages/release/bioc/html/DESeq.html) (Anders and Huber, 2010). The DEGs were filtered with FDR (false discovery rate)<0.01 and fold change≥2. GO enrichment analysis of DEGs was implemented by the topGO R packages based on Kolmogorov-Smirnov test. Blast2GO program was used to assign GO annotations to three main categories (molecular function, cellular component and biological process) based on the Nr annotation (E-value cutoff≤1e-6). The unigene amino acid sequences were predicted with HMMER software and aligned using Pfam database to obtain unigene annotation information (Finn, 2005).

2.5 qPCR to validate DEGs screened based the RNA-Seq data

Ten DEGs including odorant-binding protein (OBP), 10 kDa heat shock protein (Hsp10), trehalase (Tre), enkurin domain-containing protein (EDP), phosphofructokinase (PFK), 2-oxoglutarate dehydrogenase (2-OGD), 3-phosphoinositide-dependent protein kinase (3-PDPK), sorbitol dehydrogenase (SDH), fatty acid synthase (FAS) and hypothetical protein YQE (HP) genes screened based on the RNA-Seq data were selected to validate the RNA-Seq results by qPCR. Samples were prepared for different summer diapause stages (PD, D and TD). Total RNA as extracted by Trizol reagent and the gDNA was removed using the Prime ScriptRT Reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan). Then, 1 μg of RNA was employed to synthesize the first-strand cDNA using the PrimeScriptTM1st Strand cDNA Synthesis Kit (TaKaRa, Shiga, Japan), according to the manufacturer’s instructions. Reaction system of 20 μL, including 10 μL of GoTaq?qPCR Master Mix, 2×(Promega, Wisconsin, USA), 0.8 μL of gene-specific primers (0.2 μmol/L), 7.2 μL Nuclease-Free Water and 2 μL Template cDNA, with the following program: 95℃ 10 min; followed by 40 cycles of 95℃ 15 s, 60℃ 60 s. Primers for qPCR were designed by using Primer Premier 5.0 (PREMIER Biosoft International, Palo Alto, California, USA) (Table 1). Succinate dehydrogenase (SDHA) gene ofG.dauricawas used as an internal gene (Tan Yetal., 2017). The gene-specific primers were designed using Primer3 (http:∥bioinfo.ut.ee/primer3-0.4.0/). Each treatment included three biological replicates, and each replicate had three technical replications. The relative expression levels were measured by the comparative 2-ΔΔCtmethod (Livak and Schmittgen, 2001). The significance of difference in the expression level of DEGs among different diapause stages were compared by analysis of variance (ANOVA), followed by multiple comparison (Turkey’s test,P<0.05) using the statistical software SPSS ver. 19.0 (SPSS Inc., Chicago, IL, USA).

Table 1 Primers used in qPCR

3 RESULTS

3.1 Sequencing, assembly and annotation

A total of 202 770 198 clean reads from nine libraries ofG.dauricaadults at different summer diapause stages (PD, D and TD) were filtered, with 51 084 099 992 clean nucleotides. After quality trimming and filtering, approximately 66.9 million, 65.6 million and 70.1 million clean reads were generated in the library of PD, D and TD, respectively. The final data set was used fordenovoassembly (Table 2) accessible at NCBI’s Sequence Read Archive (SRA) under accession number PRJNA471603, resulting in 12 078 060 transcripts with a N50 of 1 545 bp and an average length of 783.59 bp. The transcripts were further assembled into 82 292 unigenes, with a mean length of 783.57 bp, consisted of 65 262 unigenes (79.31% of the total unigenes) with the length of 200-1 000 bp, 9 615 (11.68%) with the length of 1 000-2 000 bp and 7 415 (9.01%)>2 000 bp. The annotation result showed that 36 127 unigenes (43.9% of the total unigenes) were annotated in databases Nr, Swiss-Prot, Pfam, GO, KEGG, COG, and KOG. In these unigenes, 35 442 (98.1%) unigenes had significant matches in Nr databases with the species distribution showing that 11 902 unigenes (34%) fromTriboliumcastaneumhad the greatest number of matches with those ofG.daurica, followed byDendroctonusponderosae(10%).

Table 2 Summary of the RNA-Seq data from the transcriptomesof diapause adults of Galeruca daurica

3.2 Functional analysis of DEGs

The 2 395 (2 119 up-regulated and 277 down-regulated) and 62 (59 up-regulated and 3 down-regulated) DEGs were identified in the DvsPD and TDvsD comparison groups, respectively (Fig. 1). The main enriched terms for the DEGs by GO classification were listed in Table 3. In the DvsPD comparison group, four terms were significantly enriched for the biological process including protein phosphorylation, microtubule-based movement, sperm axoneme assembly, and glycolysis, four terms were significantly enriched for the molecular function including protein serine/threonine kinase activity, microtubule motor activity, ATP binding, and protein kinase activity, and four terms were significantly enriched for the cellular component including dynein complex, cytoskeletal part, microtubule cytoskeleton, and centriole. However, in the TDvsD comparison group, five terms were significantly enriched for the biological process including positive regulation of biosynthetic process, polysaccharide catabolic process, lipid biosynthetic process, regulation of lipid metabolic process, and cellular lipid metabolic process, only two terms were significantly enriched for the molecular function including oxidoreductase activity and cellulase activity, and no term was significantly enriched for the cellular component.

A total of 593 DEGs in the DvsPD comparison group were divided into 129 known pathways. Only one KEGG pathway glycolysis/gluconeogenesis was significantly enriched (Table 4). Only two DEGs were grouped into two known pathways in the TDvsD comparison group, one of which fatty acid biosynthesis was significantly enriched (Table 5).

Table 3 Main enriched GO terms of differentially expressed genes (DEGs) in the D vs PD and TD vsD comparison groups of Galeruca daurica

Table 4 Main enriched KEGG pathways of differentially expressed genes (DEGs) in the D vs PDand TD vs D comparison groups of Galeruca daurica

Table 5 Differentially expressed genes (DEGs) significantly enriched in glycolysis/gluconeogenesis(D vs PD) and fatty acid biosynthesis (TD vs D) in the transcriptomes ofdiapause adults of Galeruca daurica

3.3 Validation of DEGs screened based on RNA-Seq data by qPCR

OBP,SDH,EDP,Hsp10,PFK,2-OGD,Treand3-PDPKgenes were mainly expressed during diapause. TheFASandHPgenes were highly expressed at the pre-diapause stage. The RNA-Seq and qPCR results showed a high consistency, and both can reflect the gene expression differences at three different diapause stages (Fig. 2).

Fig. 1 Summary of differentially expressed genes (DEGs)in the transcriptomes of diapause adults of Galerucadaurica identified through pairwise comparisonsThe overlapping areas represent the genes we filtered by comparing two or three groups. D: Diapause; PD: Pre-diapause; TD: Post-diapause. The same below.

Fig. 2 Validation of the RNA-Seq results of 10 differentially expressed genes (DEGs) in the transcriptomesof diapause adults of Galeruca daurica by qPCRA: OBP; B: Tre; C: FAS; D: EDP; E: Hsp10; F: PFK; G: 2-OGD; H: HP; I: 3-PDPK; J: SDH. The qPCR results are shown with columns, and the relative expression levels are shown on the left y-axis; the RNA-Seq results are shown with lines, and the FPKM values correspond to the right y-axis.

4 DISCUSSION

The development of diapausing insect was characterized by extremely low metabolic activity (Ishikawaetal., 2000; Zhangetal., 2017). Glycolysis is a metabolic pathway that converts glucose into pyruvate and H+, and the free energy released in this process is used to form the ATP and NADH. Gluconeogenesis is a metabolic pathway that results in the generation of glucose from non-carbohydrate carbon substrates such as pyruvate, lactate, glycerol, glucogenic amino acids, and fatty acids. In this study, glycolysis/gluconeogenesis was only one significantly enriched pathway after diapause entry (DvsPD) (Table 4), and all 16 DEGs enriched in this pathway were up-regulated during diapause, indicating that this metabolic pathway provides the main source of energy required for maintaining life activities during diapause. Other researches have also shown that the glycolysis/gluconeogenesis pathway plays important roles in insect diapause (Qietal., 2015; Zhangetal., 2018).

The tricarboxylic acid cycle (TCA cycle) is a series of enzyme-catalyzed chemical reactions that form a key part of aerobic respiration in cells. The TCA cycle may be a checkpoint for developmental arrest in diverse animal groups (Xuetal., 2012). In this study, 11 DEGs were enriched in the TCA cycle during diapause. Other researchers also recently found that many genes involved in the TCA cycle were differentially expressed during diapause (Renetal., 2018; Zhaietal., 2019). These results suggested that the TCA cycle might play an important role in regulating insect diapause.

Trehalose, as blood glucose, was used for energy utilization and protection in unfavorable environment. In embryos (cysts) ofArtemia, a large amount of trehalose was stored during the diapause process not only for surviving, but also for providing energy after diapause terminated (Yangetal., 2013). Trehalase, a catalytic trehalose hydrolysis enzyme, has a special effect on the metabolism of trehalose. In this study, theTreexpression level was much higher during diapause (Fig. 2: B). Qietal. (2015) obtained the similar result inC.septempunctata. These results demonstrated that trehalase could balance the trehalose content and provide energy supply. We further found that the expression level of a sorbitol dehydrogenase gene (SDH) during diapause stage was ～50 times more than that of before diapause, and ～2 times more than that of post diapause (Fig. 2: J). InWyeomyiasmithii, theSDHexpression level was high during larval diapause, but after diapause terminated, it decreased gradually (Emersonetal., 2010). However, theSDHexpression level ofBombyxmoriduring egg diapause was low, and other studies suggested that it was related to the termination of diapause (Yaginumaetal., 1990). TheSDHexpression level was four times higher at the diapause stage ofC.septempunctatathan at other two stages, and different expression scales existed in different diapause types (Qietal., 2015). From the above analysis, theSDHexpression has the same tendency for diapause at the same developmental stage, but it is different when the developmental stages for diapause varies. We deduced that it may be because the SDH has different functions underlying diapause at different developmental stages.

Fatty acid metabolism was a significantly enriched pathway during diapause process of some insects (Qietal., 2015). Storing energy materials (i.e., sugar, fat) before diapause is essential. Lipids generate high caloric content and water, and store energy more advantage than carbohydrate-based sources (MacRae, 2010).Leptinotarsadecemlineatadoes not feed during diapause, and accumulates much lipid after emergence (Rinehartetal., 2007). In this study, the expression level of a fatty acid synthase gene (FAS) that is associated with the accumulation of lipids was much higher at the pre-diapause stage than during diapause (Fig. 2: C), and fatty acid biosynthesis was also only one significantly enriched pathway for the sole DEG (FAS) after diapause terminated (Table 5). The results indicated thatG.dauricamight accumulate a great deal of lipids at the pre-diapause stage as the main energetic materials during diapause. Several other insect species such asEctomyeloisceratoniae(Heydari and Izadi, 2014) andColaphellusbowringi(Tan QQetal., 2017) also accumulate lipids before diapause as energy resource in diapause. However, from the early to late diapause of the mosquitoCulexpipiens, a large number of genes associated with lipolysis were up-regulated (Denlinger, 2002). Two genes involved in fatty acid synthesis and modification were up-regulated in diapausingAedesalbopictuslarvae (Michaudetal., 2002).

Phosphatidylinositol signaling system plays an important role in mediating numerous physiological processes in eukaryotic organisms, such as growth, cytoskeletal rearrangement and membrane trafficking (Zhouetal., 2014). In this study, twelve DEGs were enriched in this pathway during diapause (Table 5), including five unigenes (c44629.graph_c0, c42861.graph_c0, c41117.graph_c0, c56932.graph_c0, and c42861.graph_c1) encoding calmodulin also involved in the Ca2+signaling pathway, which were up-regulated during diapause. All five calmodulin in this study were also up-regulated during diapause whereas one calmodulin was down-regulated at the protein level in our proteomic study (Maetal., 2019). However, Zhaoetal. (2017) found that the abundance of calmodulin declined at both the transcript and protein levels in the diapause stage ofTetranychusurticae. Calmodulin was also down-regulated in larval diapause initiation ofHelicoverpaarmigera(Zhangetal., 2012) andD.antiqua(Renetal., 2018). Therefore, it is necessary to further study whether the phosphatidylinositol signaling system functions in the regulation of insect diapause.