Preliminary analysis on the developmental transcriptomes of swallowtail butterfly Papilio polytes (Lepidoptera: Papilioidae)

Abstract

The butterflies are one of the insect groups which undergo typical complete metamorphosis in their life history. However, up to now, little is known about the genomic mechanism that regulates the butterfly metamorphosis. In this study, using a swallowtail butterfly P. polytes as a model organism, a high-throughput sequencing platform was employed to perform the transcriptome and gene expression analyses in order to explore the P. polytes transcriptome features during different developmental stages. The results showed that approximately 398 million useful (Q20) reads were assembled into groups of 14698 (L1), 14264 (L2), 15084 (L3), 15520 (P1), 15052 (P2), 15720 (P3) and 15709 (A1), 14668 (A2), 16152 (A3) genes, respectively, with 58.19% to 67.11% of the data successfully mapped to the reference genome; the transcriptome change analysis via the DEGs Package revealed that dramatic gene expression differences were presented among the different developmental stages, that is, totally, 1162, 891 and 1723 genes were differentially expressed between adult and pupal, adult and larval, larval and pupal stages, respectively, with a number of these differentially expressed genes associated with the functions of digestion, cuticularization, chemoreception, wing formation, and so on. These differentially expressed genes and potential candidate genes required for butterfly metamorphosis by comparative transcriptomics may shed some new insights on molecular mechanisms underlying complete metamorphosis.

Introduction

Metamorphosis as a life-history strategies of insects, is an attractive and highly successful biological adaptation (Truman and Riddiford, 1999). The insects of complete metamorphosis such as butterflies, moths, flies and bees, undergo a series of continuous changes through their lifetimes, that is, in descending order they go through egg to larval, pupal, and adult stages in their different developmental periods (Resh and Cardé, 2009). The typical pupa which spans the larval and adult periods, is undergoing a series of intensely active changes to the last formation of various adult tissues and organs, including two pairs of wings (Xu, 2009), although the surface seems statically unchanged.

Insects are the most specious and evolutionarily successful groups in animals on earth, this is caused by the interaction of various factors, and among which the complex changeable metamorphosis is undoubtedly an important one (Malmstrom, 2012). In paleontological view, the occurrence and development of complete metamorphosis on the earth is probably related to the radically climatic transition between the Carboniferous and Permian. Changes in environmental factors can alter the activity of the endocrine system of organisms which can influence a series of physiological and biochemical changes, as well as morphological and behavioral alterations of animals. Thus the holometabolism (pupae formation from larvae to adults) may have developed rapidly due to the strong resistance of pupae to various environmental changes, especially, the rapid and violent climate change (Jia and Zhang, 1999).

In recent years, next-generation high-throughput DNA sequencing techniques have provided fascinating opportunities in the life sciences and dramatically improved the efficiency and cost-effectiveness of gene discovery (Ansorge, 2009), and for that reason these techniques has been widely used in the study of diverse organisms such as human beings (Lappalainen et al., 2013), plants (Garg and Jain, 2013; Ichihashi et al., 2015), agricultural pests ( Zhang et al., 2014; Zhang et al., 2016) to clarify the dynamic changes of gene expression and regulation during different developmental staged of these model and non-model organisms (Filichkin et al., 2010; Oppenheim et al., 2015). Up to the present, dozens of lepidopteran transcriptomes have been sequenced, however, only a few of developmental transcriptomes are available (Chen et al., 2010; Chen et al., 2012; Li et al., 2013; Qi et al., 2016; Zheng et al., 2012). Li et al. (2013) compared the transcriptomes of four developmental stages of the night moth Athetis lepigone and found that the most dramatic differences in gene expression were detected in the transitions from one stage to another stage. Some of these differentially expressed genes are related to cuticle and wing formation as well as the growth and development, and these expression differences were consistent with characteristic features of each life stage, for example, cuticular protein genes were up-regulated in larvae and down-regulated in pupae. Qi et al. (2016) sequenced the transcriptomes of different development stages of the small white butterfly Pieris rapae, and the results showed that a total of 21,883 differentially expressed unigenes were detected across the developmental stages, most of which were found between the egg and first larval stages, and total 32 heat shock protein (Hsp) genes with similar expression pattern were identified. In Hymenoptera, Chen et al. (2012) contrasted the transcriptomes of full-sister queen (QL) and worker-destined (WL) larvae of the honeybee using high-throughput RNA-Seq, and found that more than 4,500 genes with different expressions were existed between the two types larvae, of which more than 70% were upregulated in tha QL larvae.

The Common Mormon, Papilio polytes (Lepidoptera: Papilioidae) is a member of swallowtail butterfly species widely distributed in Asia. They undergo a series of striking changes through its lifetime as they mature from egg to adult (Chou, 1994). P. polytes as a kind of typical complete metamorphosis insects, the life of P. polytes experienced a very complex biological process, but little is known about the mechanism that regulates the development of P. polytes. Additionally, functional transcriptome information may lead to profound insights of the biology of P. polytes. In this study, taken P. polytes as a model organism, we generated a substantial dataset of transcript reads at their different developmental stages from larvae to adults, by using Illumina NextSeq500 RNA-Seq technique, in order to reveal the gene expression characteristics and clarify the functions of related genes or gene families in their different developmental stages.

Materials and Methods

Materials

The larval, pupal and adult individuals of P. polytes were collected at the campus of Anhui Normal University (Wuhu, China) in May and June of 2017. Test samples used in this study were: (1) a tissue mixture of third instar larvae; (2) a tissue mixture of pupae and (3) a tissue mixture of adults. In this study, adult male P. polytes were used for this study. Each sample after collection was immediately frozen in liquid nitrogen and stored at -80oC until RNA extraction. Three biological repeats for each stage: larvae (L1, L2, L3), pupae (P1, P2, P3) and adults (A1, A2, A3) were set in the experiments to exclude the artificial errors.

RNA isolation, cDNA library preparation and illumina sequencing

Total RNA from tissues of larvae, pupae and adults was extracted separately according to the manufacturer’s protocol. RNA integrity was confirmed by the Agilent 2100 Bioanalyzer (Agilent Technologies), and cDNA library construction and Illumina sequencing were subsequently performed at Personal Biotechnology Co., Ltd., Shanghai China. Briefly, Magnetic beads with oligo-dT were used to combine the poly-A of the mRNA for purifying the mRNA from the total RNA. The mRNA was then mixed with fragmentation buffer to obtain 200-300 bp short fragments. The fragments were used to synthesize first-strand cDNA with random primers, and first-strand cDNA was transformed into double-strand cDNA by using RNase H and DNA polymerase I. When the second chain cDNA began to syntheze, the bases of these T should be replaced with U, so as to achieve the aim of chain specific library. Library fragments were purified with a QIAquick PCR Extraction Kit (Qiagen, Valencia, CA, USA), and the sequencing library was constructed using polymerase chain reaction (PCR). The multiplex cDNA libraries were checked using PicoGreen (Quantifluor™-ST fluorometer E6090, Promega, CA, USA) and fluoro spectrophotometry (Quant-iT PicoGreen dsDNA Assay Kit; Invitrogen, P7589) and quantified with Agilent 2100 (Agilent 2100 Bioanalyzer, Agilent, 2100; Agilent High Sensitivity DNA Kit, Agilent, 5067–4626). The synthesized cDNA libraries were then normalized to a 10 nM, and the sequencing library was gradually diluted and quantified to 4–5 pM and sequenced using the Illumina NextSeq™ 500 platform instrument with paired-end (PE) sequencing method.

Gene expression analysis

In order to minimize sequencing mistakes, adaptor sequences and low reads were filtered out prior to mapping, the reference genome index was established by Bowtie2 (Langmead et al., 2009), and then the filtered reads were compared to the reference genome using Tophat2 (Kim et al., 2013) (http://tophat.cbcb.umd.edu/). The “RPKM” method (reads per kilo bases per million reads) was used for normalization with the total mapped reads count and the gene length to quantitate the level of gene expression, and results less than 0.005 were excluded (Mortazavi et al., 2008). The significance was determined by normalizing the raw reads and calculating the P-value using DESeq (Anders, 2012). Genes with fold change>2 or <0.5 and P-value<0.05 were identified as differentially expressed gene (DEG) (Anders and Huber, 2010).

The Blast2GO (Conesa et al., 2005) program were used to compare genes with GO (Gene Ontology) database (http://www.geneontology.org/), and the GO terms for each gene were submitted to Web Gene Ontology Annotation Plot (WEGO) for further classification (Ashburner et al., 2000). The genes significantly different from the corresponding library were searched against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to determine the pathways (Kanehisa et al., 2004).

Results and Discussion

Raw data processing

In total, 404 million reads were generated from the nine libraries, and 398 million useful (Q20) reads were selected for further analysis. The useful reads were mappted to the P. polytes genome (Nishikawa et al., 2015) using tophat software and the results showed that 58.19% to 67.11% of reads were successfully mapped. Among them, 69.91% to 87.67% could be mapped to unigenes, 95.13% to 96.59% were successfully mapped to exonic regions (Table 1). Totally, 16360 genes were identified, and among these genes, 14698, 14264, 15084, 15520, 15052, 15720, 15709, 14668 and 16152 genes were identified in the L1, L2, L3, P1, P2, P3, A1, A2, A3 libraries, respectively (Table 1).

Table 1. Summary statistics of clean reads in the P. polytes transcriptomes.

The DEGs of three stages were analyzed through heat maps, and hierarchical clusters were generated to obtain a global view of the DEGs which were illustrated using a grid of colored tiles. A crosswise comparison of 9 samples were performed through cluster analysis, and the results showed that the samples from the same stage were clustered into the corresponding group; the genes detected in the different groups were clearly separated; the increase or decrease in transcript abundance in the pupae and adults was significantly different from those in the larvae; in addition, P1, P2 and P3 were basically the same, L1, L2 and L3 were similar in general, while A1, A2 and A3 were somewhat different with each other. These results indicated that the gene expressions changed dramatically during different developmental stages from larvae to adults, and overally larvae and aldults harbored higher gene expression abundances than those in the pupae (Figure 1).

Differentially expressed genes (DEGs) during P. polytes development

Differences in gene expression and DEGs were examined and identified at three stages during P. polytes development by pairwise comparisons of the three libraries. The numbers of up- and down-regulated genes were also calculated between each pair of life-stages. The statistical results were shown in Figure 2 and Figure 3, with the top 10 up and down regulated genes listed in Table 2. Apparently, over 70% of the DEGs were more highly expressed in the larvae (L) and adults (A) than in the pupae (P).

By comparison between the larval and pupal transcriptomes, 1,723 significantly differentially expressed genes were detected, including 1280 (74.3%) up-regulated and 443 (25.7%) down-regulated genes in the larval library. Seven of the top ten up-regulated genes in larval transcriptomes were matched with the orthologs of P. polytes genome: one uncharacterized protein LOC106104536, six with Predicted functions (bile salt-activated lipase-like, trypsin, alkaline C-like, membrane alanyl aminopeptidase-like, fibroin light chain-like, aminopeptidase N-like and ejaculatory bulb-specific protein 3-like precursor); the other two were matched with the sodium-dependent nutrient amino acid transporter 1 and uncharacterized protein LOC106119902 in Papilio xuthus respectively; the last one was matched with of the membrane alanyl aminopeptidase in Papilio machaon. Besides, nine of the top ten down-regulated genes in larval transcriptomes (mushroom body large-type Kenyon cell-specific protein 1 isoform X2, delta-like protein 1 isoform X1, putative cyclin-dependent serine/threonine-protein kinase DDB_G0272797/DDB_G0274007 isoform X1, uncharacterized protein LOC106100118, LOC106100129, LOC106099981, LOC106108048, LOC106109383 and LOC106108920) were matched with P. polytes, one was matched with the hypothetical protein RR46_06374 of P. xuthus (Table 2).

Table 1. The top 10 up or down-regulated genes identified among Pupae (P) vs. larvae (L), larvae (L) vs. adults (A), and pupae (P) vs. adults (A).

By comparison between the adult and pupal transcriptomes, 1,162 genes with significantly differential expression levels were detected in adults, and among them, 873 (75.1%) up-regulated genes and 289 (24.9%) down-regulated genes were found. Seven of the top ten up-regulated genes were matched with the synaptic vesicle glycoprotein 2A-like, retinoid-inducible serine carboxypeptidase-like, carboxypeptidase A1-like, pancreatic triacylglycerol lipase-like, serine proteases 1/2-like, uncharacterized protein LOC106105328 and ammonium transporter Rh type B-B in P. polytes; three genes were matched with the hypothetical protein RR46_02520 in P. xuthus, hypothetical protein RR48_11198 and uncharacterized protein LOC106709423 in P. machaon respectively. The top ten down-regulated genes in the adult library included the following Predicted functional genes: protein doublesex isoform X2 in P. xuthus; D-beta-hydroxybutyrate dehydrogenase, mitochondrial; basic juvenile hormone-suppressible protein 2-like, basic juvenile hormone-suppressible protein 1-like and actin cytoskeleton-regulatory complex protein PAN1-like, and four uncharacterized protein LOC106100118, LOC106100129, LOC106107989, LOC106111008, LOC106108048 in P. polytes (Table 2).

A total of 891 genes (including 359 (40.3%) up-regulated genes and 532 (59.7%) down-regulated genes) were detected between larval and adult transcriptomes. Among the top ten up-regulated gene in adults, eight genes were matched with P. polytes, UDP-glucuronosyltransferase 2C1-like, uncharacterized protein LOC106101520, LOC106101491, nose resistant to fluoxetine protein 6-like, mushroom body large-type Kenyon cell-specific protein 1 isoform X2, inducible metalloproteinase inhibitor protein-like, protein msta, isoform A and troponin C isoform 1 respectively; the other two up-regulated genes were matched with the female-specific doublesex isoform F1 in silkworm Antheraea mylitta and hypothetical protein RR48_14296 in P. machaon. Among the top ten down-regulated genes in the adult transcriptome, five genes were the 4-coumarate-CoA ligase 1-like, trypsin, alkaline C-like, gastrolith matrix protein-like, carboxypeptidase B-like, flexible cuticle protein 12-like, glycine-rich protein 1-like, bile salt-activated lipase-like in P. polytes; the other two genes with unknown functions were the uncharacterized proteins LOC106111207 and LOC106105319 in P. polytes; the last one was matched with the bilin-binding protein-like precursor/bilin binding protein 1 in P. polytes (Table 2).

Functional classification of DEGs during P. polytes development

WEGO annotation of DEGs: We used GO assignments to classify the functions of DEGs in pairwise comparisons of cDNA libraries between different P. polytes developmental stages. These genes were divided into three categories: molecular function (MF), cellular component (CC) and biological process (BP). In total, 67 categories were subdivided from the primary categories: 24 categories for biological process, 24 categories for cellular component and 19 categories for molecular function. In the three GO categories, the highest proportion of annotations are cellular process and localization in the category of biological process, cell and macromolecular complex in the cellular component, catalytic activity and binding in the molecular function (Figure 4). The results showed that these significantly enriched biological functions played a key role during the growth and development of the butterflies.

Pathway annotation of DEGs between developments: To understand the functions of the DEGs, we mapped them using the KEGG database for signaling pathways analysis. The top 20 most abundant differentially expressed signaling pathways from each compared group were listed in Figure 5. The results showed that in all the three groups, the cluster for metabolic pathways was the largest group, within which both of the fatty acids and cytochrome P450 (CYPs) were significantly enriched. Fatty acids (an important source of fuel for muscular contraction and general metabolism) and CYPs (the major enzymes involved in drug metabolism and participate in the synthesis of steroid hormones) are all closely related to the protein translation, processing modification, hormone synthesis and other aspects of metabolism in the development of P. polytes.

DEGs in juvenile hormone pathway

Juvenile hormone (JH) and ecdysterone are the decisive endogenous factors that induce the development of insects (Dubrovsky, 2005; Jindra & Riddiford, 2013; Truman & Riddiford, 1999, 2002). JH prevents the metamorphosis induced by the ecdysone, keeping the larvae in the larval condition (Truman and Riddiford, 2007). In the P. polytes JH pathway of this study, two enzymes were detected to have a strong effect on the JH titre, one is JHAMT (juvenile hormone acid methyltransferase) that plays a key role in regulates the last step of JH biosynthesis (Wijesekera & Dauwalder, 2016), the another is JHEH (juvenile hormone epoxide hydrolase) that is involved in primary degradation of JH. Out study showed that the JHEH gene was obviously up-regulated in larvae (L) compared with papue (P), resulting in a relatively lower JH titre in larvae (plot A, Figure 6); the JHAMT gene was not shown to be up-regulated, whereas the JHEH shown to be down-regulated in adults (A) compared with larvae (L), resulting in a higher JH titre in adults (plot B, Figure 6). These results were congruent with those of previous studies in other insect groups (Yang et al., 2011; Wijesekera and Dauwalder, 2016).

Conclusion

The P. polytes exhibits the typical developmental characteristics of holometabolous metamorphosis of insects. Our study showed that during the transitions from one to another developmental stages, the gene expressions were significantly correlated with the developmental characteristics of P. polytes. During the transition from larvae to puape, the most up-regulated genes such as membrane alanyl aminopeptidase (aminopeptidase N, APN) and ejaculatory bulb-specific protein 3-like precursor, are closely related to immune, digestive, and chemoreceptive functions respectively. For example, APNs expressed in the midgut of the insect larvae and primarily involved in dietary protein digestion (Wang and Zhang, 2005) and as receptors in Cry toxin-induced pathogenesis in insects (Bravo and Soberón, 2007); the putative cyclin-dependent serine/threonine-protein kinase was highly expressed in the pupae than in the larvae, it is a key factor in the control of cell division by modulating transcription in response to several extra- and intracellular activities in this radical transitional process (Malumbres, 2014). During the transition from puape to adults, the genes responsible for digestion, defense response, development and reproduction, such as the retinol-binding protein gene and serine proteases 1/2-like were remarkably up-regulated in the adults. For example, serine proteases are the principal digestive enzymes which provide nutrients needed for survival and fecundity (Bao et al., 2014; Zou et al., 2006), however, some genes responsible for juvenile-state keeping and metamorphosis preventing, such as the juvenile hormone-suppressible protein 2-like/1-like were highly expressed in the pupal stage, which was consistent with the results of previous studies in hemimetabolous locust Schistocerca gregaria by juvenile hormone treatment (Novak, 1969; Micciarelli, 1977; Riddiford, 2012). When larvae and adults were compared, the nose resistant to fluoxetine and the troponin C (TpnC) genes which regulates the muscle movements and carries out the flight functions were highly expressed at the adult stage (Bullard and Pastore, 2011; Herranz and Marco, 2005; Kržič et al., 2010).

In general, our study constructed nine transcriptome libraries and generated new genetic resources for the investigation of P. polytes development. This latest molecular research as well as many other in-depth studies about insect metamorphosis will shed some new lights on the underlying genomic mechanisms (Liu et al., 2017) and related evironments (Helm et al., 2017) about these complicated biological processes.

Acknowledgments

This work was supported by the Special Funds for Cultivation of Preponderant Disciplines of Anhui Normal University and funds from the State Key Laboratory of Palaeobiology and Stratigraphy (Nanjing Institute of Geology and Palaeontology, CAS) (Grant No. Y626040108).