Reconstruction of insect hormone pathways in an aquatic firefly, Sclerotia aquatilis (Coleoptera: Lampyridae), using RNA-seq

RNA extraction and sequencing

All Sclerotia aquatilis specimens used in this study were provided by Anchana Thancharoen (firefly rearing laboratory, Department of Entomology, Faculty of Agriculture, Kasetsart University). Total RNA extractions were performed according to Vongsangnak, Chumnanpuen & Sriboonlert (2016). The RNA extractions were prepared from six specimens for each developmental stage: larva (six specimens of 5th instar), pupa (six specimens), and adult (four female specimens and two male specimens). Two samples were pooled to make three replicates for cDNA library constructions. Quality and quantity of RNA samples were measured using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Preparation of cDNA library and sequencing were performed according to Vongsangnak, Chumnanpuen & Sriboonlert (2016). RNA sequencing was performed using an Illumina Hiseq2000 Sequencing System (Illumina, San Diego, CA, USA) to generate paired-end reads with 100 base pairs (bp) read length.

De novo transcriptome assembly

Low quality reads (i.e., reads with more than 20% of the base qualities being lower than 10, reads with adaptors and reads with unknown bases [N bases greater than than 5%]) were removed with FASTQC and Trimmomatic. Clean reads were assembled into Unigenes using Trinity v.2.0.6. (https://github.com/trinityrnaseq/trinityrnaseq), with parameter setting: –min_contig_length 150 (minimum 150 base of assembled contig length to report) –CPU 8 (8 CPUs to use) –min_kmer_cov 3 (minimum 3 count for K-mers to be assembled by Inchworm) –min_glue 3 (minimum 3 reads needed to glue two inchworm contigs together) –bfly_opts ’-V 5 (additional parameters to pass through to butterfly) –edget-hr=0.1 (set 0.1 as value to butterfly command) –stderr’ (standard error command to butterfly). Gene family clustering was subsequently performed using Tgicl clustering software v2.0.6.

Functional annotation of transcripts

Functional annotation was performed against NCBI non-redundant protein sequences (NR; ftp://ftp.ncbi.nlm.nih.gov/blast/db), NCBI non-redundant nucleotide sequences (NT; ftp://ftp.ncbi.nlm.nih.gov/blast/db), clusters of orthologous groups of proteins (KOG/COG; ftp://ftp.ncbi.nih.gov/pub/COG/KOG), Swiss-Prot (http://ftp.ebi.ac.uk/pub/databases/swissprot) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg), Gene Ontology (GO; http://geneontology.org) and IntroPro (v5.11–51.0; http://www.ebi.ac.uk/interpro/interproscan.html) databases using Blastn, Blastx (v.2.2.23; http://blast.ncbi.nlm.nih.gov/Blast.cgi), Diamond, Blast2GO (v2.5.0; https://www.blast2go.com) and InterPro Scan5 (v5.11–51.0; https://code.google.com/p/interproscan/wiki/Introduction). The results were filtered using an E-value of 1E-05 and a sequence identity of greater than 10%, and the best hit for each transcript was retained. Transdecoder (v.3.0.1; https://transdecoder.github.io) was used to identify the candidate coding area. The longest open reading frame (ORF) was extracted, and then the Pfam protein homologous sequences were searched by using the BLAST tool (v.2.2.23; http://blast.ncbi.nlm.nih.gov/Blast.cgi) in the SwissProt database (http://ftp.ebi.ac.uk/pub/databases/swissprot) and Hmmscan to predict the coding region. These programs were run using default parameter settings. All annotated transcripts were used to develop a Venn diagram, using jvenn (http://jvenn.toulouse.inra.fr/app/example.html) to represent the number of annotated transcripts expressed during the three stages of S. aquatilis development and other analyses.

Comparative analysis of Unigenes among firefly species

The annotated transcripts which were each assigned a SwissProt ID based on the Swiss-Prot database and were used to compare the transcripts against the SwissProt IDs from three fireflies species, Asymmetricata circumdata, Aquatica ficta, and Pyrocoelia pectoralis (Wang et al., 2017), to identify orthologous genes among fireflies species.

Identification of genes involved in insect hormone pathways

The amino acid sequences translated by Transdecoder were mapped to the KEGG database using Ghost KEGG Orthology and Links Annotation (KOALA) (https://www.kegg.jp/ghostkoala/). Those sequences that were mapped to the biosynthetic pathways for insect hormones were collected for further analysis. The expression profiles for these genes were visualised as a Heatmap, generated based on fragments per kilobase of transcript per million mapped reads (FPKM) values, using Morpheus, (https://software.broadinstitute.org/morpheus/) with a Pearson correlation hierarchical clustering, average linkage method. The significant differences in the FPKM values among the three stages were examined using a one-way analysis of variance (ANOVA) statistics and a post hoc test (Scheffe method) under the 95% confidence level, using SPSS software (IBM SPSS Statistics for Windows, Version 22.0; IBM Corp., Armonk, NY, USA).

These candidate genes were further analyzed to identify the most probable transcripts for each enzyme in the pathway. The candidate genes were BLASTed against the amino acid sequences of other insect species, i.e., T. castaneum (KEGG genes ID: tca:661876, tca:661964, tca:662004, tca:100142466, tca:656859, tca:657029, tca:657192, tca:657270, tca:657342, tca:657419, tca:657969, tca:658287, tca:658434, tca:658587, tca:658671, tca:658696, tca:659209, tca:659494, tca:103312192, tca:103314817, tca:656328, tca:659438, tca:661224, tca:661279, tca:100141982, tca:100142187, tca:662961, tca:658858, tca:658208, tca:100313955, tca:659305, tca:659375, tca:659506, tca:659572, tca:658081, tca:663127, tca:656884, tca:663098, tca:658665, tca:661451, tca:656794), D. melanogaster (sequence ID: CG15739, AF484413, AY079170, AF484414, U44753, NP_651725; Nyati et al., 2013; Warren et al., 2004; Warren et al., 2002; Petryk et al., 2003; Bassett et al., 1997), A. aegypti (sequence ID: GQ344797.1, GQ344798.1, GQ344799.1, KC243497, KC243498, KC243495, KC243496, KC243499; Mayoral et al., 2009; Rivera-Perez et al., 2013), B. mori (sequence ID: AB124839, AB113578, AY377854, AY363308, AF287267, ABB58822, AB198340; Daimon et al., 2012; Shinoda & Itoyama, 2003; Zhang et al., 2005; Li et al., 2005; Hirai et al., 2002; Rewitz, O’Connor & Gilbert, 2007; Niwa et al., 2005), D. punctata (sequence ID: AY509244 and AY509245; Helvig et al., 2004), C. fumiferana (sequence ID: AF153367; Hirai et al., 2002) and H. virescens (sequence ID: AF037196; Hirai et al., 2002). The best possible candidates for each enzyme-coding gene were selected based on the following criteria: BLASTP identity value greater than 20%; 1E-50 E-value cut off; FPKM values that were significantly different among the three developmental stages (tested by one-way ANOVA with p < 0.05); and positive detection in three other firefly species A. circumdata, A. ficta, and P. pectoralis (Wang et al., 2017). The identified candidates were then subjected to real-time PCR, RT-PCR, and sequencing.

Validation of candidate genes expression using real-time PCR

The expression profiles of the selected candidate genes were examined using real-time PCR. Three individual specimens were used for each developmental stage (three fifth instar larvae of unknown sex; three pupae of unknown sex; three adults (two females and one male). Each specimen was pulverized into a powder with liquid nitrogen. Total RNA was extracted using TRIzol reagent (Invitrogen, Waltham, MA, USA), according to manufacturer’s protocol. The total RNA was subsequently treated with DNase I (Qiagen, Hilden, Germany) as per manufacturer’s protocol. Total RNA was converted into cDNA using RevertAid first strand cDNA kit (Thermo Scientific, Waltham, MA, USA). The reverse transcription reaction included 2 µg of total RNA, 10mM dNTP, 5x reaction buffer, RiboLock RNase Inhibitor, and RevertAid M-MuLV RT as described in manufacturer’s protocol and the reaction was incubated for 60 min at 42 °C and terminated at 70 °C for 5 min. Quantitative real-time PCR was performed using an Eppendorf Mastercycler ep realplex4 S real-time RT-PCR instrument (Eppendorf, Hauppauge, NY, USA) with Maxima SYBR green qPCR master mix (Thermo Scientific, Waltham, MA, USA). The quantitative real-time PCR reaction included cDNA template 0.5 µL, 0.3 µM of forward/reverse primers (File S1), Maxima SYBR green qPCR Master Mix (2x) no ROX 5 µL and water was added up to 10 µL. Actin (CL5997) was amplified as an internal control. The PCR reaction was initial denaturation for 10 min at 95 °C, followed by denaturation for 15 s at 95 °C, annealing for 30 s at 59 °C, and the last extension for 30 s at 72 °C. Three replicates were performed for each sample. The CT value from Real-time PCR was used to calculate a relative expression using the Δ ΔCT method. Significant differences of relative expression value among stage were tested with one-way ANOVA statistics and a post hoc test (Scheffe method) at p-value <0.05.

Validation of candidate genes expression using RT-PCR and sequencing

The full-length sequences of the candidate genes were examined using RT-PCR, followed by sequencing. Complementary DNA obtained from the previously described real-time PCR experiment was used as the template during RT-PCR analysis. RT-PCR primers (File S1) were designed to amplify the longest ORFs in each candidate gene. PCR was performed using KOD-PLUS-NEO DNA polymerase (TOYOBO, Japan), according to manufacturer’s protocol. The PCR was performed using a Bio-rad T100™ Thermal Cycler (Bio-Rad, Hercules, CA, USA), with the following conditions: initial denaturation for 2 min at 94 °C, followed by 30 cycles of denaturation for 10 s at 98 °C, annealing for 30 s at 50 °C, and extension for 1 min at 68 °C. PCR products were subsequently sent for sequencing at Macrogen (South Korea).

De novo transcriptome assembly

A total of 711.9 Mb of raw reads were obtained from the Illumina HiSeq 2000 sequencing of nine specimens across the three developmental stages of S. aquatilis (Table 1). After adapter trimming and the removal of contaminating sequences and low-quality sequences, a total of 592.32 Mb of clean reads were retained. These clean reads were assembled into 383,922 high-quality contigs. Unigenes were obtained by removing the transcript abundance. The GC contents for all samples ranged from 36–38%. A total of 82,022 unigenes were retained, with a total length, average length, N50, and GC content of 117,432,935 bps, 1,394 bps, 2,666 bps, and 35.65%, respectively. The clean reads from the nine specimens of S. aquatilis used in this study were deposited into the NCBI SRA database under accession number PRJNA525262.

Gene annotation

Of these 82,022 assembled unigenes, 46,230 (56.36%) were annotated in the following seven databases: Nr (44,959), Nt (9,908), SwissProt (34,626), KEGG (35,393), KOG (33,407), InterPro (35,813) and GO (9,278) databases (File S2). These annotated transcripts were categorized into three primary functional groups: biological process, cellular component, and molecular function, according to GO annotations (Fig. 1). The functional category that contained the highest number of unigenes was the biological process (16,752). Within this biological process category, most of the unigenes were subcategorized into the metabolic (3,784) or cellular (3,932) processes. The results also showed that the number of unigenes among the three stages were distributed similarly across the various subcategories; for example, the number of unigenes in the metabolic process categories larval, pupal, and adult stages were 3,556, 3,019, and 3,012, respectively.

All of the annotated transcripts obtained from all of the databases were then combined and compared among the three different developmental stages of S. aquatilis. A total of 34,899 unigenes were found to be commonly expressed across three different developmental stages. In contrast, 4,793, 1,039, and 794 unigenes were specifically expressed during the larval, pupal, and adult stages, respectively (Fig. 2). The larval and pupal stages shared 1,253 unigenes, the pupal and adult stages shared 1,946 unigenes and the larval and adult stage shared 1,506 unigenes. These specifically expressed unigenes were then catagorized into functional catagories, based on the KEGG database (Fig. 3). During the larval stage, the highest numbers of specifically expressed unigenes were identified in the genetic information processing category (840 unigenes) and the translation subcategory (471 unigenes). During the pupal stage, the highest number of specific unigenes was identified in the metabolism category (91 unigenes). During the adult stage, the highest numbers of specific unigenes were identified in the cellular processes (63 unigenes) and the cellular community subcategory (28 unigenes). The number of specifically expressed unigenes associated with lipid metabolism (seven unigenes) was relatively high during the pupal stage compared with the larval and adult stages.

Comparative analysis of Unigenes among firefly species

Orthologous genes among four firefly species, A. circumdata, A. ficta, P. pectoralis, (Wang et al., 2017) and S. aquatilis (this study) were identified. The Swissprot IDs of the three firefly species, A. circumdata, A. ficta, and P. pectoralis were obtained from (Wang et al., 2017) and compared with S. aquatilis unigenes in the present study. The analysis showed 1,464, 1,596, 1,589, and 10,246 unigenes found specifically in A. circumdata, A. ficta, P. pectoralis, and S. aquatilis, respectively. A total of 31,150 unigenes were found in all four species (Fig. 4; File S3).

Identification of genes involved in insect hormone pathways

To identify the protein coding genes involved in the pathways of insect hormones, all unigenes identified in S. aquatilis were mapped to the KEGG pathway database. A total of 221 transcripts were identified as being involved in the biosynthetic and degradation pathways of insect hormones (Table 2; File S2) with percent identities greater than 20 and E-value of 1E-05 set as the thresholds (Figs. 5 and 6). Of these 221 transcripts, 191 transcripts were identified in the JH biosynthetic and degradation pathways and 30 transcripts were identified in the ecdysteriod biosynthetic and degradation pathways. In the JH biosynthetic pathway, 130 transcripts were mapped to five enzymes, including farnesyl diphosphate phosphatase (FPPP; nine transcripts), NADP⁺-dependent farnesol dehydrogenase (FOHSDR; 56 transcripts), aldehyde dehydrogenase (ALDH; 23 transcripts), juvenile hormone-III synthase (JHAMT; 29 transcripts), and methyl farnesoate epoxidase/farnesoate epoxidase (CYP15A1; 13 transcripts). In JH degradation pathway, 61 transcripts were mapped to three enzymes comprising of juvenile-hormone esterase (JHE; 56 transcripts), juvenile hormone epoxide hydrolase (JHEH; three transcripts), and juvenile hormone diol kinase (JHDK; two transcripts). These enzymes were similar to the enzymes identified in the model coleopteran insect, Tribolium, except for the JHDK, which was not detected in the Tribolium but was found in Bombyx mori (Li et al., 2005) and Manduca sexta (Maxwell et al., 2002; Maxwell, Welch & Schooley, 2002). The differential expression of these 191 transcripts were also compared among the three developmental stages, larval, pupal and adult. FPKM values from the 191 transcripts were analyzed using one-way ANOVA statistics and a post hoc test (Scheffe method), at a confidence level of 95%. Only 34 transcripts were found to be significantly different.

In the ecdysteriod pathway, 30 transcripts were mapped to eight enzymes; cytochrome P450 family 307 subfamily A (Spook [Spo]; 1 transcript), cytochrome P450 family 307 subfamily B (Spookiest [Spot]; 1 transcript), short-chain dehydrogenase/reductase (Shroud [Sro]; 3 transcripts), ecdysteroid 25-hydroxylase (Phantom [Phm]; 1 transcript), ecdysteroid 22-hydroxylase (Disembodied [Dib]; 16 transcripts), ecdysteroid 2-hydroxylase (Shadow [Sad]; 1 transcript), ecdysone 20-monooxygenase (Shade [Shd]; 4 transcripts), and 26-hydroxylase (Cyp18a1; three transcripts) (Fig. 6). Almost all members of the insect ecdysterol Halloween gene family (Spo, Spot, Phm, Dib, Sad, and Shd) were identified in S. aquatilis, except for Neverland (Nvd). According to the KEGG pathway, this Nvd enzyme is also absent from the model coleopteran insect, T. castaneum. The expression profiles of these 30 transcripts were analyzed using their FPKM values. The significant differences among the differential expression profiles were analyzed using one-way ANOVA. However, the expression levels of all these transcripts were not significantly different among the three stages, except for Sro (Unigene21593).

Of these 191 transcripts identified in the JH pathway and the 30 transcripts in the ecdysone pathway, only twenty most probable candidates from the JH (eleven transcripts) and the ecdysone (nine transcripts) pathways were selected for further analyses. These twenty candidates were chosen based on at least one of the following criteria: (i) the protein coding gene transcript was detected in all four firefly species (S. aquatilis [this study], A. circumdata, A. ficta and P. pectoralis (Wang et al., 2017)); (ii) the gene was returned as a best hits when BLASTed against the Tribolium hormone enzymes; and (iii) the gene was differentially expressed among the three stages of development.

Validation of differential gene expression by quantitative real-time PCR, RT-PCR, and sequencing

The expression profiles of these twenty potential candidates were validated using real-time PCR. Actin (CL5997) was used as a reference gene, due to its relatively stable expression in comparison with the other house keeping genes identified in this study. Real-time PCR results for most of the selected genes coincided with the expression data obtained from the transcriptome analysis (Fig. 7; File S4). However, the real-time PCR expression patterns for FPPP (Unigene19505), JHEH (Unigene9157), Phm (Unigene12481), and Cyp18a1 (CL2776) among the three stages of S. aquatilis development slighlty deviated fom the expression patterns obtained from the transcriptome analysis. Moreover, RT-PCR and Sanger sequencing were also performed to analyze the full-length coding sequences of these transcripts. RT-PCR products of an expected sizes were successfully amplified for all of the selected candidates, except for Unigene12798 (Spot) (Figs. 8A and 8B). The nucleotide sequences for most of these candidates were nearly identical to the sequences obtained from the transcriptome analysis (File S5). However, the RT-PCR product for Unigene12798 (Spot) showed two bands of approximately 2,000 and 750 bp, instead of the expected 1,476 bp band. Sanger sequencing results for both Unigene12798 (Spot) and Unigene19886 (Spo) showed mixed peaks and were unable to be read. Thus, the primers for these two genes were redesigned and, subsequently, the partial coding sequences of spo (1,038bp) and spot (1,011bp) were able to be amplified and sequenced (Fig. 8C).