Introduction
The oriental armyworm, Mythimna separata (Walker) (Lepidoptera: Noctuidae), is typical of widespread migratory pests, and is widely distributed in Asia, Europe, and Oceania (Sharma et al., Reference Sharma, Sullivan and Bhatnagar2002; Wang et al., Reference Wang, Zhang, Ye and Luo2006; Li et al., Reference Li, Li, Jiang, Zhao, Xu and Wu2019a). The agricultural economy suffers giant losses due to damage by M. separate of several crops including maize, wheat, and rice (Sharma et al., Reference Sharma, Sullivan and Bhatnagar2002; Li et al., Reference Li, Li, Wang, He, Tan and Yang2019b). In China, there were about 12 M. separata outbreaks in the period 1970–1989 (Ye, Reference Ye1993); in 2012, the actual losses of grain production caused by M. separata were 9.92 × 105 tons, and in 2013, 3.93 × 105 tons (Jiang et al., Reference Jiang, Li, Zeng and Liu2014). Effective control of M. separata is essential for the successful harvest of grain (Sharma et al., Reference Sharma, Sullivan and Bhatnagar2002; Jiang et al., Reference Jiang, Li, Zeng and Liu2014; Li et al., Reference Li, Li, Jiang, Zhao, Xu and Wu2019a). At present, the main methods for the control of armyworms are ecological regulation and integrated use of physical, biological, and chemical measures (Jiang et al., Reference Jiang, Li, Zeng and Liu2014). Although these methods are to a degree effective in the control of armyworms, there are significant associated challenges such as insecticide resistance (Li et al., Reference Li, Li, Wang, He, Tan and Yang2019b).
Mythimna separata is a holometabolous insect, and distinguished by the four developmental stages of egg, larva, pupa, and adult (Truman and Riddiford, Reference Truman and Riddiford1999). The pupal stage is characteristic in Lepidoptera and occurs between vegetative (larval stage) and reproductive growth (adult stage) (Barbora et al., Reference Barbora, Vlastimil and Marek2011). After M. separata larva mature, they drill 1–2 cm into the soil near the roots of host plants and construct a soil chamber for pupation. This transformation is characterized by dramatic morphological and structural changes and is an extremely complicated biological process regulated by an array of genes, hormones, and nutritional components (Dubrovsky, Reference Dubrovsky2004). Manipulation of insect metamorphosis by regulating genes involved in the deformation process has the potential in new pest control strategies.
The epidermis of insects not only protects against pathogens and harmful aspects of the environment but also shapes and maintains activity abilities during normal development (Delon and Payre, Reference Delon and Payre2004; Moussian et al., Reference Moussian, Schwarz, Bartoszewski and Nüsslein-Volhard2005). Cuticular proteins (CPs) are structural proteins of the insect epidermis. The number of CPs is usually more than 1% of the total insect coding proteins in the genome (Futahashi et al., Reference Futahashi, Okamoto, Kawasaki, Zhong, Iwanaga, Mita and Fujiwara2008), which indicated that CPs may play an important role in growth and development, environmental adaptation, and innate immunity (Andersen et al., Reference Andersen, Hojrup and Roepstorff1995; Willis, Reference Willis2010). The expression of CP genes is closely related to the rhythms of insect molt and metamorphosis (Andersen, Reference Andersen2010). Particularly, morphological characteristics change markedly during prepupal–pupal metamorphosis (Fraenkel and Rudall, Reference Fraenkel and Rudall1940; Gu et al., Reference Gu, Huang, Gong, Zheng, Liu, Huang and Feng2013). Guan et al. (Reference Guan, Middlebrooks, Alexander and Wasserman2006) found that deletion of the epidermal protein gene (Tweedle D1) shortens the body of larvae and pupae of Drosophila melanogaster, which implicates epidermal protein genes in body formation. In Anopheles gambiae, five genes of the CPs family were found to be expressed before the molting of pupae or adults, suggesting that they may be involved in the formation of the outer epidermis of these stages (Togawa et al., Reference Togawa, Dunn, Emmons and Willis2007). Therefore, CPs are the prime candidates for inclusion in the mechanistic models of insect molting and metamorphosis, and a key to improving the understanding of the biological and physical chemistry, and structural modification that occurs during the development of insects (Shahin et al., Reference Shahin, Iwanaga and Kawasaki2016).
Juvenile hormone (JH) and 20-hydroxyecdysone (20E) are two major hormones in insects and coordinate in the orchestration of insect growth and development, including features of growth, molting, and reproduction (Liu et al., Reference Liu, Dai, Guo, Li, Ma, Tian, Cao, Zhang, Palli and Li2015; Jia et al., Reference Jia, Liu, Wen, Cheng, William, Wang and Li2017). The two hormones appear to complement each other; JH maintains insect growth in larvae, whereas 20E induces metamorphosis of the mature larvae (Truman and Riddiford, Reference Truman and Riddiford2002; Jindra et al., Reference Jindra, Palli and Riddiford2013). JH is synthesized and secreted by insect corpora allata, which regulates the growth of larvae, promotes ovarian maturation, and prevents larvae from entering the next instar stage prematurely (Qu et al., Reference Qu, Bendena, Tobe and Hui2018). JH is also involved in the physiological and biochemical processes of the insect body, molting, diapause, and stimulation of the yolk in female adults (Sugahara et al., Reference Sugahara, Tanaka and Shiotsuki2017). 20E is synthesized and secreted by insect prothoracic gland cells. At low concentration, it mainly controls growth, while high concentration mainly regulates the molting between insect larvae and the metamorphosis during development (Nijhout and Callier, Reference Nijhout and Callier2015). Regulation of metamorphosis by hormones has been studied in many insects, including Bombyx mori (Zhang et al., Reference Zhang, Liu, Shiotsuki, Wang, Xu, Huang, Li, Li and Tan2017), Aedes aegypti (Liu et al., Reference Liu, Fu and Zhu2018), and Laodelphax striatellus (Zhai et al., Reference Zhai, Zhang, Gao, Chen, Sun, Zhang, Yu and Zheng2017). From these studies, it appears that JH and 20E co-mediate the prepupal–pupal metamorphosis process of insects by regulating the expression patterns of CPs in various tissues and at various developmental stages.
Previous studies on metamorphosis have been mainly focused on the gene expression patterns of the cytochrome P450s, very high-density lipoproteins, chitinase, chitin deacetylase, serine protease, CPs, and 20E, and which have been studied in insects including Spodoptera litura (Gu et al., Reference Gu, Huang, Gong, Zheng, Liu, Huang and Feng2013), Bactrocera dorsalis (Chen et al., Reference Chen, Hou, Dou, Wei, Yue, Yang, Yu, De Schutter, Smagghe and Wang2018), B. mori (Zhang and Zheng, Reference Zhang and Zheng2017), and Helicoverpa armigera (Zhang et al., Reference Zhang, Huang, Zhou, Zhang, Liu, Miao and Huang2009; Cai et al., Reference Cai, Zhao, Jing, Song, Zhang, Wang and Zhao2016). However, little is known about the molecular mechanisms of metamorphosis in M. separata. In this study, we elucidate the larval prepupal–pupal transition by performing a comparative transcriptome analysis. Using RNA-seq, we discovered genes that were differentially expressed at different stages of development, results which were confirmed by real-time quantitative PCR (RT-qPCR). The results of this study provide a foundation for further understanding the molecular mechanism of metamorphosis in M. separata, and also open an avenue for establishing RNAi or CRISPR-based pest control through regulating the development of this pest.
Materials and methods
Insects and sample preparation
For these studies, we used a susceptible strain of M. separata that has been raised in the laboratory for more than 50 generations without exposure to any insecticides. The feeding method used is described in Li et al. (Reference Li, Li, Wang, He, Tan and Yang2019b). Mature larvae (ML), wandering stage (W), and three stages of pupation (1, 5, and 10 days after pupation, designated P1, P5, and P10, respectively) were collected (fig. 1). The ML stage of M. separata is black or gray black and have the longest body length about 40–50 mm, while the body length of W stage of M. separata is shorter than ML and often bent. The P1 are yellow after that the pupae change from yellow to brown. Three biological replicates were conducted, each with three individuals sampled for each developmental stage.
Figure 1. The morphology of five developmental stages including mature larvae (ML), wandering (W), and different stages of pupation (1, 5, and 10 days after pupation were designated P1, P5, and P10) during the pupariation of M. separata. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
RNA isolation, library construction, and sequencing
We used the Eastep® Super Total RNA Extraction Kit (Promega, Shanghai, China) to extract the total RNA of the whole body of M. separata samples, according to the manufacturer's instruction. RNA degradation and contamination were monitored on 1% agarose gels. The concentration and quality of RNA were assessed using a spectrophotometer NanoDrop 2000c (Thermo Fisher Scientific, Waltham, MA, USA).
After total RNA was treated with DNase I, mRNA was enriched with magnetic beads with Oligo (dT). Next, a disruption reagent was added to break the mRNA into short fragments, and we used the disrupted mRNA as a template to synthesize the first-strand cDNA with random hexamers. Next, buffer, dNTPs, and DNA polymerase I were added to synthesize the second-strand cDNA. After purification and recovery, repair of sticky ends, the addition of a base A at the 3′ end, and ligation of a sequencing adapter, the resulting fragments were selected for size and enriched by PCR. The constructed library was qualified by the Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System and was sequenced using the Illumina sequencing platform at Shenzhen Hengchuang Gene Technology Co., Ltd. (TGS, Shenzhen, China).
Raw data quality control
The Illumina platform converts sequenced image signals into text signals and stores them in FASTQ format as raw data (Cock et al., Reference Cock, Rice, Goto, Heuer and Fields2010). The original sequences were filtered to eliminate low-quality reads. Generally, the accuracy of sequencing decreases with the read length of reads, and reads are considered reliable when the Q20 of Reads is above the set threshold; Q20 indicates that the sequencing error rate is 1%, and Q20 reflects the quality of sequencing. Reads were filtered under the threshold of containing more than 5% N bases, or 50% bases with mass values less than 10. The filtered reads, or clean reads, were used in subsequent analysis.
De novo assembly
The Trinity software (Grabherr et al., Reference Grabherr, Haas, Yassour, Levin, Thompson, Amit, Adiconis, Fan, Raychowdhury, Zeng, Chen, Mauceli, Hacohen, Gnirke, Rhind, Palma, Birren, Nusbaum, Lindblad-Toh, Friedman and Regev2011) was used to assemble clean reads. Sequences were de-replicated and clustered using Tgicl (Pertea et al., Reference Pertea, Huang, Liang, Antonescu, Sultana, Karamycheva, Lee, White, Cheung, Parvizi, Tsai and Quackenbush2003). Unigenes obtained by Tgicl dereplication and splicing are divided into clusters and singletons. In a given cluster, there can be several unigenes (starting with CL, followed by CL family number) with a high degree of similarity (>70%). Others are singletons (starting with unigene), which are individual unigenes.
Functional annotation of unigene set
The obtained unigene set was annotated by Blast (Altschul et al., Reference Altschul, Gish, Miller, Myers and Lipman1990) search against several databases: non-redundant nucleic acid sequence (NT), non-redundant protein sequences (NR), Gene Ontology (GO), Clusters of Orthologous Groups of proteins (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and manual annotation and expected non-redundant protein database (SwissProt). Blast2GO (Conesa et al., Reference Conesa, Gotz, Garcia-Gomez, Terol, Talon and Robles2005) and NR annotation results were used to annotate the unigene set.
Prediction of the coding sequences of unigenes
According to the functional annotation results, we selected the best-aligned segment as the coding sequences (CDS) of the unigene. The coding sequence was translated into amino acid sequences according to the standard codon table. Unigenes of the un-annotated model were predicted using the CDS predicted in the previous step along with ESTScan (Iseli et al., Reference Iseli, Jongeneel and Bucher1999).
Differential expression analysis of unigenes
RNA-Seq by Expectation-Maximization (Dewey and Li, Reference Dewey and Li2011) was used to calculate the abundance of each transcript for each library. The expression of unigene was standardized by the Fragment Per Kilobase of exon model per Million mapped reads (FPKM) method. Differentially expressed genes (DEGs) between different developmental stages were calculated using the DESeq2 software (Love et al., Reference Love, Huber and Anders2014). The normalization factors were calculated using the trimmed mean of M-values (TMM) method. The threshold false discovery rate (FDR) < 0.05 was adjusted to identify the DEGs by fold change (≥2). Significance was adjusted with the threshold FDR < 0.05 and log2fold change (|Log2FC|) > 1. Heat maps were drawn in the R programming language. A detailed description of the calculation formula of FPKM and differential expression analysis of unigene is given in the Supplementary Information.
Phylogenetic analysis of DEGs
We constructed a phylogenetic tree to analyze the sequence homology and evolutionary relationships of DEGs of the CPs, 20E, and JH biosynthesis and signaling pathway, including protein sequences from closely related insects S. litura, Athetis dissimilis, Trichoplusia ni, Spodoptera frugiperda, Galleria mellonella, B. mori, H. armigera, and Spodoptera exigua, which were obtained from the NCBI repository. Sequences were aligned using Clustalx 1.81. A phylogenetic tree was constructed by MEGA5 with the neighbor-joining method and 1000 bootstrap replicates.
Quantitative real-time PCR
Ten DEGs of the CP genes, and genes in the 20E and JH biosynthesis and signaling pathway, were selected and further verified using quantitative real-time PCR (qRT-PCR). Total RNA of samples from ML, W, P1, P5, and P10 were extracted using RNAiso Plus (Takara, Dalian, China) according to the manufacturer's instructions. The concentration and quality of the extracted RNA were then detected with a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific). First-strand cDNA was synthesized from 1 μg total RNA using a PrimeScript™ RT reagent Kit with gDNA Eraser (Takara) according to the manufacturer's instructions. Primer pairs (table S1) were designed using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/). The β-Actin (Accession number: GQ856238) and GAPDH of M. separata were used as reference genes (Li et al., Reference Li, Dai, Zhang, He, Ran and Chen2018; Li et al., Reference Li, Li, Wang, He, Tan and Yang2019b). qRT-PCR was performed on a Bio-Rad CFX96 (BioRad, Hercules, USA) thermal cycler. The relative expression of candidate genes was calculated using the 2−ΔΔCT method (Livak and Schmittgen, Reference Livak and Schmittgen2001). A detailed description of the RT-qPCR is given in the Supplementary Information.
Data analysis
All experiments were performed with three biological replications, and data were analyzed with SPSS Statistics 20 (IBM, Chicago, IL, USA). The differences in mRNA expression in different developmental stages were assessed by one-way analysis of variance and the Duncan-test (P ≤ 0.05).
Results
Transcriptome sequencing and assembly
Illumina transcriptomes of total RNA were sequenced and assembled for five different developmental stages of M. separata. The RNA-seq raw sequence data have been submitted to SRA at NCBI with SRA ID of PRJNA604659 and the accession number SRR11011830. The assembly results totaled 140,562 contigs, with an average length of 1201.07 bp, N50 of 2159 bp, and GC of 41.94%. A detailed description of the assembly results is shown in table 1.
Table 1. Statistics derived from the transcriptomes of the prepupal–pupal transition of M. separata.
ML, mature larvae; W, wandering; P1, P5, and P10: 1, 5, and 10 days after pupation. ML_1, ML_2, and ML_3 represent three replicates.
Through Blast search of NCBI (e value < e–3), more than half of genes were annotated (74,059) according to the target sequences. Among them, 6120, 351,851, 47,632, 48,973, 27,817, and 31,569 unigenes were annotated in NR, NT, Swiss-Prot, KEGG, COG, and GO, respectively (fig. S1), more than in B. dorsalis (Chen et al., Reference Chen, Hou, Dou, Wei, Yue, Yang, Yu, De Schutter, Smagghe and Wang2018) and Chilo suppressalis (Sun, Reference Sun2015). A total of 19,733 unigenes were collectively annotated over the five databases (fig. S2).
Species homology analysis gave the most paired species as Amyelois transitella 12,958 (21.17%), B. mori 10,635 (17.38%), Papilio xuthus 5457 (8.92%), Papilio machaon 4205 (6.87%), Danaus plexippus 3085 (5.04%), Operophtera brumata 2920 (4.77%), Papilio polytes 2348 (3.84%), Plutella xylostella 2335 (3.82%), H. armigera 1148 (1.88%). Other insect species were hit at <1% (fig. S3).
Global analysis of differentially expressed unigenes
There were 22,884 (10,990 up-regulated and 11,894 down-regulated), 23,534 (11,534 up-regulated and 12,000 down-regulated), 26,643 (12,718 up-regulated and 13,925 down-regulated), and 33,238 (18,001 up-regulated and 15,237 down-regulated) DEGs identified between ML and W, W and P1, P1 and P5, and P5 and P10, respectively (fig. 2a). The number of DEGs between other groups is shown in fig. S4. A total of 3135 (0.76%) DEGs were found to be common over the five developmental stages (fig. 2b).
Figure 2. Histogram (a) and Venn diagram (b) analysis of the number of differentially expressed genes (DEGs) between samples. ML, mature larvae; W, wandering; P1, P5, and P10: 1, 5, and 10 days after pupation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
GO distribution analysis of DEGs of the prepupal–pupal transition
GO distribution analysis results showed that these DEGs were categorized into 59 sub-categories belonging to the three main GO categories (biological processes, cellular components, and molecular functions). These categories consisted of 17, 22, and 21 sub-categories respectively (fig. 3).
Figure 3. Gene ontology (GO) enrichment analysis of differentially expressed genes (DEGs). (a–d) Represent the GO term for the ML-vs.-W, W-vs.-PP, P1-vs.-P5, and P5-vs.-P10, respectively. Unigene statistical map of ‘up-down’ differential expression of enriched GO term. The x-axis represents the GO term, and the y-axis represents the number of up-down genes corresponding to the GO term. A1–A30 represent the carbohydrate metabolic process, carbohydrate transport, chitin metabolic process, defense response, immune system process, innate immune response, metabolic process, oxidation-reduction process, proteolysis, transmembrane transport, chorion, collagen, extracellular region, extracellular space, integral component of membrane, membrane, mitochondrial inner membrane, mitochondrial membrane, mitochondrial proton-transporting ATP synthase complex, coupling factor F(o), troponin complex, chitin binding, heme binding, hydrolase activity, iron ion binding, monooxygenase activity, oxidoreductase activity, oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular, peptidase activity, serine-type endopeptidase activity, serine-type peptidase activity GO term, respectively. B1–B30 represent the carbohydrate metabolic process, carbohydrate transport, chitin metabolic process, defense response to bacterium, immune system process, innate immune response, metabolic process, oxidation–reduction process, proteolysis, transmembrane transport, chorion, collagen, cytoskeleton, extracellular region, extracellular space, integral component of membrane, membrane, mitochondrial membrane, mitochondrial proton-transporting ATP synthase complex, coupling factor F(o), peroxisome, carboxylic ester hydrolase activity, chitin binding, heme binding, hydrolase activity, monooxygenase activity, oxidoreductase activity, peptidase activity, serine-type endopeptidase activity, serine-type peptidase activity, structural constituent of cuticle GO term, respectively. C1–C30 represent the carbohydrate metabolic process, chitin metabolic process, defense response, defense response to bacterium, lipid metabolic process, metabolic process, oxidation–reduction process, proteolysis, steroid hormone-mediated signaling pathway, transport, Golgi membrane, Golgi stack, LUBAC complex, chorion, extracellular region, extracellular space, integral component of membrane, membrane, postsynaptic membrane, troponin complex, carboxylic ester hydrolase activity, chitin binding, hydrolase activity, oxidoreductase activity, peptidase activity, sequence-specific DNA binding transcription factor activity, serine-type endopeptidase activity, serine-type peptidase activity, structural constituent of cuticle, transporter activity GO term, respectively. D1–D10 represent the carbohydrate metabolic process, glycolysis, innate immune response, ion transport, metabolic process, oxidation–reduction process, proteolysis, translation, transmembrane transport, transport, cytoplasm, cytosolic large ribosomal subunit, extracellular region, integral component of membrane, intracellular, membrane, mitochondrial inner membrane, mitochondrion, ribonucleoprotein complex, ribosome, calcium ion binding, catalytic activity, hydrolase activity, metal ion binding, nucleotide binding, oxidoreductase activity, peptidase activity, structural constituent of cuticle, structural constituent of ribosome, transferase activity GO term, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
For biological processes, the oxidation–reduction process was the most significant prepupal–pupal transition-related GO term in ML-vs.-W and P5-vs.-P10, with 700 (6.06%, 193 up-regulated and 507 down-regulated genes) and 1023 (4.77%, 585 up-regulated and 438 down-regulated genes) DEGs, respectively, while the most significant prepupal–pupal transition-related GO term was found in the proteolysis in W-vs.-P1, P1-vs.-P5, with 629 (5.95%, 233 up-regulated and 396 down-regulated genes) and 622 (5.83%, 215 up-regulated and 407 down-regulated genes) DEGs, respectively (fig. 3).
GO enrichment analyses of DEGs showed that the membrane and integral component of the membrane in the cellular component category were significantly enriched. A total of 2388 (8.11%, 1082 up-regulated and 1306 down-regulated genes), 2245 (9.49%, 901 up-regulated and 1344 down-regulated genes), 2404 (9.41%, 1115 up-regulated and 1289 down-regulated genes), and 3615 (8.03%, 2024 up-regulated and 1591 down-regulated genes) DEGs were enriched in the membrane sub-category in ML-vs.-W, W-vs.-P1, P1-vs.-P5, and P5-vs.-P10, respectively. Integral component of membrane was enriched with 2346 (7.97%, 1064 up-regulated and 1282 down-regulated genes), 2208 (9.33%, 892 up-regulated and 1316 down-regulated genes), 2372 (9.28%, 1093 up-regulated and 1279 down-regulated genes), and 3514 (7.80%, 1975 up-regulated and 1539 down-regulated genes) DEGs in ML-vs.-W, W-vs.-P1, P1-vs.-P5, and P5-vs.-P10, respectively (fig. 3). These DEGs were mainly associated with the solute carrier family, MFS transporter, low-density lipoprotein receptor, glucuronosyltransferase, ATP-binding cassette (ABC), P450, alcohol-forming fatty acyl-CoA reductase, ATPase, NADH dehydrogenase, stearoyl-CoA desaturase, chitin, fatty acids, and E3 ubiquitin-protein (table S3). Also, the CPs, 20E, and JH-related genes, which are involved in the metamorphic development of Lepidoptera insects (Willis, Reference Willis2010; Liu et al., Reference Liu, Dai, Guo, Li, Ma, Tian, Cao, Zhang, Palli and Li2015), were also identified in the membrane and integral component of membrane sub-categories (table S4), suggesting these genes are likely to play a role in the prepupal–pupal transition of M. separata.
Hydrolase activity was the most significantly affected GO term under molecular function in ML-vs.-W, W-vs.-P1, P1-vs.-P5, with 784 (6.79%, 286 up-regulated and 498 down-regulated genes) and 727 (9.49%, 201 up-regulated and 526 down-regulated genes), 689 (6.46%, 271 up-regulated and 418 down-regulated genes) DEGs, respectively. Metal ion binding was significantly enriched in P5-vs.-P10, with 1262 (5.88%, 596 up-regulated, and 666 down-regulated genes) (fig. 3).
KEGG pathway analysis of DEGs potentially involved in the prepupal–pupal transition
KEGG pathway analysis indicated that the DEGs were enriched in the metabolic pathways during the prepupal–pupal transition of M. separata. These DEGs were mainly associated with starch and sucrose metabolism, pyrimidine metabolism, protein digestion, and absorption (fig. 4). In the metabolic pathways, there were 2218 (16.85%, 745 up-regulated and 1473 down-regulated genes), 2053 (16.26%, 890 up-regulated and 1163 down-regulated genes), 2046 (14.64%, 914 up-regulated and 1132 down-regulated genes), and 3037 (16.14%, 1890 up-regulated and 1147 down-regulated genes) DEGs in ML-vs.-W, W-vs.-P1, P1-vs.-P5, and P5-vs.-P10, respectively.
Figure 4. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of specific differentially expressed genes (DEGs) in different stages of M. separata. (a–d) Represent the KEGG pathway for the ML-vs.-W, W-vs.-PP, P1-vs.-P5, and P5-vs.-P10, respectively. RichFactor indicates the percentage of the ratio of genes in each pathway vs. the total genes in the term. A1–A20 represent steroid hormone biosynthesis, starch and sucrose metabolism, retinol metabolism, protein digestion and absorption, phenylalanine metabolism, pentose and glucuronate interconversions, metabolism of xenobiotics by cytochrome P450, metabolic pathways, malaria, lysosome, glycosaminoglycan degradation, galactose metabolism, ECM–receptor interaction, drug metabolism-other enzymes, drug metabolism-cytochrome P450, chemical carcinogenesis, bile secretion, ascorbate and aldarate metabolism, amoebiasis, amino sugar and nucleotide sugar metabolism KEGG pathway, respectively. B1–B20 represent tyrosine metabolism, riboflavin metabolism, renin-angiotensin system, protein digestion and absorption, phenylalanine metabolism, phagosome, peroxisome, other glycan degradation, neuroactive ligand–receptor interaction, metabolic pathways, malaria, lysosome, glycosphingolipid biosynthesis-ganglio series, glycosaminoglycan degradation, fatty acid elongation, fat digestion and absorption, ECM–receptor interaction, complement and coagulation cascades, amoebiasis, amino sugar and nucleotide sugar metabolism KEGG pathway, respectively. C1–C20 represent steroid hormone biosynthesis, retinol metabolism, RNA polymerase, pyrimidine metabolism, protein digestion and absorption, prion diseases, primary immunodeficiency, pentose and glucuronate interconversions, herpes simplex infection, glycosaminoglycan degradation, galactose metabolism, fatty acid metabolism, fatty acid elongation, fat digestion and absorption, ECM–receptor interaction, dorso-ventral axis formation, complement and coagulation cascades, biosynthesis of unsaturated fatty acids, amoebiasis, ABC transporters KEGG pathway, respectively. D1–D20 represent Vibrio cholerae infection, ribosome, renin secretion, RNA polymerase, pyruvate metabolism, phagosome, pathogenic Escherichia coli infection, Parkinson's disease, oxidative phosphorylation, non-alcoholic fatty liver disease (NAFLD), metabolic pathways, hypertrophic cardiomyopathy (HCM), Huntington's disease, glyoxylate and dicarboxylate metabolism, ECM–receptor interaction, dilated cardiomyopathy, cardiac muscle contraction, carbon metabolism, amoebiasis, Alzheimer's disease KEGG pathway, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Besides, there were only 35 (0.27%), 26 (0.21%), 28 (0.20%), and 37 (0.2%) DEGs enriched in insect hormone biosynthesis, which was closely associated with JH and 20E, and 78 (0.59%), 73 (0.58%), 79 (0.57%), 143 (0.76%) and 60 (0.46%), 50 (0.40%), 70 (0.50%), and 73 (0.39%) DEGs enriched in CPs belongs to the GnRH signaling pathway and the ErbB signaling pathway in ML-vs.-W, W-vs.-P1, P1-vs.-P5, and P5-vs.-P10, respectively (fig. S5).
We further investigated the gene families which were of particular interest as potential targets for the prepupal–pupal transition of M. separate below.
The number of DEGs detected in the CPs, 20E, and JH biosynthesis
A total of 30 genes encoding CPs, and 17 and seven genes involved in 20E and JH biosynthesis, respectively, were detected during the prepupal–pupal transition of M. separata (table 2). The number of genes in CPs, 20E, and JH in the five different development stages of M. separata is similar (table S5). We tested whether those genes were clustered consistently in a phylogenetic context. Results showed that CP (fig. S6), 20E (fig. S7), and JH (fig. S8) genes from M. separata were well clustered into their relevant phylogenetic branches based on the phylogenetic analysis with eight lepidopteran insects including S. litura, A. dissimilis, T.ni, S. frugiperda, G. mellonella, B. mori, H. armigera, and S. exigua.
Table 2. Differentially expressed genes in the prepupal–pupal transition of M. separata
ML, mature larvae; W, wandering; P1, P5, and P10: 1, 5, and 10 days after pupation.
The second to fifth columns of numbers represent |log2Ratio|.
+ is up-regulated, − is down-regulated and – represents no DEGs. Differential gene screening conditions are FDR ≤ 0.01 and |log2Ratio|≥1.
Expression levels of CP genes
There were 11 up-regulated and two down-regulated CP genes with a maximum |log2Ratio| value of 3.3 (CL1279.Contig108_All) and 4.9 (CL7899.Contig3_All) in ML-vs.-W (table 2). From P1 to P5, one gene (CL9103.Contig5_All) significantly increased (|log2Ratio| = 2), and eight genes decreased (the maximum |log2Ratio| = 3.8) in expression. From W to P1 and P5 to P10, three and five genes were up-regulated, respectively (table 2).
Five genes (CL12980.Contig3_All, CL1585.Contig16_All, CL2596.Contig51_All, CL14397.Contig3_All, and CL16856.Contig2_All) were least-expressed in ML (fig. 5a). The expression levels of seven genes (unigene12087_All, CL16856.Contig2_All, unigene8645_All, CL18975.Contig2_All, CL10379.Contig10_All, CL9076.Contig3_All, and CL1337.Contig49_All) in W and P1 were higher than in other stages, while three genes (CL1585.Contig16_All, CL2596.Contig51_All, and CL1279.Contig108_All) were highly expressed in W (fig. 5a). The expression levels of three genes (Unigene15668_All, CL10689.Contig4_All, and CL1337.Contig49_All) in P5 were higher than in other stages. Three genes, including CL12980.Contig3_All, CL4397.Contig3_All, and CL9739.Contig5_All, were highly expressed in P10 while the expression level of unigene 15668_All was lower in P10 than in other stages (fig. 5a).
Figure 5. The heatmaps of differentially expressed genes (DEGs) were identified at different developmental stages during the prepupal–pupal transition of M. separata using Fragment Per Kilobase of exon model per Million mapped reads (FPKM) values. The gene IDs of cuticular proteins (a), juvenile hormone (b), and 20-hydroxyecdysone (c) are shown on the left side of the heatmaps and the developmental stages are given at the top of the heatmaps. ML, mature larvae; W, wandering; P1, P5, and P10: 1, 5, and 10 days after pupation. ML_1, ML_2, and ML_3 represent three biological replicates of ML stage. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Expression levels of genes involved in JH biosynthesis
There were ten and six down-regulated genes from ML to W and P1 to P5, respectively, and both of them had two up-regulated genes. Compared to W, four genes were significantly increased (the maximum |log2Ratio| = 3.8) and five genes were decreased (the maximum |log2Ratio| = 4.3) in P1. There were six up-regulated and three down-regulated genes in P5-vs.-P10; of these genes, CL18477.Contig4_All and CL9391.Contig2_All showed the maximum up-regulated and down-regulated |log2Ratio| values of 3.8 and 2.3, respectively (table 2).
Six genes (CL17942.Contig2_All, CL14984.Contig2_All, CL15275.Contig2_All, CL817.Contig4_All, CL4540.Contig4_All, and CL9391.Contig2_All) in the JH pathway were mainly expressed in ML (fig. 5b). The expression levels of CL18477.Contig4_All and CL11190.Contig1_All in W was higher than in other stages. Among five stages, CL5146.Contig76_All was highest expressed in P1 (fig. 5b). While CL17942.Contig2_All and CL817.Contig4_All were least expressed in P5, the expression level of CL16665.Contig4_All in P5 was higher than in other stages (fig. 5b). The expression levels of CL4540.Contig4_All and CL9391.Contig2_All in P10 were lower than in other stages. Meanwhile, unigene19996_All was mainly expressed in P10 (fig. 5b).
Expression levels of genes involved in 20E biosynthesis
Except for P1-vs.-P5 which has one up-regulated gene (CL17221.Contig1_All, |log2Ratio| = 2.1), there are no DEGs in the other three groups. A total of two (Unigene4431_All and CL17221.Contig1_All), one (CL14010.Contig2_All), and three (CL14010.Contig2_All, CL17221.Contig1_All, and CL19241.Contig2_All) genes were down-regulated in W-vs.-P1, P1-vs.-P5, and P5-vs.-P10, respectively (table 2). The expression of DEGs in other groups is shown in table S4. The expression levels of CL17221.Contig1_All in ML, W, and P5 were higher than in P1 and P10 (fig. 5c). Besides, CL19241.Contig2_All was mainly expressed in the P5 stage, and the expression level of CL14010.Contig2_All was lowest in P10 than in other stages (fig. 5c).
Other DEGs related to metamorphosis
In addition to the DEGs mentioned above, we also identified other 265 DEGs related to metamorphosis with a log2ratio value (|log2Ratio| ≥ 10). Among them, 28 and 11 genes were up- and down-regulated from ML to W, respectively (table S6); 13 and four genes were up- and down-regulated from W to P1, respectively (table S7); 12 and 18 genes were up- and down-regulated from P1 to P5, respectively (table S8); 147 and 42 genes were up- and down-regulated from P5 to P10, respectively (table S9).
Validation of the expression of DEGs using qRT-PCR
Correlation analysis revealed that the relative expression levels of 10 selected CP, 20E, and JH-related genes obtained by RT-qPCR corresponded well with their FPKM values derived from RNA-seq (fig. 6).
Figure 6. RT-qPCR validation of differentially expressed genes identified by RNA-Seq. Relative expression levels of ten genes, including four cuticular proteins (Unigene12087_All, CL2596.Contig51_All, Unigene8645_All, and CL10379.Contig10_All), three 20E genes (CL19241.Contig2_All, CL17221.Contig_All, and CL14010.Contig2_All), and three JH biosynthesis genes (CL17942.Contig2_All, CL4540.Contig4_All, and CL817.Contig4_All) among different stages in M. separata were analyzed. Yellow columns and blue polylines indicate the relative expression levels of differentially expressed genes identified by RNA-Seq and real-time PCR quantitative (RT-qPCR), respectively. The expression was normalized to transcript levels of reference genes β-actin and GAPDH. Letters above columns and polylines represent statistically significant differences by one-way analysis of variance (ANOVA) with Duncan's test (P < 0.05). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Discussion
Insect metamorphism is a series of changes in the external form, internal structure, physiological functions, living habits, and behavioral instincts of an insect during its development from larva to adult (Truman, Reference Truman2019). The pupal stage, characteristic of holometabolous insects, is a transitional morphology between larva and adult, represents the decomposition and reconstruction of tissue and organs (Jindra, Reference Jindra2019). Research on these metamorphic processes might provide molecular targets that can be utilized in pest control and to confer advantages to beneficial insects.
In this study, RNA-seq was used to analyze changes in expression levels of genes during the prepupal–pupal transition in the oriental armyworm M. separata. A total of 74,059 unigenes were identified and annotated in the transcriptome. Using equivalent insect developmental stages such as mature larvae, wandering, and pupation stages, a total of 11,301 unigenes were annotated in B. dorsalis (Chen et al., Reference Chen, Hou, Dou, Wei, Yue, Yang, Yu, De Schutter, Smagghe and Wang2018), and 22,197 in C. suppressalis (Sun, Reference Sun2015). These results highlight non-trivial differences in the number of genes expressed during prepupal–pupal transition across insect species. There are varied potential causes of this expression variance, such as functional differentiation of insect genes during evolution and environmental adaptation.
The epidermis covers the entire insect body surface, and not only defends against pathogen attacks and adverse environmental damage but also maintains the shape of the body and enables any required activities during development (Moussian et al., Reference Moussian, Schwarz, Bartoszewski and Nüsslein-Volhard2005). A previous study has shown that 50 CP genes were significantly expressed in the prepupal–pupal transition stage of B. mori (Shahin et al., Reference Shahin, Iwanaga and Kawasaki2016). In this study, 30 DEGs of CPs were identified, which is somewhat less than reported by Shahin et al. (Reference Shahin, Iwanaga and Kawasaki2016). The expression patterns of CP genes in different developmental stages or tissues of insects can give clues on functionality (Andersen et al., Reference Andersen, Hojrup and Roepstorff1995). For example, the expression of CP genes from larvae to pupae in B. mori was significantly different (Liang et al., Reference Liang, Zhang, Xiang and He2010). In this study, we found that eight CP genes in W and P1 were highly expressed, while two were highly expressed in W. The expression levels of three genes in P5 were higher than in other stages, and three genes were highly expressed in P10 (fig. 5a). Togawa et al. (Reference Togawa, Dunn, Emmons and Willis2007) found that for A. gambiae, five CP genes were only expressed before the pupa stage or adult molting (Togawa et al., Reference Togawa, Dunn, Emmons and Willis2007). These results suggest that CP genes are involved in the formation of the epidermis during the prepupal–pupal transition in insects.
JH is produced by the corpus allatum, and its main function is regulating the growth of insects by inhibiting the action of ecdysone (Zhang et al., Reference Zhang, Song, Li, Qian, Wei, Yang, Wang, Zhou, Meng, Peng, Xia, Perrimon and Cheng2018). JH plays an essential role in the growth of insect larvae, pupal metamorphosis, and adult reproduction (Raikhel et al., Reference Raikhel, Brown and Belles2005). In this study, 17 DEGs in the JH pathway were detected during the prepupal–pupal transition in M. separata. We found that the expression levels of six genes in the JH pathway decreased during the metamorphosis of mature larvae to pupae (fig. 5b). Our results are similar to those of Chen et al. (Reference Chen, Hou, Dou, Wei, Yue, Yang, Yu, De Schutter, Smagghe and Wang2018), who found that ten genes in the JH pathway were down-regulated during the process of larval metamorphosis in B. dorsalis. These results suggest that the expression of these genes in the JH pathway is to regulate deformation processes during these stages. The regulation of metamorphosis is partly the effect of JH titers (Niimi and Sakurai, Reference Niimi and Sakurai1997). Research on insect molting and metamorphic regulation has suggested that the metamorphosis of insects occurs when the JH titer decreases continually: the concentration of JH is higher in the larval stage, and drops to a mid- or low-level during pupation, and elevates subsequently in the adult stage (Bownes and Rembold, Reference Bownes and Rembold1987). Therefore, it needs to be determined whether the JH titer changes during the prepupal–pupal transition of M. separata.
Ecdysone plays a regulatory role in the growth, development, and reproduction of insects (Truman and Riddiford, Reference Truman and Riddiford2007). Ecdysone is a steroid compound that is secreted into the hemolymph and oxidized to 20E, and which can induce insect molting (Warren et al., Reference Warren, Petryk, Marques, Jarcho, Parvy, Dauphin-Villemant, O'Connor and Gilbert2002). In our study, a total of seven DEGs were highly expressed in the 20E pathway. This is similar to the result of Chen (Reference Chen2019) for B. dorsalis, who found that seven key genes of the 20E pathway were significantly overexpressed during the wandering stages. Zhao (Reference Zhao2019) showed that the expression levels of two ecdysone receptor genes (EcR and USP) in the 6th instar larvae of M. loreyi was higher than that in the 5th instar. Similar results were found in our study for CL17221.Contig1_All and CL19241.Contig2_All, which is higher in W than ML (fig. 5c). We found that the expression levels of CL17221.Contig1_All and CL19241.Contig2_All was higher in P5 than in P1 or P10 (fig. 5c). Previous research showed that the expression levels of four Halloween genes and three ecdysone receptor genes increased first and then decreased, in male pupae of M. loreyi (Zhao, Reference Zhao2019). These studies suggest that the expression of genes encoding ecdysone may co-regulate deformation processes during the insect larval–pupal transition. The titer of 20E in insects also changes along with the developmental age, and the change of insect 20E titer is key to the regulation of larval molting, pupal metamorphosis, reproduction, and embryo development of adults (Warren et al., Reference Warren, Petryk, Marques, Jarcho, Parvy, Dauphin-Villemant, O'Connor and Gilbert2002). Previous research has documented that there is a peak of 20E titer in the mid-stage of the embryo and another peak in the early stage of larval ecdysis; the higher concentration of 20E at the end of the larvae stage may be related to the initiation of transition to the pupal stage, while the higher levels of 20E in the early pupae may promote the transformation of larvae to pupae during the pre-pupal stage in Drosophilid (Thummel, Reference Thummel2001). How 20E titer changes during larval to pupal stages needs to be examined to further reveal the deformation processes of larvae–pupa in M. separata.
In this study, we used RNA-seq analysis in conjunction with RT-qPCR to provide a global analysis of the genes expressed from mature larvae to pupae of M. separata. Most genes in the prepupal–pupal transition transcriptome have been annotated according to reference information, and a large number of DEGs and pathways, especially the CPs, 20E, and JH pathways, were identified and their expression profiles characterized. These results suggest that M. separata prepupal–pupal transition is coordinated by hormonal regulation and is a result of a transcriptional shift of CP, 20E, and JH-associated genes (fig. 7). Collectively, the obtained transcriptome data herein not only provide a foundation for further understanding of the molecular mechanism of metamorphosis in M. separata, but also potentially help us to establish RNAi or CRISPR-based pest control strategies, which could operate through regulating the developmental processes of this pest.
Figure 7. A proposed model of 20E and JH biosynthesis in M. separata. ML, mature larvae; W, wandering; P1, P5, and P10: 1, 5, and 10 days after pupation. ML_1, ML_2, and ML_3 represent three biological replicates of ML stage. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Differentially expressed genes (DEGs), which are up- and down-regulated during the prepupal–pupal transition in M. separata in the pathway, are indicated with red and blue font followed by red and blue arrows, respectively.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485321000171
Acknowledgements
We thank Douglas Chesters (Institute of Zoology, Chinese Academy of Sciences, Beijing, China) for the language editing of this manuscript. This research was supported by the China Postdoctoral Science Foundation (2017M621160), National Major Science and Technology Projects of China (2019ZX08012004-011), the Innovative Talents Support Program of Liaoning Province (LCR2018024), and the Liaoning Revitalization Talents Program (XLYC1907097).
