Long-Read Genome Sequencing of Abscondita cerata (Coleoptera: Lampyridae), the Endemic Firefly of Taiwan

© 2023 Academia Sinica, Taiwan Open Access Long-Read Genome Sequencing of Abscondita cerata (Coleoptera: Lampyridae), the Endemic Firefly of Taiwan Tzi-Yuan Wang1, King-Siang Goh2, Liang-Jong Wang3, Li-Ling Wu4, Feng-Yu Wang5,* , and Yu-Wei Wu6,7,8,* 1Biodiversity Research Center, Academia Sinica, Taipei 115, Taiwan. E-mail: [email protected] (TY Wang) 2Genomics Research Center, Academia Sinica, Taipei 115, Taiwan. E-mail: [email protected] (Goh) 3Forest Protection Division, Taiwan Forestry Research Institute, Taipei 100, Taiwan. E-mail: [email protected].tw (LJ Wang) 4Department and Institute of Physiology, College of Medicine, National Yang Ming Chiao Tung University, Taipei 112, Taiwan. E-mail: [email protected] (LL Wu) 5Taiwan Ocean Research Institute, National Applied Research Laboratories, Kaohsiung 852, Taiwan. *Correspondence: E-mail: [email protected] (FY Wang) 6Graduate Institute of Biomedical Informatics, College of Medical Science and Technology, Taipei Medical University, Taipei 106, Taiwan. *Correspondence: Email: [email protected] (YW Wu) 7Clinical Big Data Research Center, Taipei Medical University Hospital, Taipei 235, Taiwan 8TMU Research Center for Digestive Medicine, Taipei Medical University Taipei 110, Taiwan Received 27 August 2022 / Accepted 17 February 2023 / Published 26 May 2023 Communicated by Chih-Ming Hung Abscondita cerata is the most abundant and widely distributed endemic firefly species in Taiwan and is considered a key environmental and ecological indicator organism. In this study, we report the first longread genome sequencing of Abs. cerata sequenced by Nanopore technology. The draft genome size, 967 Mb, was measured through a hybrid approach that consisted of assembling using 11.25-Gb Nanopore long reads and polishing using 9.47-Gb BGI PE100 short reads. The drafted genome was assembled into 4,855 contigs, with the N50 reaching 325.269 kb length. The assembled genome was predicted to possess 55,206 protein-coding genes, of which 20,862 (37.78%) were functionally annotated with public databases. 47.11% of the genome sequences consisted of repeat elements; among them DNA transposons accounted for the largest proportion (26.79%). A BUSCO (Benchmarking Universal Single Copy Orthologs) evaluation demonstrated that the genome and gene completeness were 84.8% and 79%, respectively. The phylogeny constructed using 1,792 single copy genes was consistent with previous studies. The comparative transcriptome between adult male head and lantern tissues revealed (1) the vision of Abs. cerata is primarily UV-sensitive to environmental twilight, which determines when it begins its nocturnal activity, (2) the major expressed OR56d receptor may be correlated to suitable humidity sensing, and (3) Luc1-type luciferase is responsible for Abs. cerata’s luminescent spectrum. Key words: Nanopore, firefly, Lampyridae, Abscondita, Genome, Vision, Olfaction, Luciferase. BACKGROUND Abscondita cerata (Fig. 1) is an abundant Abscondita species (Ballantyne et al. 2013 2019) that ranges widely from sea level to 1,500 m above sea level and is endemic to the country side of Taiwan island (Ohba and Yang 2003). Its mating season is mainly from March to May (Wu et al. 2010). Abs. cerata is most active during the nighttime around 18:30–23:00 (Goh et al. 2022a; Ohba and Yang 2003). Several cohabitated species also can be found during its mating season, such as the genus Luciola, Aquatica, Citation: Wang TY, Goh KS, Wang LJ, Lee CM, Wu LL, Wang FY, Wu YW. 2023. Long-Read Genome Sequencing of Abscondita cerata (Coleoptera: Lampyridae), the endemic firefly of Taiwan. Zool Stud 62:25. doi:10.6620/ZS.2023.62-25. Zoological Studies 62:25 (2023) doi:10.6620/ZS.2023.62-25 1 © 2023 Academia Sinica, Taiwan Abscondita, Pyrocoelia, Curtos (Goh et al. 2022ab; Jeng et al. 1998 1999; Ohba and Yang 2003). Among the cohabitated species in Taiwan, male Abs. cerata produce the brightest flashes (luminescence emissions λmax = 563.6 ± 0.3 nm), measuring up to 14 lux (Goh et al. 2022b). Abs. cerata starts its nocturnal activity when the environmental light intensity is lower than 6.49 lux (Goh et al. 2022b). Comparative morphological characters provided evidence to support the transfer of six Luciola species (including Abs. cerata) into Abscondita (Ballantyne et al. 2013). Two new species of Abscondita and Luc. pallescens were also included in a later study (Ballantyne et al. 2019). In addition, the mitochondrial cytochrome c subunit I (COI) barcode with 3.27%–12.3% “intraspecific” variations revealed the uniqueness of Taiwanese fireflies, such as Abs. chinensis, Aquatica ficta, Luc. curtithorax, Luc. filiformis and Curtos costipennis (Goh et al. 2022b). Thus, more genome comparison could help to solve the phylogenetic placements of these cryptic Asian fireflies. A total evidence phylogenetic approach revealed the evolution of adult luminescence in fireflies (Martin et al. 2017). There are two molecular bases for firefly luminescence communication: one gene family (opsins) controls the detection of the signal, and another gene family (luciferase) controls the production of the signal (Martin et al. 2015). The major insect opsin class (SW) that typically confers sensitivity to “bluelight” was lost before the origin of modern beetles, whereas duplications within the UV opsin class have likely overcome the loss of blue sensitivity (Sharkey et al. 2017). Thus, the color vision of fireflies are mainly controlled by two family genes, the 480–600 nm longwavelength-sensitive (LWS) opsin and the 300–400 nm ultraviolet-sensitive (UVS) opsin (Owens and Lewis 2018). Previous studies also supported the idea that shortand mid-wavelength artificial light may influence the flash signal of fireflies, such as Aq. ficta, which can sense these light wavelengths (Owens et al. 2018). On the other hand, two paralogous luciferases (Luc1 and Luc2) were identified in fireflies, in which Luc1-type luciferase was responsible for producing the “yellowish” luminescence of adult firefly lanterns while Luc2-type luciferase was used to produce the distinct dim, greenish glow of eggs and the whole pupa (Bessho-Uehara et al. 2017; Oba et al. 2013). High-throughput genome sequencing is an important milestone for enhancing the quality and speed of traditional genome research, in which both nextand third-generation sequencing platforms contribute to accelerating the study of draft genome sequences of non-model organisms. This study demonstrates a hybrid genome assembly of the firefly Abs. cerata, which is an endemic species of Abscondita in Lampyridae, by using low coverage of long-read sequences followed by error correction using reasonable coverage of shortread sequences. The phylogenetic and comparative genomic analyses provided insights into the evolution of luciferase, opsin genes and olfactory receptor (OR) genes in fireflies. The draft genome also serves as a valuable resource for the molecular basis of firefly luminescence and ecology conservation. Fig. 1. Abscondita cerata male photographed by Shih-Chieh Huang in Chienshih Township, Hsinchu County, Taiwan. page 2 of 14 Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan MATERIALS AND METHODS Specimen Seventeen adult male specimens were collected from the habitat in Nankang, Taipei City, Taiwan (25°01'40.4"N 121°38'02.6"E) in April 2019. After they were collected from the field, each living specimen was transferred and stored in a 50 mL sampling tube with moist paper tissues in the laboratory, under a 12L:12D photoperiod at 25°C. To obtain the head and lantern tissues used for extracting RNA and genomic DNA, each specimen was sacrificed as previously described (Goh and Li 2011). The isolated head or lantern was immediately preserved via Invitrogen™ RNAlater™ Stabilization Solution (Cat No. AM7021) in a -80°C freezer for further RNA extraction. The remaining thorax and body tissues were freshly preserved in a -20°C freezer for further genomic DNA extraction. Genome size estimated by propidium iodide (PI) DNA staining with flow-cytometer The fresh head tissues were dissected in 2 mL filtered Galbraith buffer (dissolve 4.26 g MgCl2, 8.84 g sodium citrate, 4.2 g 3-[N-morpholino] propane sulfonic acid (“MOPS”), 1 mL Triton X-100, and 1 mg boiled ribonuclease A into 1 L of ddH2O with a final pH of 7.2), and the debris was filtered via 5-mL round bottom tubes (Falcon Cat# 352235). The filtered cell solution was centrifuged at a low rate of 1500–3000 rpm at 4°C for 10 minutes. After removing the supernatant, remaining cells were re-suspended in 2 mL of 2% paraformaldehyde for 10–30 minutes then centrifuged again and resuspended with 1 mL of 1X PBS (Phosphate Buffered Saline) for PI staining (Orgogozo and Rockman 2011). 10,000 events of stained cells were analyzed on BD FACS Verse flow cytometer. During analysis, aggregates were removed by applying a gate for doublet elimination. The percentage of cells in G0/ G1, S, G2/M and Sub G1 phases of the cell cycle was determined using BD FACSDiva 8.0.2 software (Becton, Dickinson and Company BD Biosciences, USA). The genome size estimation was conducted following the protocol as previously described (Orgogozo and Rockman 2011) and normalized using a Danio rerio AB strain (genome size = 1.412 Gb). Whole-genome sequencing The genomic DNA of thorax tissue was extracted via the ZR Tissue & Insect DNA MicroPrep™ kit (D6015) following the supplier’s instructions. The extracted genomic DNA was checked in 0.8–1% agarose gel and concentrated by the size selection of KAPA Pure Beads (Cat# KR1245). For long-read sequencing, the concentrated genomic DNA (1–3 μg) was endrepaired and ligated with Nanopore sequencing adaptors (AMX in Nanopore Ligation Sequencing Kit, SQKLSK109) via a KAPA Hyper Prep Kit (Cat#KR0961, Kapa Biosystems, Wilmington, MA), following the manufacturer’s instructions. The genomic DNA library was premixed with the LB and SQB of the Nanopore Ligation Sequencing Kit (SQK-LSK109), loaded in a Flow Cell (R9.4.1; FLO-MIN106), and sequenced by MinION devices for 48–72 h in laboratory. The BGISEQ-500RS short reads data were referred to in a previous study (Wang et al. 2021) and can be accessed at NCBI SRA (accession no. SRR13626958). Transcriptomic sequencing The RNA of each tissue was grinded in Invitrogen™ TRIzol™ Reagents (Cat. No. 15596026) and then extracted using a Qiagen RNeasy Mini Kit (Cat. No. 74104) with in-silico DNase I incubation, following the manufacturer’s instructions. The library of purified mRNA was constructed using a TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA), following the manufacturer’s instructions. The constructed libraries were sequenced on an Illumina NovaSeq 6000 platform with 150 bp paired-end reads generated by Genomics, BioSci & Tech Co., New Taipei City, Taiwan. Genome Assembly and Annotation Nanopore reads were trimmed using Porechop v0.2.4 (Wick et al. 2017) and assembled using canu v2.2 (Koren et al. 2017) (parameter: -nanopore-raw). The assembled contigs were corrected by the following procedure. The assembly was first corrected by Racon v1.4.3 (Vaser et al. 2017) for four rounds (parameters: -m 8 -x -6 -g -8 -w 500; alignment was made using minimap2 v2.22-r1101 (Li 2018)) followed by one round of Medaka v1.4.4 correction (using medaka_ consensus) (https://github.com/nanoporetech/medaka). The BGI paired-end reads (Wang et al. 2021) were trimmed using SOAPnuke v2.0.7 (Chen et al. 2018) with parameters “-l 10 -q 0.5 -n 0.1 -M 2 --adaMR 0.5 -f AAGTCGGAGGCCAAGCGGTCTTAGGAAGA CAA -r AAGTCGGATCGTAGCCATGTCGTTCTG TGAGCCAAGGAGTTG” and mapped against the racon/medaka-corrected contigs using bowtie2 v2.2.3 (Langmead and Salzberg 2012). After sorting and indexing the mapped reads into BAM format using samtools v1.9-45-ge415187 (Li et al. 2009), the mapped reads alignment was then used to correct the medakapage 3 of 14Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan corrected contigs using Pilon v1.23 (Walker et al. 2014) with default parameters. Since the genome was assembled using long reads and polished using both long reads and short reads, no additional scaffolding procedure was conducted explicitly, as the scaffolding information was already embedded in the long reads assembly process. The genome completeness was evaluated using BUSCO v4.0.6 (Manni et al. 2021) with insecta_odb10 library. The genes were annotated using maker v2.31.10 (Holt and Yandell 2011). The annotation steps were briefly introduced as follows. We first identified the repetitive elements in the assembled genome using RepeatModeler v1.0.11 (http://www.repeatmasker.org/ RepeatModeler/), re-annotated the “Unknown” elements into known repeat types using TEClass v2.1.3 (Abrusan et al. 2009), and quantified the genomic repeat content using RepeatMasker v4.1.1 (https://www.repeatmasker. org/RepeatMasker/). The repeat library was then specified in the maker configuration file along with the transcriptome assembled using Trinity v2.9.0 (Grabherr et al. 2011). Four runs of maker were then proceeded, in which the first run was executed with transcriptome data as a gene prediction guide, the second and third runs were conducted with trained gene models using SNAP v2.48.3 (Korf 2004), and the fourth run was performed with updated Augustus v3.3.2 gene profiles (Stanke et al. 2006). Gene functions were predicted by mapping the predicted genes against RefSeq (Pruitt et al. 2012) (with -evalue 1e-5 cutoff) and extracting the functional annotation of the best-matched entities for the predicted genes. “Hypothetical genes” were assigned to genes that could not be mapped to any RefSeq sequences. Genome contaminations were also checked in this step, in which only contigs with the majority (> 50%) of genes mapped to Arthropoda phylum were kept for further analysis. Data Records Raw sequencing libraries and genome and transcriptome assemblies were deposited into NCBI SRA as part of the BioProject PRJNA699421. The nanopore reads were available under accession number SRR18531391. The BGISEQ-500RS short reads data were referred to the associated SRA and Bio-Sample numbers are SRR13626958 and SAMN17776035, respectively (Wang et al. 2021). The genome assembly with gene and transcript annotations was deposited into GenBank under the accession number JALBCR000000000. The Illumina RNA sequencing reads were available under accession number SRR20723813 (lantern) and SRR20723814 (head). Nuclear phylogeny Five firefly-related genomes (see Table 1) were downloaded from NCBI and other data sources. The genome of Tribolium castaneum (Family: Tenebrionidae) and Ignelater luminous (Family: Elateridae) were included as outgroups. Single copy genes were identified using orthoMCL v2.0.9 (with mcl --abc -I 1.5) (Li et al. 2003). The single copy genes were then aligned separately using Muscle v3.8.31 (Edgar 2004); poorly-aligned regions were trimmed using Gblocks v0.91b (Castresana 2000). The alignments of the single copy genes were concatenated using a custom Perl script, and the phylogenetic tree was built from the alignment using FastTree v2.1.9 (Price et al. 2010) with default parameters and visualized using MEGA X (Kumar et al. 2018). The synonymous and non-synonymous trees were Table 1. Species and data sources used for building the phylogenetic tree. The putative LWS gene candidates from each genome were identified and reannotated in the table Species Data source Accession no. or website LWS (hypothetic protein) Reference Lamprigera yunnana NCBI GCA_013368075.1 FQA39_LY08143_1402397-1401056 (Zhang et al. 2020) Photinus pyralis Fireflybase http://www.fireflybase.org/ Ppyr1.3_LG5_10739432-10740780 (Fallon et al. 2018) Aquatica lateralis Fireflybase http://www.fireflybase.org/ Alat1.3_scaffold_703_90391-91758 (Fallon et al. 2018) Abscondita terminalis NCBI GCA_013368085.1 FQR65_LT10610_1659332-1660681 (Zhang et al. 2020) Pyrocoelia pectoralis Gigadb http://gigadb.org/dataset/100376 maker-scaffold382-augustus-gene1.421-mRNA (Fu et al. 2017) Outgroup: Ignelater luminosus Fireflybase http://www.fireflybase.org/ Ilumi1.2_Scaffold13785_69352-65621 (Fallon et al. 2018) Tribolium castaneum NCBI GCA_000002335.3 NA (Tribolium Genome Sequencing et al. 2008) page 4 of 14Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan constructed using the codeml program of PAML v4.10.6 package (Yang 2007). The dN/dS was also estimated using PAML. Comparative transcriptome of adult male head and lantern To reveal tissue-specific expressed genes, we compared transcriptomes between the adult male head and lantern. The transcriptome short reads were trimmed using Trimmomatic v0.39 (Bolger et al. 2014) with the following parameters: ILLUMINACLIP:TruSeq3PE.fa:2:30:10 LEADING:12 TRAILING:12 SLIDINGWINDOW:4:15 MINLEN:36. The trimmed reads were then mapped against the assembled genome using bowtie2 v2.2.3 (Langmead and Salzberg 2012). The differential gene expression analysis was carried out by initially quantifying gene counts from the mapped reads using featureCounts v1.6.1 (Liao et al. 2014) followed by differential analysis using NOISeq v2.40.0 (Robinson and Oshlack 2010; Tarazona et al. 2015; Zhao et al. 2021) in terms of RPKM (reads per kilobases per million reads). Differentially-expressed genes were defined as those with at least 2 fold-of-changes. The expressed OR, luciferase and opsin genes in adult males were manually identified from gene annotations. Since some genes such as the LWS opsin genes were incomplete in the assembled genome, an additional transcriptome quantification was conducted based on transcriptome assembly. The Trinity-assembled transcriptome was first clustered at 95% amino acid identity using CD-HIT v4.7 (Li and Godzik 2006) followed by RNA-Seq reads mapping using bowtie2 2.2.3 (Langmead and Salzberg 2012). RSEM v1.2.28 (Li and Dewey 2011) was used to quantify the expression of genes predicted by Transdecoder v5.3.0 (https://github. com/TransDecoder/TransDecoder) in terms of FPKM (fragments per kilo-bases per million reads). Photic environment measurement The changes of the environmental spectrum before and after firefly activities were measured using an Ocean optic USB2000+ spectrometer (Ocean Optics Inc., Dunedin, Florida, USA). The measurements were performed in an open place without the shade of tall vegetation located at the studied habitat in Nangang. During the measurement, the spectrometer probe was directly pointing towards the sky to collect the spectral range within (350–1,000 nm). Five readings were collected for each of the following time points: 1) before the fireflies started activity, 2) when the fireflies started blinking and 3) when the fireflies started flying. LWS gene amplicon validation Crude DNA was extracted from thoracic muscles via the ZR Tissue & Insect DNA MicroPrep™ kit (D6015). The specific primers (AcGR_fwd: 5'-ATATCAGAATGTCGGTaTTgGGT-3' and AcGR_ rev: 5'-TTATGCAGTTGCTTTTTCTTCTGA-3') were designed based on available LWS genes of fireflies to amplify a 1348-bp segment. Polymerase chain reactions (PCRs) in 25-μL volumes were performed with a dNTP concentration of 200 μM and a primer concentration of 0.3 μM, with 50 ng of genomic DNA, one unit of TaKaRa Taq™ DNA Polymerase, and the buffer supplied by the manufacturer. The PCR was run for 30 cycles under the following conditions: denaturation at 95°C for 20 s, annealing at 55–58°C for 1 min, extension at 72°C for 2 min, and a final extension at 72°C for 5 min. The product mixture was used as a template for DNA sequencing (Mission Biotech CO., LTD, Taipei, Taiwan). The re-sequenced LWS gene was deposited into GenBank under accession numbers OQ225001. RESULTS De novo genome assembly Table 2 summarized the genome assembly of Abs. cerata. Both long and short reads were sequenced from the Abs. cerata tissue samples, in which 2,121,726 reads with genome coverage of 10.38X were yielded from one FLO-MIN106 flowcell while 47,335,845 short reads (total size 9,467,169,000 bps) with genome coverage of 9.79X were obtained from BGISEQ-500RS. The assembly and polishing of the genome using canu v2.2 yielded 4,855 contigs (N50 = 325,269; N90 = 99,159; L50 = 907; L90 = 2,979) after the decontamination of non-firefly contigs, such as microbial genomic sequences. The assembly genome size was 967.14 Mb, with a longest contig of 2.115 Mb. 47.11% of the genome sequences consisted of repeat elements, of which DNA transposons accounted for more than half (26.83% of the genomic sequences, or 56.9% among all identified repeat). Alternatively, the genome size could be estimated via flow cytometer (see MATERIALS AND METHODS). Table 3 shows that the genome size was about 715 Mb. The size variation between PI staining and de novo genome assembly might be due to insufficient staining of nuclear DNA or copy number variation of repeat elements. page 5 of 14Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan Genome annotation and BUSCO analysis 127,704,450 raw reads were generated through RNA sequencing, and 126,260,884 (98.87%) clean reads were obtained after removing low quality reads. With RNA-seq based annotation to predict proteincoding genes, 55,206 protein-coding genes were predicted, and 20,862 genes could be annotated (Table S1). The genome of Abs. cerata is about 85% complete (BUSCO4 reported C:81.7% [S:80.5%, D:1.2%], F:3.1%, M:15.2%, n:1367 for Insecta lineage). Also, 79% of BUSCO gene completeness could be predicted by the transcriptome (C:73.9% [S:70.9%, D:3.0%], F:5.1%, M:21.0%, n:1367). Differentially expressed genes between adult male head and lantern We identified 12,747 out of the 20,863 annotated genes that were expressed in head and/or lantern tissues (Table S2). Among the 12,747 expressed genes, 2,076 and 2,190 of the differentially expressed genes were above two fold-of-changes in head and lantern tissues, respectively (Table S3). We wanted to know if the opsin and luciferase genes are related to nocturnal activity and mating behavior. There is one UVS opsin gene and one LWS opsin gene within the Abs. cerata genome. The transcriptome of the head revealed that the UVS opsin gene is the eighth highly expressed (FPKM = 97,932 in head; 12 in lantern) in male adults (Table 4a). However, the LWS opsin gene is incomplete in the genome draft due to the low coverage reads assembly and/or sequence error, with only a partial 5'-end fragment (377 bps) revealed in the 3'-end of tig00652085_segment0_pilon (84,107 kb) (Fig. 2). We then designed the specific primers to amplify the LWS opsin gene amplicon using genomic DNA. The complete LWS opsin gene amplicon (1348 bps) could be re-sequenced by Sanger sequencing. The results showed 5 exons (129–378 bps) with 4 short introns (46–55 bps). The short introns with A/T homopolymers may be the noise in the genome assembly. Since we cannot calculate those incomplete gene expressions from genome, we re-calculated the UVS and LWS opsin gene expressions via de novo transcriptomic assembly (Table 4b), in which the FPKM for the LWS and UVS opsin gene expression were 7,617 and 464 in head and 0 for both in lantern, respectively. On the other hand, two ultra-high expressed genes could be found in the transcriptome of the lantern (Table 4a). The FPKM value of Luciferase 1 (Luc1; snap_masked-tig00000382_segment0_pilon-processedgene-2.24) is 1,037,989 while the other FPKM value of the Cytochrome P450 4G15-like gene (CYP4G15; maker-tig00064854_segment0_pilon-snap-gene-0.23) is 805,540 in the lantern. The transcriptome of the lantern revealed highly expressed Luciferase 1 (FPKM = 2 in head; 1,037,989 in lantern) and less expressed Luciferase 2 (Luc2; maker-tig00008619_segment0_ pilon-augustus-gene-2.43; FPKM = 18 in head; 77 in lantern) genes in male adults. The top two ultra-high expressed genes in the transcriptome of the head were arrestin-like and a small heat-shock protein (Table S3). The FPKM value of arrestin (maker-tig00014644_ segment0_pilon-augustus-gene-1.39) is 929,529 Table 2. Genome assembly summary of Abscondita cerata Abscondita cerata Assembled genome size (bp) 967,144,706 Genome coverage (fold) 10.38x Total no of contigs 4,855 Max contig length (bp) 2,115,774 Contig N50 (bp) 325,269 Contig L50 907 Contig N90 (bp) 99,159 Contig L90 2,979 Number of genes predicted 55,206 Mean exon length (bp) 268.16 Average no of exons/gene 3.82 Insecta BUSCO gene number 1,367 % complete BUSCOs* 81.7% % complete and single-copy BUSCOs 80.5% % complete and duplicated BUSCOs 1.2% % fragmented BUSCOs 3.1% % missing BUSCOs 15.2% * BUSCO version is 4.0.4; The lineage dataset is eukaryota_odb10 (creation date: 2019-11-20, number of species: 70, number of BUSCOs: 255). Table 3. Genome size determination using flow cytometer of propidium iodide (PI)-stained nuclei Species Genome size (Mb) Mean_PI (2N) SD_PI (2N) Danio rerio* 1,412 102.991 7.999 Abscondita cerata (BM1) 717 ± 64 52.272 4.647 Abscondita cerata (BM3) 712 ± 71 51.938 5.158 * zebrafish AB strain and the genome size is referred to 1,412 Mb from NCBI. Standard deviation (SD). page 6 of 14Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan Table 4. Selected gene candidates expressed in head or lantern tissues. The expressed olfactory receptor, luciferase and opsin genes in adult males were specifically addressed in the discussion (a) Predicted genes from genome: Gene id Putative gene Head Lantern Fold_of Pfam Protein family (RPKM) (RPKM) Change snap_masked-tig00064992_segment0_pilon-processed-gene-1.66 UVS opsin 97,932 21 12.2 PF00001 7 transmembrane receptor tig00652085_segment0_pilon (gDNA)* LWS opsin NA NA NA NA NA snap_masked-tig00000382_segment0_pilon-processed-gene-2.24 Luciferase 1 21,037,989 -19.02 PF00501 AMP-binding enzyme maker-tig00008619_segment0_pilon-augustus-gene-2.43 Luciferase 2 18 77 -2.14 PF00501 AMP-binding enzyme maker-tig00064854_segment0_pilon-snap-gene-0.23 CYP4G15 7.78 805,540 -16.659 PF00501 Cytochrome P450 maker-tig00651154_segment0_pilon-snap-gene-0.67 OR56d 18,671 213.13 PF01395 PBP/GOBP family snap_masked-tig00010855_segment0_pilon-processed-gene-0.13 OR19d 7,857 1,193 2.72 PF01395 PBP/GOBP family maker-tig00651459_segment0_pilon-snap-gene-0.34 OR56d 7,784 4 10.87 PF01395 PBP/GOBP family maker-tig00010855_segment0_pilon-augustus-gene-0.4 OR19d 4,225 353 3.58 PF01395 PBP/GOBP family maker-tig00064231_segment0_pilon-augustus-gene-2.19 OR83a 4,073 407 3.32 PF01395 PBP/GOBP family maker-tig00010855_segment0_pilon-snap-gene-1.0 OR19d 3,146 346 3.18 PF01395 PBP/GOBP family maker-tig00005308_segment0_pilon-augustus-gene-2.8 OR19d 2,433 0 12.29 PF01395 PBP/GOBP family genemark-tig00004518_segment0_pilon-processed-gene-4.15 OR99b 1,408 1,011 0.48 PF01395 PBP/GOBP family maker-tig00005221_segment0_pilon-augustus-gene-3.75 OR59a 1,328 136 3.29 maker-tig00065369_segment0_pilon-augustus-gene-0.64 OR83a 1,223 19 6.03 PF01395 PBP/GOBP family maker-tig00012100_segment0_pilon-augustus-gene-0.3 OR56d 1,197 0 11.26 PF01395 PBP/GOBP family maker-tig00652072_segment0_pilon-augustus-gene-0.6 OR19d 965 0 10.95 PF01395 PBP/GOBP family * only 377 bp fragment could be found in genome sequence. (b) Predicted genes from transcriptome: Gene id Putative gene Head Lantern Pfam Protein family (count, FPKM) (count, FPKM) TRINITY_DN10655_c0_g1_i2 UVS opsin (17,753, 464) (0, 0) PF00001 7 transmembrane receptor TRINITY_DN15410_c3_g7_i1 LWS opsin (3,78,270, 7,617) (0, 0) PF10324.8 Serpentine type 7TM GPCR chemoreceptor TRINITY_DN16196_c1_g1_i1 Luciferase 1 (1, 0.02) (465,842, 7,269) PF00501.27 AMP-binding enzyme TRINITY_DN51863_c0_g1_i1 Luciferase 2 (860, 9.40) (1,207, 10.6) PF00501.27 AMP-binding enzyme Fig. 2. Long wavelength-sensitive (LWS) and ultraviolet-sensitive (UVS) opsin gene structure. The opsin genes were amplified by degenerated primers and validated via the sanger sequencing. The exons and introns were predicted by the NCBI orffinder and supported by transcriptomic reads. page 7 of 14Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan while the other FPKM value of the small Heat shock protein-like gene (Hsp20; snap_masked-tig00003087_ segment0_pilon-processed-gene-3.6) is 866,534. In addition, there are 111 putative OR genes, classified into 27 OR gene families, within the Abs. cerata genome (Table S1). The transcriptome of the head revealed that 36 OR genes were expressed in male adults (Table 4a). The 36 OR genes from 15 gene families were predicted and belonged to two Pfam domains (PF01395: PBP/GOBP family and PF02949: 7tm odorant receptor). Six 7tm odorant receptors were less expressed (FPKM = 2–540 in head). For the rest of the 30 PBP/GOBP family genes, only 11 genes were specifically expressed in the head (higher FPKM > 1000 and/or FPKM = 0 in lantern), which were one OR59a, three OR56d, four OR19d, one OR99b, two OR83a family genes (Table 4a). “maker-tig00651154_ segment0_pilon-snap-gene-0.67”, an OR56d family gene, was the highest expressed OR gene, with FPKM = 18,671. The other one OR56d and two OR19d genes also showed a higher expression level than other OR genes. Molecular phylogeny of single copy nuclear genes We reconstructed the molecular phylogeny of 6 species of Lampyridae (see Table 1) based on 1,792 orthologous single copy genes (12,867,000 amino acids after alignment and Gblocks trimming) from the genomic data (Fig. 3a). Bootstrap values shown at the branch nodes indicate every clade within Lampyridae is well supported (100%). As shown in figure 3a, the family Lampyridae is a monophyletic group. The monophyly of the subfamily Luciolinae contains Aquatica lateralis and (Abs. cerata + Abs. terminalis) while Lamprigera yunnana is sister to the clade of (Aq. lateralis + (Abs. cerata + Abs. terminalis)). The Photinus pyralis and Pyrocoelia pectoralis form another monophyly of the subfamily Lampryriinae. Among studied species, Abs. cerata is the endemic island species while others are all continental species. The phylogenetic tree suggested that the Abs. cerata has a slightly faster rate in proteins than that of relatively close species Abs. terminalis. By building both synonymous and non-synonymous trees (Fig. 3b and 3c) and estimated the dN/dS value, we estimated that the dN/dS was 0.0922, suggesting negative selection. The comparison between both the synonymous and non-synonymous tree branch lengths of Abs. cerata and Abs. terminalis also shows that the mutations in Abs. cerata has ~1.42 fold discrepancy than those in Abs. terminalis. Photic environment Figure 4 shows how environmental short wavelength light (< 500 nm) dramatically decreases during the twilight period, especially before and after the starting time of firefly nocturnal activity. The Abs. cerata start to flash or fly 5–10 minutes later when the short wavelength light disappears. In addition, the mid wavelength light (~400–600 nm) could be still detectable at this moment. Fig. 3. Phylogenetic tree of fireflies (Lampyridae) constructed by 1,792 single copy orthologous genes. (a) protein tree, (b) synonymous tree and (c) non-synonymous tree between two Abscondita species. The species Tribolium castaneum and Ignelater luminous were used as the outgroups for the protein tree; the Aquatica lateralis species was used as the species for the synonymous and non-synonymous tree. The assembled genome size is taken from table 5. page 8 of 14Zoological Studies 62:25 (2023) © 2023 Academia Sinica, Taiwan DISCUSSION This study demonstrated a hybrid genome assembly of an endemic firefly (Abs. cerata in Lampyridae) using low coverage of long-read sequences and error correcting by reasonable coverage of short-read sequences. Previous study indicated the large variations of firefly COI genes between Taiwan and adjacent regions (Goh et al. 2022b). Thus, the endemic Abs. cerata is a suitable genome reference for the molecular basis of firefly luminescence and ecology conservation. Genome evolution of firefly Previous studies (Liu et al. 2017; Lower et al. 2017) revealed the genome size of fireflies ranged from 410.76 to 2,154.37 Mb, which contains repeats from 10.3% to 56.6% (Lower et al. 2017). Figure 5 suggested the genome size variation may be due to different proportion of repeats (Zhang et al. 2020) and/or (tandem) gene duplications, such as acyl-CoA synthetase (ACS), luciferase, etc (Oba and Schultz 2022). A previous study indicated the ACS gene underwent a tandem duplication on a single chromosome while the neofunctionalized luciferase evolved from adjacent ACS and duplicated to another chromosome later (Oba and Schultz 2022). Since the Abs. cerata genome is still fragmented, these genes remained to separate into different contigs (data not shown). Further chromosomal-level genome will give more insight on the endemic genome evolution. Consistent phylogeny between nuclear and mitochondrial genes The nuclear phylogeny, constructed using 1,792 single copy genes, agreed with the other phylogenetic trees using mitochondrial genes (Chen et al. 2019; Goh et al. 2022b; Wang et al. 2021) and 436 nuclear genes (Martin et al. 2019). Figure 3a revealed the monophyly Fig. 5. Genome size and total repeats length of studied fireflies. Fig. 4. Photic environment before and after the starting time of firefly nocturnal activity. After twilight begins, there is a 5–10 minute time interval between the start of firefly flash and fly conditions on two independent days (5/17 and 5/18). page 9 of 14Zoological Studies 62:25 (2023)