Firefly genomes illuminate parallel origins of bioluminescence in beetles

Fireflies and their luminous courtships have inspired centuries of scientific study. Today firefly luciferase is widely used in biotechnology, but the evolutionary origin of bioluminescence within beetles remains unclear. To shed light on this long-standing question, we sequenced the genomes of two firefly species that diverged over 100 million-years-ago: the North American Photinus pyralis and Japanese Aquatica lateralis. To compare bioluminescent origins, we also sequenced the genome of a related click beetle, the Caribbean Ignelater luminosus, with bioluminescent biochemistry near-identical to fireflies, but anatomically unique light organs, suggesting the intriguing hypothesis of parallel gains of bioluminescence. Our analyses support independent gains of bioluminescence in fireflies and click beetles, and provide new insights into the genes, chemical defenses, and symbionts that evolved alongside their luminous lifestyle.


Introduction
Fireflies (Coleoptera: Lampyridae) represent the best-studied case of bioluminescence. The coded language of their luminous courtship displays ( Figure 1A; Video 1) has been long studied for its role in mate recognition (Lloyd, 1966;Lewis and Cratsley, 2008;Stanger-Hall and Lloyd, 2015), while non-adult bioluminescence is likely a warning signal of their unpalatable chemical defenses (De Cock and Matthysen, 1999), such as the cardiotoxic lucibufagins of Photinus fireflies (Meinwald et al., 1979). The biochemical understanding of firefly luminescence: an ATP, Mg 2+ , and O 2 -dependent luciferase-mediated oxidation of the substrate luciferin (Shimomura, 2012), along with the cloning of the luciferase gene (de Wet et al., 1985;Ow et al., 1986), led to the widespread use of luciferase as a reporter with unique applications in biomedical research and industry (Fraga, 2008). With >2000 species globally, fireflies are undoubtedly the most culturally appreciated bioluminescent group, yet there are at least three other beetle families with bioluminescent species: click beetles (Elateridae), American railroad worms (Phengodidae) and Asian starworms (Rhagophthalmidae) (Martin et al., 2017). These four closely related families (superfamily Elateroidea) have homologous luciferases and structurally identical luciferins (Shimomura, 2012), implying a single origin of beetle bioluminescence. However, as Darwin recognized in his 'Difficulties on Theory' (Darwin, 1872), the eLife digest Glowing fireflies dancing in the dark are one of the most enchanting sights of a warm summer night. Their light signals are 'love messages' that help the insects find a mate -yet, they also warn a potential predator that these beetles have powerful chemical defenses. The light comes from a specialized organ of the firefly where a small molecule, luciferin, is broken down by the enzyme luciferase.
Fireflies are an ancient group, with the common ancestor of the two main lineages originating over 100 million years ago. But fireflies are not the only insects that produce light: certain click beetles are also bioluminescent.
Fireflies and click beetles are closely related, and they both use identical luciferin and similar luciferases to create light. This would suggest that bioluminescence was already present in the common ancestor of the two families. However, the specialized organs in which the chemical reactions take place are entirely different, which would indicate that the ability to produce light arose independently in each group.
Here, Fallon, Lower et al. try to resolve this discrepancy and to find out how many times bioluminescence evolved in beetles. This required using cutting-edge DNA sequencing to carefully piece together the genomes of two species of fireflies (Photinus pyralis and Aquatica lateralis) and one species of click beetle (Ignelater luminosus). The genetic analysis revealed that, in all species, the genes for luciferases were very similar to the genetic sequences around them, which code for proteins that break down fat. This indicates that the ancestral luciferase arose from one of these metabolic genes getting duplicated, and then one of the copies evolving a new role.
However, the genes for luciferase were very different between the fireflies and the click beetles. Further analyses suggested that bioluminescence evolved at least twice: once in an ancestor of fireflies, and once in the ancestor of the bioluminescent click beetles.
More results came from the reconstituted genomes. For example, Fallon, Lower et al. identified the genes 'turned on' in the bioluminescent organ of the fireflies. This made it possible to list genes that may be involved in creating luciferin, and enable flies to grow brightly for long periods. In addition, the genetic information yielded sequences from bacteria that likely live inside firefly cells, and which may participate in the light-making process or the production of potent chemical defenses.
Better genetic knowledge of beetle bioluminescence could bring new advances for both insects and humans. It may help researchers find and design better light-emitting molecules useful to track and quantify proteins of interest in a cell. Ultimately, it would allow a detailed understanding of firefly populations around the world, which could contribute to firefly ecotourism and help to protect these glowing insects from increasing environmental threats.
light organs amongst the luminous beetle families are clearly distinct ( Figure 1B), implying independent origins. Thus, whether beetle bioluminescence is derived from a single or multiple origin(s) remains unresolved.
To address this long-standing question, we sequenced and analyzed the genomes of three bioluminescent beetle species. To represent the fireflies, we sequenced the widespread North American 'Big Dipper Firefly', P. pyralis ( Figure 1A,C) and the Japanese 'Heike-botaru' firefly Aquatica lateralis ( Figure 1B). P. pyralis was used in classic studies of firefly bioluminescent biochemistry (Bitler and McElroy, 1957) and the cloning of luciferase (de Wet et al., 1985), while A. lateralis, a species with specialized aquatic larvae, is one of the few fireflies that can be reliably cultured in the laboratory (Oba et al., 2013a). These two fireflies represent the two major firefly subfamilies, Lampyrinae and Luciolinae, which diverged from a common ancestor over 100 Mya ( Figure 1B) (Misof et al., 2014;  patrol flashes over a field in Homer Lake, Illinois. Females cue in on these species-specific flash patterns and respond with their own species-specific flash (Lloyd, 1966). Photo credit: Alex Wild. Inset: male and female P. pyralis in early stages of mating. Photo credit: Terry Priest. (B) Cladogram depicting the hypothetical phylogenetic relationship between P. pyralis and related bioluminescent and non-bioluminescent taxa with Tribolium castaneum and Drosophila melanogaster as outgroups. Numbers at nodes give approximate dates of divergence in millions of years ago (mya) (Misof et al., 2014;Mckenna et al., 2015). Right: Dorsal and ventral photos of adult male specimens. Note the well-developed ventral light organs on the true abdominal segments 6 and 7 of P. pyralis and A. lateralis. In contrast, the luminescent click beetle, I. luminosus, has paired dorsal light organs at the base of its prothorax (arrowhead) and a lantern on the anterior surface of the ventral abdomen (not visible). (C) Empirical range of P. pyralis in North America, extrapolated from 541 reported sightings (Appendix 1.2). Collection sites of individuals used for genome assembly are denoted with circles and location codes. Cross hatches represent areas which likely have P. pyralis, but were not sampled. Diagonal hashes represent Ontario, Canada. DOI: https://doi.org/10.7554/eLife.36495.003 Mckenna et al., 2015). To facilitate evolutionary comparisons, we also sequenced the 'Cucubano', Ignelater luminosus ( Figure 1B), a Caribbean bioluminescent click beetle, and member of the 'Pyrophorus' used by Raphaë l Dubois (1849Dubois ( -1929 to first establish the enzymatic basis of bioluminescence in the late 1800s (Dubois, 1885;Dubois, 1886). Comparative analyses of the genomes of these three species allowed us to reconstruct the origin(s) and evolution of beetle bioluminescence.

Results
Sequencing and assembly of firefly and click-beetle genomes Photinus pyralis adult males were collected from the Great Smoky Mountains National Park, USA (GSMNP) and Mercer Meadows New Jersey, USA (MMNJ) ( Figure 1C), and sequenced using shortinsert, mate-pair, Hi-C, and long-read Pacific Biosciences (PacBio) approaches (Appendix 4-table 1). These datasets were combined in a MaSuRCA (Zimin et al., 2013) hybrid genome assembly (Appendix 1.5). The Aquatica lateralis genome was derived from an ALL-PATHs (Butler et al., 2008) assembly of short insert and mate-pair reads from a single adult female from a laboratory-reared population, whose lineage, dubbed 'Ikeya-Y90', was first collected 25 years ago from a now extinct population in Yokohama, Japan (Appendix 2.5). A single Ignelater luminosus adult male, collected in Mayagü ez Puerto Rico, USA, was used to produce a high-coverage Supernova (Weisenfeld et al., 2017) linked-read draft genome (Appendix 3.5), which was further manually scaffolded using lowcoverage long-read Oxford Nanopore MinION sequencing (Appendix 3.5.4).
The gene completeness and contiguity statistics of our P. pyralis (Ppyr1.3) and A. lateralis (Alat1.3) genome assemblies are comparable to the genome of the model beetle Tribolium castaneum ( Figure 2F; Appendix 4.1). The I. luminosus genome assembly (Ilumi1.2) is less complete, but is comparable to other published insect genomes ( Figure 2F; Appendix 4.1). Protein-coding genesets for our study species were produced via an EvidenceModeler-mediated combination of homology alignments, ab initio predictions, and de novo and reference-guided RNA-seq assemblies followed by manual gene curation for gene families of interest (Appendix 1. 10; 2.8; 3.8). These coding gene annotation sets for P. pyralis, A. lateralis, and I. luminosus are comprised of 15,773, 14,285, and 27,557 genes containing 94.2%, 90.0%, and 91.8% of the Endopterygota Benchmarking Universal Single-Copy Orthologs (BUSCOs) (Simão et al., 2015), respectively. Protein clustering via predicted orthology indicated 77% of genes were found in orthogroups with at least one other species ( Figure 2E; Appendix 4-figure 1). We found the greatest orthogroup overlap between the P. pyralis and A. lateralis genesets, as expected given the more recent phylogenetic divergence of these species. Remaining redundancy in the P. pyralis assembly and annotation, as indicated by duplicates of the BUSCOs and the assembly size ( Figure 2F;  is likely due to the heterozygosity of the outbred input libraries (Appendix 1). The higher BUSCO completeness of the assemblies as compared to the genesets (Appendix 4-table 3), suggests that future manual curation efforts will lead to improved annotation completeness.
To enable the characterization of long-range genetic structure, we super-scaffolded the P. pyralis genome assembly into 11 pseudo-chromosomal linkage groups using a Hi-C proximity-ligation linkage approach (Figure 2A; Appendix 1.5.3). These linkage groups contain 95% of the assembly (448.8 Mbp). Linkage group LG3a corresponds to the X-chromosome based on expected adult XO male read coverage and gene content (Appendix 1.6.4.1) and its size (22.2 Mbp) is comparable to the expected X-chromosome size based on sex-specific genome size estimates using flow cytometry (~26 Mbp) (Lower et al., 2017). Homologs to T. castaneum X-chromosome genes were enriched on LG3a over every other linkage group, suggesting that the X-chromosomes of these distantly related beetles are homologous, and that their content has been reasonably conserved for >200 MY (Appendix 1.6.4.1) (Mckenna et al., 2015). We hypothesized that the P. pyralis orthologs of known bioluminescence genes, including the canonical luciferase Luc1 (de Wet et al., 1985) and the
(C) Genome schematic of P. pyralis mitochondrial genome (mtDNA). Like other firefly mtDNAs, it has a tandem repetitive unit (TRU) (Appendix 1.8). (D) mCG is enriched across gene bodies of P. pyralis and shows methylation levels that are at least two times higher than other holometabolous insects (Appendix 1.12). (E) Orthogroup (OGs) clustering analysis of genes with Orthofinder (Emms and Kelly, 2015) shows a high degree of overlap of the P. pyralis, A. lateralis, and I. luminosus genesets with the geneset of Tribolium castaneum. Numbers within curved brackets (colored by species) represent gene count from specific species within the shared orthogroups. Numbers with square brackets (black color) represent total gene count amongst shared orthogroups. OGs = orthogroups, *=Not fully filtered to single isoform per gene. See Appendix 4.2.1 for more detail. Intermediate scripts and species-specific overlaps are available as  4, 2.4, 3.4, 4.1). DOI: https://doi.org/10.7554/eLife.36495.005 The following source data is available for figure 2: Source data 1. Figure 2E. Orthogroup clustering analysis. DOI: https://doi.org/10.7554/eLife.36495.006 Source data 2. Excel file of Figure 2F AlatACS5 PpyrPACS2  P p y r M G S T P p y r A C S 3 P p y r A C S 6 P p y r P A C S 7 28,350 P P Y R _ 0 1 2 0 8 P P Y R _ 0 1 2 0 9 P p y r A C S 2 P p y r A C S 1 P p y r L u c 1 P p y r P A C S 2 P p y r P A C S 1 P p y r P A C S 3 P p y r P A C S 4 P p y r A C S 4 P p y r P A C S 5 Color gradients indicate the transcript per million (TPM) values of whole body in each sex/stage (grey to blue) and in the prothorax or abdominal lantern (grey to orange to green). Tree and annotation visualized using iTOL (Letunic and Bork, 2016). Prothorax and abdominal lantern expression values for I. luminosus are from whole prothorax plus head, and metathorax plus the two most anterior abdominal segments. Fluc = firefly luciferases, Eluc = elaterid luciferases, R/ PLuc = rhagophthalmid/phengodid luciferases. (Appendix 4.3.2) Gene tree, gene accession numbers, annotation, and expression values are available as The genomic context of firefly luciferase evolution Two luciferase paralogs have been previously described in fireflies (Oba et al., 2013a;Bessho-Uehara et al., 2017). P. pyralis Luc1 was the first firefly luciferase cloned (de Wet et al., 1985), and its direct orthologs have been widely identified from other fireflies (Oba, 2014). The luciferase paralog Luc2 was previously known only from a handful of Asian taxa, including A. lateralis Figure 3 continued scaffold 228, while 4 of the 13 P. pyralis PACS/ACS genes are close neighbors of PpyrLuc1 on LG1, with a further seven genes 2.4 Mbp and 39.1 Mbp away on the same linkage-group. Although the Luc1 loci in P. pyralis and A. lateralis are evidently derived from a common ancestor, the relative positions of the most closely related flanking PACS/ACS genes have diverged between the two species. IlumLuc was captured on a separate scaffold (Ilumi1.2_Scaffold13255) from its most most closely related PACSs (IlumPACS8, IlumPACS9) on Ilumi1.2_Scaffold9864, although three more distantly related PACS genes (IlumiPACS1, IlumiPACS2, IlumiPACS4) are co-localized with IlumLuc. In contrast, a different scaffold (Ilumi1.2_Scaffold9654) shows orthology to the firefly Luc1 locus. The full Ilumi1.2_Scaffold13255 was produced by a manual evidence-supported merge of two scaffolds (Appendix 3.5.4). Genes with a PTS1 are indicated by a dark outline, except for the genes with white interiors, which instead represent non-PACS/ACS genes without an identified homolog in the other scaffolds. Co-orthologous genes are labeled in the same color in the phylogenetic tree and are connected with corresponding color bands in synteny diagram. Genes and genomic regions are to scale (Scale bar = 25 Kbp). Gaps excluded from the figure are shown with dotted lines and are annotated with their length in square brackets. Scaffold ends are shown with rough black bars. MGST = Microsomal glutathione S-transferase, IMP = Inositol monophosphatase, PRNT = Polyribonucleotide nucleotidyltransferase. Figure produced with GenomeTools 'sketch' (v1.5.9) (Gremme et al., 2013). Figure production scripts available as Figure 3-source data 2. DOI: https://doi.org/10.7554/eLife.36495.008 The following source data is available for figure 3: Source data 1. Gene tree, gene accession numbers, annotation, and expression values for Figure 3C.     Figure 4. Parallel evolution of elaterid and firefly luciferase. (A) Ancestral state reconstruction recovers at least two gains of luciferase activity in bioluminescent beetles. Luciferase activity (top right figure key; black: luciferase activity, white: no luciferase activity, shaded: undetermined) was annotated on extant firefly luciferase homologs via literature review or inference via direct orthology. The ancestral states of luciferase activity within the putative ancestral nodes were then reconstructed with an unordered parsimony framework and a maximum likelihood (ML) framework (bottom left figure key; Appendix 4.3.3). Two gains ('G') of luciferase activity, annotated with black arrows and yellow stars, are hypothesized. These hypothesized gains occurred once in a gene within the common ancestor of fireflies, rhagophthalmid, and phengodid beetles, and once in a gene within the common ancestor of bioluminescent elaterid beetles. Scale bar is substitutions per site. Numbers adjacent to nodes represents node support. NEXUS and newick files available as Figure 4-source data 1 (B) Molecular adaptation analysis supports independent neofunctionalization of click beetle luciferase. We tested the molecular adaptation of elaterid luciferase using the adaptive branch-site REL test for episodic diversification (aBSREL) method (Smith et al., 2015) (Appendix 4.3.4). The branch leading to the common ancestor of elaterid luciferases (red star) was one of three branches (red and blue stars) recovered with significant (p<0.01) evidence of positive selection, with 35% of sites showing strong directional selection (w or max d N /d S = 3.98), which we interpret as signal of the initial neofunctionalization of elaterid ancestral luciferase (EAncLuc) from an ancestor without luciferase activity. As the selected branches with blue stars are red-shifted elaterid luciferases (Oba et al., 2010a;Stolz et al., 2003), they may represent the post-neofunctionalization selection of a few key sites via sexual selection of emission colors. Specific sites identified as under selection using Mixed Effect Model of Evolution (MEME) and Phylogenetic Analysis by Maximum Likelihood (PAML) methods are described in Appendix 4.3.4. The tree and results from the full adaptive model are shown. Branch length, with the exception of the PpyrLuc1 branch which was shortened, reflects the number of substitutions per site. Numbers adjacent to nodes represents node support. Figure was produced with iTOL (Letunic and Bork, 2016). Gene tree, metadata, and coding nucleotide multiple sequence alignment available as Figure 4-source data 2. DOI: https://doi.org /10.7554/eLife.36495.011 The following source data is available for figure 4: Source data 1. NEXUS and Newick files for luciferase ancestral state reconstruction in Figure 4A.   (Oba et al., 2013a;Bessho-Uehara et al., 2017). Previous investigations of these Asian taxa have shown that Luc1 is responsible for light production from the lanterns of adults, larvae, prepupae and pupae, whereas Luc2 is responsible for the dim glow of eggs, ovaries, prepupae and the whole pupal body . From our curated genesets (Appendix 1.10; 2.8), we unequivocally identified two firefly luciferases, Luc1 and Luc2, in both the P. pyralis and A. lateralis genomes. Our RNA-Seq data further show that in both P. pyralis and A. lateralis, Luc1 and Luc2 display expression patterns consistent with previous reports. While Luc1 is the sole luciferase expressed in the lanterns of both larvae and adults, regardless of sex, Luc2 is expressed in other tissues and stages, such as eggs ( Figure 3C). Notably, Luc2 expression is detected in RNA libraries derived from adult female bodies (without head or lantern), suggesting detection of ovary expression as described in previous studies . Together, these results support that since their divergence via gene duplication prior to the divergence of Lampyrinae and Luciolinae, Luc1 and Luc2 have established different, but conserved roles in bioluminescence throughout the firefly life cycle.
Firefly luciferase is hypothesized to be derived from an ancestral peroxisomal fatty acyl-CoA synthetase (PACS) ( Figure 3A) (Oba et al., 2003;Oba et al., 2006a). We found that, in both firefly species, Luc1 is genomically clustered with its closely related homologs, including PACSs and nonperoxisomal acyl-CoA synthetases (ACSs), enzymes which can be distinguished by the presence/ absence of a C-terminal peroxisomal-targeting-signal-1 (PTS1). We also found nearby microsomal glutathione S-transferase (MGST) family genes ( Figure 3D) that are directly orthologous between  The following source data is available for figure 5: both species, Genome-wide phylogenetic analysis of the luciferases, PACSs and ACSs genes indicates that Luc1 and Luc2 form two orthologous groups, and that the neighboring PACS and ACS genes near Luc1 form three major clades ( Figure 3C): Clade A, whose common ancestor and most extant members are ACSs, and Clades B and C whose common ancestors and most extant members are PACSs. Luc1 and Luc2 are highly conserved at the level of gene structure-both are composed of seven exons with completely conserved exon/intron boundaries (Appendix 4- figure 4; Appendix 4-figure 5), and most members of Clades A, B, and C also have seven exons. The exact syntenic and orthology relationships of the ACS and PACS genes adjacent to the Luc1 locus remains unclear, likely due to subsequent gene divergence and shuffling ( Figure 3C,D). Luc2 is located on a different linkage-group from Luc1 in P. pyralis and on a different scaffold from Luc1 in A. lateralis, consistent with the interpretation that Luc1 and Luc2 lie on different chromosomes in both firefly species. No PACS or ACS genes were found in the vicinity of Luc2 in either species. These data support that tandem gene duplication in a firefly ancestor gave rise to several ancestral PACS paralogs, one of which neofunctionalized in place to become the ancestral luciferase (AncLuc) ( Figure 3B). Prior to the divergence of the firefly subfamilies Lampyrinae and Luciolinae around 100 Mya (Appendix 4.3), this AncLuc duplicated, possibly via a long-range gene duplication event (e.g. transposon mobilization), and then subfunctionalized in its transcript expression pattern to give rise to Luc2, while the original AncLuc subfunctionalized in place to give rise to Luc1 ( Figure 3B). From the shared Luc gene clustering in both fireflies, we infer the structure of the pre Luc1/Luc2 duplication AncLuc locus contained one or more ACS genes (Clade A), one or more PACS genes (Clade B/C), and one or more MGST family genes ( Figure 3B).

Independent origins of firefly and click beetle luciferase
To resolve the number of origins of luciferase activity, and therefore bioluminescence, between fireflies and click beetles, we first identified the luciferase of I. luminosus luciferase (IlumLuc), and compared its genomic context to the luciferases of P. pyralis and A. lateralis ( Figure 3D). Unlike some other described bioluminescent Elateridae, which have separate luciferases expressed in the dorsal prothorax and ventral abdominal lanterns (Oba et al., 2010a), we identified only a single luciferase in the I. luminosus genome which was highly expressed in both of the lanterns ( Figure 3C; Appendix 3.8). The exon number and exon-intron splice junctions of IlumLuc are identical to those of firefly luciferases, but unlike the firefly luciferases which have short introns less than <100 bp long, IlumLuc has two long introns (Appendix 4-figure 4). We found several PACS genes in the I. luminosus genome which were related to IlumLuc and formed a clade (Clade D) specific to the Elateridae ( Figure 3C,D). IlumLuc lies on a 366 Kbp scaffold containing 18 other genes, including three related Clade D PACS genes (Scaffold 13255; Figure 3D; Figure 4); however, the Clade D genes that are most closely related to IlumLuc are found on a separate 650 Kbp scaffold (Scaffold 9864; Figure 3D). We infer that the IlumLuc locus is not orthologous to the extant firefly Luc1 locus, as IlumLuc is not physically clustered with Clade A, B or C ACS or PACS genes ( Figure 3C,D). We instead identified a different scaffold in I. luminosus that is likely orthologous to the firefly Luc1 locus (Scaffold 9654; Figure 3D). This assessment is based on the presence of adjacent Clade A and B ACS and PACS genes, as well as orthologous exoribonuclease family (PRNT) and inositol monophosphatase family (IMP) genes, both of which were found adjacent to the A. lateralis Luc1 locus, but not the P. pyralis Luc1 locus ( Figure 3D). Interestingly, IlumPACS11, the most early-diverging member of Clade D, was also found on Scaffold 9654 ( Figure 3D). This finding is consistent with an expansion of Clade D following duplication of the IlumPACS11 syntenic ancestor to a distant site. Overall, these genomic structures are consistent with independent origins of firefly and click beetle luciferases.
We then carried out targeted molecular evolution analyses including the known beetle luciferases and their closely related homologs. Ancestral state reconstruction of luminescent activity on the gene tree using Mesquite (Maddison and Maddison, 2017) recovered two independent gains of luminescence as the most parsimonious and likely scenario: once in click beetles, and once in the common ancestor of firefly, phengodid, and rhagophthalmid beetles ( Figure 4A; Appendix 4.3.3). In an independent molecular adaptation analysis utilizing the coding nucleotide sequence of the elaterid luciferases and their close homologs within Elateridae, 35% of the sites of the branch leading to the ancestral click beetle luciferase showed a statistically significant signal of episodic positive selection with d N /d S > 1 (w or max d N /d S = 3.98) as compared to the evolution of its paralogs using the aBSREL branch-site selection test (Smith et al., 2015) ( Figure 4B; Appendix 4.3.4). This implies that the common ancestor of the click beetle luciferases (EAncLuc) underwent a period of accelerated directional evolution. As the branch under selection in the molecular adaptation analysis ( Figure 4B) is the same branch of luciferase activity gain via ancestral reconstruction ( Figure 4A), we conclude that the identified selection signal represents the relatively recent neofunctionalization of click beetle luciferase from a non-luminous ancestral Clade D PACS gene, distinct from the more ancient neofunctionalization of firefly luciferase. Based on the constraints from our tree, we determine that this neofunctionalization of EAncLuc occured after the divergence of the elaterid subfamily Agrypninae. In contrast, we cannot determine if the original neofunctionalization of AncLuc occurred in the ancestral firefly, or at some point during the evolution of 'cantharoid' beetles, an unofficial group of beetles including the luminous Rhagophthalmidae, Phengodidae and Lampyridae among other non-luminous groups, but not the Elateridae (Branham and Wenzel, 2003). There is evidence for a subsequent luciferase duplication event in phengodids, but not in rhagophthalmids, that is independent of the duplication event that gave rise to Luc1 and Luc2 in fireflies ( Figures 3C and  4). Altogether, our results strongly support the independent neofunctionalization of luciferase activity in click beetles and fireflies, and therefore at least two independent gains of luciferin-utilizing luminescence in beetles.

Metabolic adaptation of the firefly lantern
Beyond luciferase, we sought to characterize other metabolic traits which might have co-evolved in fireflies to support bioluminescence. Of particular importance, the enzymes of the de novo biosynthetic pathway for firefly luciferin remain unknown (Oba et al., 2013b). We hypothesized that bioluminescent accessory enzymes, either specialized enzymes with unique functions in luciferin metabolism or enzymes with primary metabolic functions relevant to bioluminescence, would be highly expressed (HE: 90th percentile; Appendix 4.2.2) in the adult lantern, and would be differentially expressed (DE; Appendix 4.2.2) between luminescent and non-luminescent tissues. To determine this, we performed RNA-Seq and expression analysis of the dissected P. pyralis and A. lateralis adult male lantern tissue compared with a non-luminescent tissue (Appendix 4.2.2). We identified a set of predicted orthologous enzyme-encoding genes conserved in both P. pyralis and A. lateralis that met our HE and DE criteria ( Figure 5). Both luciferase and luciferin sulfotransferase (LST), a specialized enzyme recently implicated in luciferin storage in P. pyralis , were recovered as candidate genes using four criteria (HE, DE, enzymes, direct orthology across species), confirming the validity of our approach. While a direct ortholog of LST is present in A. lateralis, it is absent from I. luminosus, suggesting that LST, and the presumed luciferin storage it mediates, is an exclusive ancestral firefly or cantharoid trait. This finding is consistent with previous hypotheses of the absence of LST in Elateridae , and with the overall hypothesis of independent evolution of bioluminescence between the Lampyridae and Elateridae.
Moreover, we identified several additional enzyme-encoding HE and DE lantern genes that are likely important in firefly lantern physiology ( Figure 5). For instance, adenylate kinase likely plays a critical role in efficient recycling of AMP post-luminescence, and cystathionine gamma-lyase supports a key role of cysteine in luciferin biosynthesis (Oba et al., 2013b) and recycling (Okada et al., 1974). We also detected a combined adenylyl-sulfate kinase and sulfate adenylyltransferase enzyme (ASKSA) among the lantern-enriched gene list (Appendix 4-figure 8), implicating active biosynthesis of 3'-phosphoadenosine-5'-phosphosulfate (PAPS), the cofactor of LST, in the lantern. This finding highlights the importance of LST-catalyzed luciferin sulfonation for bioluminescence. These firefly orthologs of ASKSA are the only members amongst their paralogs to contain a PTS1 (Appendix 4figure 8), suggesting specialized localization to the peroxisome, the location of the luminescence reaction. This suggests that the levels of sulfoluciferin and luciferin may be actively regulated within the peroxisome of lantern cells in response to luminescence. Overall, our findings of several directly orthologous enzymes that share expression patterns in the light organs of both P. pyralis and A. lateralis suggests that the enzymatic physiology and/or the gene expression patterns of the photocytes were already fixed in the Luciolinae-Lampyrinae ancestor.
We also performed a similar expression analysis for genes not annotated as enzymes, yielding several genes with predicted lysosomal function (Appendix 4-table 6; Appendix 4.4). This suggests that the abundant but as yet unidentified 'differentiated zone granule' organelles of the firefly light organ (Ghiradella and Schmidt, 2004) could be lysosomes. Interestingly, we found a HE (TPM value~300) and DE opsin, Rh7, in the light organ of A. lateralis, but not P. pyralis (Appendix 4-figure 9; Appendix 4.5), suggesting a potential light perception role for Rh7 in the A. lateralis lantern, akin to the light perception role described for Drosophila Rh7 (Ni et al., 2017).

Genomic insights into firefly chemical defense
Firefly bioluminescence is postulated to have first evolved as an aposematic warning of larval chemical defenses (Branham and Wenzel, 2003). Lucibufagins are abundant unpalatable defense steroids described from certain North American firefly species, most notably in the genera Photinus (Meinwald et al., 1979), Lucidota (Gronquist et al., 2005), and Ellychnia (Smedley et al., 2017), and hence are candidates for ancestral firefly defense compounds. To test whether lucibufagins are widespread among bioluminescent beetles, we assessed the presence of lucibufagins in P. pyralis, A. lateralis, and I. luminosus by liquid-chromatography high-resolution accurate-mass mass-spectrometry (LC-HRAM-MS). While lucibufagins were found in high abundance in P. pyralis adult hemolymph, they were not observed in A. lateralis adult hemolymph, nor in I. luminosus metathorax extract ( Figure 6B; Appendix 4.6). Since chemical defense is presumably most critical in the long- The following source data is available for figure 6: Source data 1. CYP303 multiple sequence alignment and gene tree for Figure 6C. DOI: https://doi.org /10.7554/eLife.36495.017 lived larval stage, we next tested whether lucibufagins are present in all firefly larvae even if they are not present in the adults of certain species. We found lucibufagins in P. pyralis larval extracts; however, they were not observed in A. lateralis larval extracts ( Figure 6B; Appendix 4.6). Together, these results suggest that the lucibufagin biosynthetic pathway is either a derived trait only found in particular firefly taxa (e.g. subfamily: Lampyrinae), or that lucibufagin biosynthesis was an ancestral trait that was lost in A. lateralis. Consistent with the former hypothesis, the presence of lucibufagins in non-North-American Lampyrinae has been previously reported (Tyler et al., 2008), but to date there are no reports of lucibufagins in the Luciolinae. The lucibufagin biosynthetic pathway is currently unknown. However, their chemical structure suggests a biosynthetic origin from cholesterol followed by a series of hydroxylations, -OH acetylations, and the side-chain oxidative pyrone formation ( Figure 6A) (Meinwald et al., 1979). We hypothesized that cytochrome P450s, an enzyme family widely involved in metabolic diversification of organic substrates (Hamberger and Bak, 2013), could underlie several oxidative reactions in the proposed lucibufagin biosynthetic pathway. We therefore inferred the P450 phylogeny among our three bioluminescent beetle genomes to identify any lineage-specific genes correlated with lucibufagin presence. Our analysis revealed a unique expansion of one P450 family, the CYP303 family, in P. pyralis. While 94/97 of currently sequenced winged-insect genomes on OrthoDB (Zdobnov et al., 2017), as well as the A. lateralis and I. luminosus genomes, contain only a single CYP303 family gene, the P. pyralis genome contains 11 CYP303 genes and two pseudogenes ( Figure 6C), which expanded via tandem duplication on the same linkage group ( Figure 6D). The CYP303 ortholog of D. melanogaster, CYP303A1, has been shown to play a role in mechanosensory bristle development (Willingham and Keil, 2004). Although the exact biochemical function and substrate of D. melanogaster CYP303A1 is unknown, its closely related P450 families operate on an insect steroid hormone ecdysone (Willingham and Keil, 2004). As ecdysone and lucibufagins are structurally similar, CYP303 may operate on steroid-like compounds. Therefore, the lineage-specific expansion of the CYP303 family in P. pyralis is a compelling candidate in the metabolic evolution of lucibufagins as chemical defenses associated with the aposematic role of bioluminescence. Alternatively, this CYP303 expansion in P. pyralis may be associated with other lineage-specific chemical traits, such as pheromone production.

Symbionts of bioluminescent beetles
Given the increasingly recognized contributions of symbionts to host metabolism (Newman and Cragg, 2015), we characterized the hologenome of all three beetles as potential contributors to metabolic processes related to bioluminescence. Whole genome sequencing of our wild-caught and laboratory reared fireflies revealed a rich microbiome. Amongst our firefly genomes, we found various bacterial genomes, viral genomes, and the complete mtDNA for a phorid parasitoid fly, Apocephalus antennatus, the first mtDNA reported for genus Apocephalus. This mtDNA was inadvertently included in the P. pyralis PacBio library via undetected parasitization of the initial specimens, and was assembled via a metagenomic approach (Appendix 5.2). Independent collection of A. antennatus which emerged from field-collected P. pyralis adults and targeted COI sequencing later confirmed the taxonomic origin of this mtDNA (Appendix 5.3). We also sequenced and metagenomically assembled the complete circular genome (1.29 Mbp, GC: 29.7%;~50x coverage) for a P. pyralis-associated mollicute (Phylum: Tenericutes), Entomoplasma luminosum subsp. pyralis (Appendix 5.1). Entomoplasma spp. were first isolated from the guts of North American fireflies (Hackett et al., 1992) and our assembly provides the first complete genomic assembly of any Entomoplasma species. Broad read coverage for the E. luminosus subsp. pyralis genome was detected in 5/6 of our P. pyralis DNA libraries, suggesting that Entomplasma is a highly prevalent, possibly vertically inherited, P. pyralis symbiont. It has been hypothesized that these Entomoplasma mollicutes could play a role in firefly metabolism, specifically via contributing to cholesterol metabolism and lucibufagin biosynthesis (Smedley et al., 2017).
Lastly, we detected two species of novel orthomyxoviridae-like ssRNA viruses, which we dub Photinus pyralis orthomyxo-like virus 1 and 2 (PpyrOMLV1/2), that were highly prevalent across our P. pyralis RNA-Seq datasets, and showed multi-generational transovarial transmission in the laboratory (Appendix 5.4). We also found several endogenous viral elements (EVEs) for PpyrOMLV1/2 in P. pyralis (Appendix 5.5). These viruses are the first reported in any firefly species, and represent only the second report of transgenerational transfer of any Orthomyxoviridae virus (Marshall et al., 2014), and the second report of Orthomyxoviridae derived EVEs (Katzourakis and Gifford, 2010). Together, these genomes from the firefly holobiont provide valuable resources for the continued inquiry of the symbiotic associates of fireflies and their biological and ecological significance.

Discussion
Here, we generated genome assembles, diverse tissue and life-stage RNA-Seq data, and LC/MS data for three evolutionarily informative and historically well-studied bioluminescent beetles, and used a series of comparative analyses to illuminate long-standing questions on the origins and evolution of beetle bioluminescence. By analyzing the genomic synteny and molecular evolution of the beetle luciferases and their extant and inferred-ancestral homologs, we found strong support for the independent origins of luciferase, and therefore bioluminescence, between fireflies and click beetles. Our approaches and analyses lend molecular evidence to the previous morphology-phylogeny based hypotheses of parallel gain proposed by Darwin and others (Darwin, 1872;Branham and Wenzel, 2003;Costa, 1975;Sagegami-Oba et al., 2007;Bocakova et al., 2007;Oba, 2009;Day, 2013). While our elaterid luciferase selection analysis strongly supports an independent gain, we did not perform an analogous selection analysis of luciferase homologs across all bioluminescent beetles, due to the lack of genomic data from key related beetle families. Additional genomic information from early-diverged firefly lineages, other luminous beetle taxa (e.g. Phengodidae and Rhagophthalmidae), and non-luminous elateroid taxa (e.g. Cantharidae and Lycidae), will be useful to further develop and test models of luciferase evolution, including the hypothesis that bioluminescence also originated independently in the Phengodidae and/or Rhagophthalmidae. As some phylogenetic relationships of fireflies and other lineages of superfamily Elateroidea remain uncertain, continued efforts to produce reference phylogeny for these taxa are required (Martin et al., 2017;Bocak et al., 2018). Toward this goal, the recently published Pyrocoelia pectoralis Lampyrinae firefly genome is an important advance which will contribute to future phylogenetic and evolutionary studies (Fu et al., 2017). The independent origins of the firefly and click beetle luciferases provide an exemplary natural model system to understand enzyme evolution through parallel mutational trajectories and the evolution of complex metabolic traits generally. The abundance of gene duplication events of PACSs and ACSs at the ancestral luciferase locus in both fireflies and I. luminosus suggests that ancestral promiscuous enzymatic activities served as raw materials for the selection of new adaptive catalytic functions (Weng, 2014). But while parallel evolution of luciferase implies evolutionary independence of bioluminescence overall, the reality may be more complex, and the other subtraits of bioluminescence amongst the bioluminescent beetles likely possess different evolutionary histories from luciferase. While subtraits presumably dependent on an efficient luciferase, such as specialized tissues and neural control, almost certainly arose well after luciferase specialization, and thus can be inferred to also have independent origins between fireflies and click beetles, luciferin, which was presumably a prerequisite to luciferase neofunctionalization, may have been present in their common ancestor. Microbial endosymbionts, such as the tenericutes detected in our P. pyralis and A. lateralis datasets, are intriguing candidate contributors to luciferin metabolism and biosynthesis. Alternatively, recent reports have shown that firefly luciferin is readily produced non-enzymatically by mixing benzoquinone and cysteine (Kanie et al., 2016), and that a compound resulting from the spontaneous coupling of benzoquinone and cysteine acts as a luciferin biosynthetic intermediate in A. lateralis (Kanie et al., 2018). Benzoquinone is known to be a defense compound of distantly related beetles (Dettner, 1987) and other arthropods (e.g. millipedes) (Shear, 2015). Therefore, the evolutionary role of sporadic low-level luciferin synthesis through spontaneous chemical reactions, either in the ancestral bioluminescent taxa themselves, or in non-bioluminescent taxa, and dietary acquisition of luciferin by either the ancestral or modern bioluminescent taxa, should be considered. To decipher between these alternative evolutionary possibilities, the discovery of genes involved in luciferin metabolism in fireflies and other bioluminescent beetles will be essential. Here, as a first step toward that goal, we identified conserved, enriched and highly expressed enzymes of the firefly lantern that are strong candidates in luciferin metabolism and the elusive luciferin de novo biosynthetic pathway. Ultimately focused experimentation will be needed to decipher the biochemical function of these enzymes.
The early evolution of firefly bioluminescence was likely associated with an aposematic role. The adaptive light production of the primordial firefly (or alternatively, a primordial bioluminescent cantharoid beetle) that enabled the selection and neofunctionalization of luciferase was perhaps linked to a response to predators by a primitive whole-body oxygen-gated luminescence, where a startleresponse mediated increase in hemolymph oxygenation through spiracle opening and escape locomotion caused a concomitant increase in luminescence (Buck and Case, 2002;Case, 2004). Alternatively, an early role for firefly luminescence in mate attraction has not been ruled out (Buck and Case, 2002). The presence of particular unpalatable defense compounds in all extant fireflies would be consistent with an ancestral role and the former hypothesis, and the chemical analysis of tissues across species and life stages presented in this work provides new insights into the evolutionary occurrence of lucibufagins, the most well-studied defense compounds associated with fireflies. Our results reject lucibufagins as ancestral defense compounds of fireflies, but rather suggest them as a derived metabolic trait associated with Lampyrinae. Additional chemical analyses across more lineages of fireflies are needed, however, to further support or falsify this hypothesis. Toward this goal, the high sensitivity of our LC-HRAM-MS and MS 2 molecular networking-based lucibufagin identification approach is particularly well suited to broadened sampling in the future, including those of rare taxa and possibly museum specimens. Combined with genomic data showing a concomitant expansion of the CYP303 gene family in P. pyralis, we present a promising path toward elucidating the biosynthetic mechanism underlying these potent firefly toxins.
Overall, the resources and analyses generated in this study shed valuable light on the evolutionary questions Darwin first pondered, and will enable future studies of the ecology, behavior, and evolution of bioluminescent beetles. These resources will also accelerate the discovery of new enzymes from bioluminescent beetles that could enhance biotechnological applications of bioluminescence. Finally, we hope that the genomic resources shared here will facilitate the development of effective population genomic tools to monitor and protect wild bioluminescent beetle populations in the face of changing climate and habitats.

Materials and methods
Detailed materials and methods are available in the Appendices. Methods relating to P. pyralis are given in Appendix 1, while methods relating to A. lateralis and I. luminosus are given in Appendix 2 and Appendix 3, respectively. Methods for comparative genomic analyses are given in Appendix 4, while methods for microbiome characterization are given in Appendix 5. References to relevant sections of the Appendices are placed in-line throughout the maintext.

Data and materials availability
Genomic assemblies (Ppyr1.3,Alat1.3,and Ilumi1.2), associated official geneset data, a Sequence-Server (Priyam et al., 2015) BLAST server, and a JBrowse (Skinner et al., 2009) (Lloyd, 1966;Lloyd, 2008). It inspired extensive work on the biochemistry and physiology of firefly bioluminescence in the early 20th century, and the first luciferase gene was cloned from this species (de Wet et al., 1985). A habitat generalist, P. pyralis occurs in fields, meadows, suburban lawns, forests, and woodland edges, and even urban environments. For example, the authors have observed P. pyralis flashing in urban New York City and Washington D.C. Adults rest on vegetation during the day and signaling begins as early as 20 min before sunset (Lloyd, 1966). Male flashing is cued by ambient light levels, thus shaded or unshaded habitats can show up to a 30 min difference in the initiation of male flashing (Lloyd, 1966). Males can be cued to flash outside of true twilight if exposed to light intensities simulating twilight (Case, 2004). P. pyralis were also reported to flash during totality of the total solar eclipse of 2017 (Personal communication: L.F. Faust, M.A. Branham). Courtship activity lasts for 30-45 min and both sexes participate in a bioluminescent flash dialog, as is typical for Photinus fireflies. Males initiate courtship by flying low above the ground while repeating a single~300 ms patrol flash at~5-10 s intervals (Case, 2004). Males emit their patrol flash while dipping down and then ascending vertically, creating a distinctive J-shaped flash gesture (Lloyd, 1966;Case, 2004) ( Figure 1A). During courtship, females perch on vegetation and respond to a male patrol flash by twisting their abdomen toward the source of the flash and giving a single response flash given after a 2-3 s delay (Video 1). Receptive females will readily respond to simulated male flashes, such as those produced by an investigator's penlight. Females have fully developed wings and are capable of flight. Both sexes are capable of mating several times during their adult lives. During mating, males transfer to females a fitness-enhancing nuptial gift consisting of a spermatophore manufactured by multiple accessory glands (van der Reijden et al., 1997); the molecular composition of this nuptial gift has recently been elucidated for P. pyralis (Al-Wathiqui et al., 2016). In other Photinus species, male gift size decreases across sequential matings (Cratsley et al., 2003), and multiple matings are associated with increased female fecundity (Rooney and Lewis, 2002).
Adult P. pyralis live 2-3 weeks, and although these adults are typically considered nonfeeding, both sexes have been reported drinking nectar from the flowers of the milkweed Asclepias syriaca (Faust and Faust, 2014). Mated females store sperm and lay~30-50 eggs over the course of a few days on moss or in moist soil. The eggs take 2-3 weeks to hatch. Larval bioluminescence is thought to be universal for the Lampyridae, where it appears to function as an aposematic warning signal. Like other Photinus, P. pyralis larvae are predatory, live on and beneath the soil, and appear to be earthworm specialists (Hess, 1920). In the northern parts of its range, slower development likely requires P. pyralis to overwinter at least twice, most likely as larvae. Farther south, P. pyralis may complete development within several months, achieving two generations per year (Faust, 2017), which may be possibly be observed in the South as a 'second wave' of signalling P. pyralis in September-October.
Anti-predator chemical defenses of male P. pyralis include several bufadienolides, known as lucibufagins, that circulate in the hemolymph (Meinwald et al., 1979). Pterins have also been reported to be abundant in P. pyralis (Goetz et al., 1981); however, the potential defense role of these compounds has never been tested (Personal communication: J. Meinwald). When attacked, P. pyralis males release copious amounts of rapidly coagulating hemolymph and such 'reflex-bleeding' may also provide physical protection against small predators (Blum and Sannasi, 1974;Faust et al., 2012).

Species distribution
Although Photinus pyralis is widely distributed in the Eastern United States, published descriptions of its range are limited, with the notable exception of Lloyd's 1966 monograph (Lloyd, 1966) which addresses the range of many Photinus species. We therefore sought to characterize the current distribution of P. pyralis in order to produce an updated map to inform our experimental design and enable future population genetic studies. Four sources of data were used to produce the presented range map of P. pyralis: (i) Field surveys by the authors (ii) Published (Lloyd, 1966;Luk et al., 2011) and unpublished sightings of P. pyralis at county level resolution, provided by Dr. J. Lloyd (University of Florida), (iii) coordinates and dates of P. pyralis sightings, obtained by targeted e-mail surveys to firefly field biologists, (iv) citizen scientist reports of P. pyralis through the iNaturalist platform (iNaturalist, 2017). iNaturalist sightings were manually curated to only include reports which could be unambiguously identified as P. pyralis from the photos, and also that also included GPS geotagging to <100 m accuracy. A spreadsheet of these sightings is available on FigShare (DOI: 10.6084/m9.figshare.5688826). QGIS (v2.18.9, OSG Foundation, 2017) was used for data viewing and figure creation. A custom Python script (Fallon, 2018e; copy archived at https://github.com/elifesciencespublications/2017_misc_scripts) within QGIS was used to link P. pyralis sightings to counties from the US census shapefile (United States Census Bureau, 2017). Outlying points that were located in Desert Ecoregions of the World Wildlife Fund (WWF) Terrestrial Ecoregions shapefile (Olson et al., 2001;World Wildlife Fund, 2017) or the westernmost edge of the range were manually removed, as they are likely isolated populations not representative of the contiguous range. For Figure 1B, these points were converted to a polygonal range map using the 'Concave hull' QGIS plugin ('nearest neighbors = 19') followed by smoothing with the Generalizer QGIS plugin with Chaiken's algorithm (Level = 10, and Weight = 3.00). Below (Appendix 1-figure 1), red circles indicate county-centroided presence records.
Appendix 1-figure 1. Detailed geographic distribution map for P. pyralis. P. pyralis sightings (red circles show county centroided reports) in the United States and Ontario, Canada (diagonal hashes). The World Wildlife Fund Terrestrial Ecoregions (Olson et al., 2001;World Wildlife Fund, 2017) are also shown (colored shapes). The P. pyralis sighting dataset shown is identical to that used to prepare Figure 1B.
In our field surveys, we found that the range of P. pyralis was notably extended from the range reported by Lloyd, specifically we found P. pyralis in abundance to the west of the Mill river in Connecticut. P. pyralis is found with confidence roughly from Connecticut to Texas, and possibly as far south as Guatemala (Personal communication: A. Catalá n). These possible southern populations require further study.

Specimen collection and identification
Adult male P. pyralis specimens for Illumina short-insert and mate-pair sequencing were collected at sunset on June 13th, 2011 near the Visitor's Center at Great Smoky Mountains National Park (permit to Dr. Kathrin Stanger-Hall). Specimens were identified to species and sex via morphology (Green, 1956), flash pattern and behavior (Lloyd, 1966), and cytochromeoxidase I (COI) similarity (partial sequence: primers HCO, LCO [Stanger-Hall and Lloyd, 2015]) when blasted against an in-house database of firefly COI nucleotide sequences. Collected fireflies were stored in 95% ethanol at À80˚C until DNA extraction.
Adult male P. pyralis specimens for Pacific Biosciences (PacBio) RSII sequencing were captured during flight at sunset on June 9th, 2016, from Mercer Meadows in Lawrenceville, NJ (40.3065 N 74.74831 W), on the basis of the characteristic 'rising J' flash pattern of P. pyralis (permit to TRF via Mercer County Parks Commission). Collected fireflies were sorted, briefly checked to be likely P. pyralis by the presence of the margin of ventral unpigmented abdominal tissue anterior to the lanterns, flash frozen with liquid N 2 , lyophilized, and stored at À80˚C until DNA extraction. A single aedeagus (male genitalia) was dissected from the stored specimens and confirmed to match the P. pyralis taxonomic key (Green, 1956) (Appendix 1figure 2).

A B
Appendix 1-figure 2. P. pyralis aedeagus (male genitalia). (A) Ventral and (B) side view of a P. pyralis aedeagus dissected from specimens collected on the same date and locality as those used for PacBio sequencing. Note the strongly sclerotized paired ventro-basal processes ('mickey mouse ears') emerging from the median process, characteristic of P. pyralis (Green, 1956).

Collection and rearing of P. pyralis larvae
We intended to survey the lucibufagin content of P. pyralis larvae ( Figure 6B; Appendix 4.6), and as well as the transovarial transmission of Photinus pyralis orthomyxo-like viruses from parent to larvae (Appendix 5.4), but as P. pyralis larvae are subterranean and extremely difficult to collect from the wild, we reared P. pyralis larvae from eggs laid from mated pairs. It is important to note that these P. pyralis larval rearing experiments were unexpectedly successful. Although there has been some success in laboratory rearing and domestication of Asian Aquatica spp. (Chiang and Yang, 2010), including the A. lateralis Ikeya-Y90 strain described in this manuscript, rearing of North American fireflies is considered extremely difficult with numerous unpublished failures for unclear reasons (Lloyd, 1996), and limited reports of successful rearing of mostly non-Photinus genera, including Photuris sp. (McLean et al., 1972), Pyractomena angulata (Buschman, 1988), and Pyractomena borealis (Personal communication: Scott Smedley). The below protocol for P. pyralis larval rearing is presented in the context of disclosure of the methods of this manuscript, and should be considered a preliminary, unoptimized rearing protocol. A full description of the P. pyralis larvae and it's life history and behavior will be presented in a separate manuscript. Four adult female P. pyralis were collected from the Bluemont Junction Trail in Arlington, VA from June 12th through June 18th 2017 (collection permission obtained by TRF from Arlington County Parks and Recreation department). The females were mated to P. pyralis males collected either from the same locality and date, or to males collected from Kansas in late June. Mating was performed by housing one to two males and one female in small plastic containers for~1-3 days with a wet kimwipe to maintain humidity. Mating pairs were periodically checked for active mating, which in Photinus fireflies takes several hours. Successfully mated females were transferred to Magenta GA-7 plastic boxes (Sigma-Aldrich, USA), and provided a~4 cm x 4 cm piece of locally collected moss (species diverse and unknown) as egg deposition substrate, and allowed to deposit eggs until their death in~1-4 days. Deceased females were removed, artificial freshwater (AFW; 1:1000 diluted 32 PSU artificial seawater) was sprayed into the box to maintain high humidity, and eggs were kept for 2-3 weeks at room temperature and periodically checked until hatching. Like other firefly eggs, the eggs of P. pyralis were observed to be faintly luminescent imaging using a cooled CCD camera (Appendix 1- figure 3); however, this luminescence was not visible to the darkadapted eye, indicating that this luminescence is less intense than other firefly species such as Luciola cruciata (Harvey, 1952).
Upon hatching, first instar larvae were mainly fed~1 cm cut pieces of Canadian Nightcrawler earthworms (Lumbricus terrestris; Windsor Wholesale Bait, Ontario, Canada), and occasional live White Worms (Enchytraeus albidus; Angels Plus, Olean, NY). Although P. pyralis first instar larvae were observed to attack live Enchytraeus albidus, an experiment to determine if this would be suitable as a single food source was not performed. Uneaten and putrefying earthworm pieces were removed after 1 day, and the container cleaned. Once the larvae had been manually fed for~2 weeks and deemed sufficiently strong, they were transferred to plastic shoeboxes (P/N: S-15402, ULINE, USA) which were intended to mimic a soil ecosystem. In personal discussions of unpublished firefly rearing attempts by various firefly researchers, we noted that a common theme was the difficulty of preventing the uneaten prey of these predatory larvae from putrifying. Therefore, we sought to create ecologically inspired 'eco-shoeboxes', where fireflies would prey on live organisms, and other organisms would assist in cleanup of uneaten or partially eaten prey that had been fed to the firefly larvae, to prevent the growth of pathogenic microorganisms on uneaten prey.
First, these shoeboxes were filled with 1L of mixed 50% (v/v) potting soil, and 50% coarse sand (Quikrete, USA) that had been washed several times with distilled water to remove silt and dust. The soil-sand mix was wet well with AFW, and live Enchytraeus albidus (50+), temperate springtails (50+; Folsomia candida; Ready Reptile Feeders, USA), and dwarf isopods (50+; Trichorhina tomentosa; Ready Reptile Feeders, USA) were added to the box, and several types of moss, coconut husk, and decaying leaves were sparingly added to the corners of the box. The non-firefly organisms were included to mimic a primitive detritivore (Enchytraeus albidus and Trichorhina tomentosa) and fungivore (Folsomia candida) system. About 50 firefly larvae were included per box. No interactions between the P. pyralis larvae and the additional organisms were observed. Predation on Enchytraeus albidus seems likely, but careful observations were not made. Distilled water was sprayed into the box every~2 days to maintain a high humidity. Throughout this period, live Lumbricus terrestris (~10-15 cm) were added to the box every 2-3 days as food. These earthworms were first prepared by washing with distilled water several times to remove attached soil, weakened and stimulated to secrete coelemic fluid and gut contents by spraying with 95% ethanol, washed several times in distilled water, and left overnight in~2 cm depth distilled water at 4˚C. Anecdotally this precleaning and preparation process reduced the rate and degree that dead earthworms putrefied. Young P. pyralis larvae were observed to successfully kill and gregariously feed on these live earthworms (Appendix 1- figure 4). The possibility that firefly larvae possess a paralytic venom used to stun or kill prey has been noted by other researchers (Hess, 1920;Williams, 1917). In our observations, an earthworm would immediately react to the bite from a single P. pyralis larvae, thrashing about for several minutes, but would then become seemingly paralyzed over time, supporting the role of a potent, possibly neurotoxic, firefly venom. The P. pyralis larvae would then begin extra-oral digestion and gregarious feeding on the liquified earthworm. Once the earthworm had been killed and broken apart by firefly larvae, Enchytraeus albidus would enter through gaps in the cuticle and begin to feed in large numbers throughout the interior of the earthworm. The other detritivores were observed at later stages of feeding. Between the combined action of the P. pyralis larvae, and the other detritivores, the live earthworm was completely consumed within 1-2 days, and no manual cleanup was required.
Compared to the initial manual feeding and cleaning protocol for P. pyralis 1st instar larvae, the 'eco-shoebox' rearing method was low-input and convenient for large numbers of larvae. The feeding and cleanup process was efficient for~2 months (July through September), leading to a large number of healthy 3-4th instar larvae (Appendix 1- figure 5). However, after that point, P. pyralis larvae, possibly in preparation for a winter hibernation, seemingly became quiescent, and were less frequently seen patrolling throughout the box. At the same time, the Enchytraeus albidus earthworms were observed to become less abundant, either due to continual predation by P. pyralis, or due to population collapse from insufficient fulfillment of nutritional requirements from feeding of Enchytraeus albidus on Lumbricus terrestris alone.
At this point, earthworms were not consumed within 1-2 days, and became putrid, and P. pyralis which had been feeding on these earthworms were frequently found dead nearby, and themselves quickly putrefied. Generally after this point P. pyralis larvae were more frequently found dead and partially decayed, indicating the possibility of pathogenesis from microorganisms from putrefying earthworms. At this stage, it was observed that mites (Acari), probably from the soil contained in the guts of the fed earthworms, became abundant, and were observed to act as ectoparasitic on P. pyralis larvae. An attempt to simulate hibernation of P. pyralis larvae was made by storing them at 4˚C for~3 weeks, however a large proportion (~30%) of larvae died during this hibernation to a seeming fungal infection. Other larvae revived quickly when returned to room temperature, but all Trichorhina tomentosa were killed by even transient exposure to 4˚C. To date, a smaller number of fifth and sixth instar P. larvae have been obtained, but pupation in the laboratory has not occured. The lack of pupation is unsurprising as it is likely occurs in the wild after 1-2 years of growth, is likely under temperature and photoperiodic control, and may require a licensing stage of cold temperature hibernation for several weeks. Overall, manual feeding of first1 st instar larvae followed by the 'eco-shoebox' method was unexpectedly successful approach for the maintenance and growth of P. pyralis larvae. Appendix 1-figure 4. Gregarious predation of young P. pyralis larvae on a live Lumbricus terrestris. Both P. pyralis larvae (red arrows), and Enchytraeus albidus (yellow arrows), were observed to feed on the paralyzed earthworms.
The karyotype of P. pyralis was previously reported to be 2n = 20 with XO sex determination (male, 18A + XO; female, 18A + XX) (Wasserman and Ehrman, 1986). The genome sizes of four P. pyralis adult males were previously determined to be 422 ± 9 Mbp (SEM, n = 4), whereas the genome sizes of five P. pyralis adult females were determined to be 448 ± 7 (SEM, n = 5) by nuclear flow cytometry analysis (Lower et al., 2017). From these analyses, the size of the X-chromosome is inferred to be~26 Mbp. Genome size inference via kmer spectral analysis of the P. pyralis short-insert Illumina data from a single adult P. pyralis male estimated a genome size of 343 Mbp (Appendix 1-figure 6).

Library preparation and sequencing
See Appendix 4-table 1 for a overview of all sequence libraries. Library specific construction methods are detailed below.

Illumina
DNA was extracted from sterile-water-washed thorax of Great Smoky Mountains National Park collected specimens using phenol-chloroform extraction with RNAse digestion, checked for quality via gel electrophoresis, and quantified by Nanodrop or Qubit (Thermo Scientific, USA). To obtain sufficient DNA for both short insert and mate-pair library construction, libraries were constructed separately from DNA from each of two individual males and pooled DNA of three males, all from the same population. Males were selected for sequencing as they are more easily found in the field than females. In addition, as P. pyralis males are XO (Dias et al., 2007), differences in sequencing coverage could inform localization of scaffolds to the X chromosome. Illumina TruSeq short insert (average insert size: 300 bp) and Nextera mate-pair libraries (insert size: 3 Kbp, 6 Kbp) were constructed at the Georgia Genomics Facility (Athens, GA) and subsequently sequenced on two lanes of Illumina HiSeq2000 100 Â 100 bp PE reads (University of Texas; Appendix 4-table 1).

A B
Appendix 1 . len = inferred haploid genome length, uniq = percentage non-repetitive sequence, het = overall rate of genome heterozygosity, kcov = mean kmer coverage for heterozygous bases, err = error rate of the reads, dup: average rate of read duplications. These results are consistent with the genome size of a XO male, when possible systematic error of kmer spectral analysis and flow cytometry genome size estimates is considered. The heterozygosity is somewhat low when compared to some other arthropods.

PacBio
High-molecular-weight DNA (HMW DNA) was extracted from four pooled lyophilized adult male P. pyralis (dry mass 90.8 mg) from the MMNJ field site. These specimens were first externally washed using 95% ethanol, after which DNA extraction proceeded with a 100/G Genomic Tip plus Genomic Buffers kit (Qiagen, USA). DNA extraction followed the manufacturer's protocol, with the exception of the final precipitation step, where HMW DNA was pelleted with 40 mg RNA grade glycogen (Thermo Scientific, USA) and centrifugation (3000 x g, 30 min, 4˚C) instead of spooling on a glass rod. Although increased genomic heterozygosity from four pooled males and a resulting more complicated genome assembly was a concern for a wild population like P. pyralis, four males were used in order to extract enough DNA for workable coverage using 15 Kbp+ size selected PacBio RSII sequencing. All extracted DNA was used for library preparation, and all of the final library was used for sequencing. Adult males, being XO, were chosen over the preferable XX females, as adult males are much more easily captured because they signal during flight, whereas females are typically found in the brush below and generally only flash in response to authentic male signals. Precipitated HMW DNA was redissolved in 80 mL Qiagen QLE buffer (10 mM Tris-Cl, 0.1 mM EDTA, pH 8.5) yielding 17.1 mg of DNA (214 ng/mL) and glycogen (500 ng/mL). Final DNA concentration was measured with a Qubit fluorometer (Thermo Scientific) using the Qubit Broad Range kit. Manipulations hereafter, including HMW DNA size QC, fragmentation, size selection, library construction, and PacBio RSII sequencing, were performed by the Broad Technology Labs of the Broad Institute (Cambridge, MA).
First, the size distribution of the HMW DNA was confirmed by pulsed-field-gelelectrophoresis (PFGE). In brief, 100 ng of HMW DNA was run on a 1% agarose gel (in 0.5x TBE) with the BioRad CHEF DRIII system. The sample was run out for 16 hr at six volts/cm with an angle of 120 degrees with a running temperature of 14˚C. The gel was stained with SYBRgreen dye (Thermo Scientific -Part No. S75683). 1 mg of 5 Kbp ladder (BioRad, part no 170-3624) was used as a standard. These results demonstrated the HMW DNA had a mean size of >48 Kbp (Appendix 1- figure 7). This pool of HMW DNA is designated 1611_PpyrPB1 (NCBI BioSample SAMN08132578).
Next, HMW DNA (17.1 mg) was sheared to a targeted average size of 20-30 Kbp by centrifugation in a Covaris g-Tube (part no. 520079) at 2500 x g for 2 min. SMRTbell libraries for sequencing on the PacBio platform were constructed according to the manufacturer's recommended protocol for 20 Kbp inserts, which includes size selection of library constructs larger than 15 Kbp using the BluePippin system (Sage Science, Beverly, MA). Two separate cassettes were run. In each cassette, two lanes were used in which there was 1362 ng/lane (PAC20kb kit). Constructs 15 Kbp and above were eluted over a period of 4 hr. An additional damage repair step was carried out post size-selection. Insert size range for the final library was determined using the Fragment Analyzer System (Advanced Analytical, Ankeney, IA). The size-selected SMRTbell library was then sequenced over 61 SMRT cells on a PacBio RSII instrument of the Broad Technology Labs (Cambridge, MA), using the P6 v.2 polymerase and the v.4 DNA Sequencing Reagent (P6-C4 chemistry; part numbers 100-372-700, 100-612-400). PacBio sequencing data is available on the NCBI Sequence Read Archive (Bioproject PRJNA378805).
Appendix 1-figure 7. PFGE of P. pyralis HMW DNA used for PacBio sequencing. Lane 1 was used for further library prep and sequencing, Lanes 2-5 represent separate batches of P. pyralis HMW DNA that was not used for PacBio sequencing. Lane 1 was used as it had the highest DNA yield, and an equivalent DNA size distribution to the other samples.

Hi-C library preparation
Two adult P. pyralis MMNJ males were flash frozen in liquid nitrogen, stored at À80˚C, and shipped on dry-ice to Phase Genomics (Seattle, WA). Manipulations hereafter occurred at Phase Genomics, following previously published protocols (Lieberman-Aiden et al., 2009;Burton et al., 2013;Bickhart et al., 2017). Briefly, a streamlined version of the standard Hi-C protocol (Lieberman-Aiden et al., 2009) was used to perform a series of steps resulting in proximity-ligated DNA fragments, in which physically proximate sequence fragments are joined into linear chimeric molecules. First, in vivo chromatin was cross-linked with formaldehyde, fixing physically proximate loci to each other. Chromatin was then extracted from cellular material and digested with the Sau3AI restriction enzyme, which cuts at the GATC motif. The resulting fragments were proximity ligated with biotinylated nucleotides and pulled down with streptavidin beads. These chimeric sequences were then sequenced with 80 bp PE sequencing on the Illumina NextSeq platform, resulting in Hi-C read pairs.

Ppyr1.1: MaSuRCA hybrid assembly
Several genome assembly approaches were evaluated with the general goal of maximizing conserved gene content and contiguity. The highest quality P. pyralis assembly was generated by a hybrid assembly approach using a customized MaSuRCA (v3.2.1_01032017) (Zimin et al., 2013;Zimin et al., 2017) pipeline that combined both Illumina-corrected PacBio reads (Mega-reads) and synthetic long reads constructed from short-insert reads alone (Super-reads) using a custom small overlap length (59 bp).
We first applied MaSuRCA (v3.2.1_01032017) (Zimin et al., 2013;Zimin et al., 2017) to correct our long reads (38x coverage; Library ID 1611_PpyrPB1; Appendix 4-table 1) using our short-insert and mate-pair reads (Libraries: 8369, 375_3K, 8375_6K, 83_3K, 83_6K; Appendix 4-table 1). No pre-filtering of reads was performed, as Illumina adaptors are automatically removed within the MaSuRCA pipeline. We modified the pipeline to assemble the genome using both corrected long reads (Mega-reads) and synthetic long reads (Superreads) with a custom smaller overlap length (59 bp). All reads (short-insert, mate-pair and PacBio) were then used within the MaSuRCA pipeline to call a genomic consensus.

Ppyr1.2: Scaffolding with Hi-C
The Hi-C read pairs were applied in a manner similar to that originally described here (Burton et al., 2013) and later expanded upon (Bickhart et al., 2017). Briefly, Hi-C reads were mapped to Ppyr1.1 with BWA (v1.7.13) (Li and Durbin, 2009), requiring perfect, unique mapping locations for a read pair to be considered usable. The number of read pairs joining a given pair of contigs is referred to as the 'link frequency' between those contigs, and when normalized by the number of restriction sites in the pair of contigs, is referred to as the 'link density' between those contigs.
A three-stage scaffolding process was used to create the final scaffolds, with each stage based upon previously described analysis of link density (Burton et al., 2013;Bickhart et al., 2017). First, contigs were placed into chromosomal groups. Second, contigs within each chromosomal group were placed into a linear order. Third, the orientation of each contig is determined. Each scaffolding stage was performed many times in order to optimize the scaffolds relative to expected Hi-C linkage characteristics.
In keeping with previously described methods (Burton et al., 2013;Bickhart et al., 2017), the number of chromosomal scaffolds to create-10-was an a priori input to the scaffolding process derived from the previously published chromosome count of P. pyralis (Wasserman and Ehrman, 1986). However, to verify the correctness of this assumption, scaffolds were created for haploid chromosome numbers ranging from 5 to 15. A scaffold number of 10 was found to be optimal for containing the largest proportion of Hi-C linkages within scaffolds, which is an expected characteristic of actual Hi-C data.
Mapping of Tribolium castaneum X chromosome proteins (NCBI Tcas 5.2) to the Ppyr1.2 assembly using both tblastn (v2.6.0) (Camacho et al., 2009) and Exonerate(v2.2.0) (Slater and Birney, 2005) based 'protein2genome' alignment through the MAKER pipeline revealed a relative enrichment on LG3a only. Taken together, this data suggested that the half-coverage section of LG3 (LG3a) corresponded to the X-chromosome of P. pyralis, and that it was misassembled onto an autosome. Therefore, we manually split LG3 into LG3a and LG3b in the final assembly. 1.6.4.2 Taxonomic annotation filtering Given the recognized importance of filtering genome assemblies to avoid misinterpretation of the data (Koutsovoulos et al., 2016), we sought to systematically remove assembled nonfirefly contaminant sequence from Ppyr1.2. Using the blobtools toolset (v1.0.1) (Laetsch and Blaxter, 2017), we taxonomically annotated our scaffolds by performing a blastn (v2.6.0+) nucleotide sequence similarity search against the NCBI nt database, and a diamond (v0.9.10.111) (Buchfink et al., 2015) translated nucleotide sequence similarity search against the of Uniprot reference proteomes (July 2017). Using this similarity information, we taxonomically annotated the scaffolds with blobtools using parameters '-x bestsumorder -rank phylum'. A tab delimited text file containing the results of this blobtools annotation are available on FigShare (DOI: 10.6084/m9.figshare.5688982). We then generated the final genome assembly by retaining scaffolds that either contained annotated features (genes or non-simple/low-complexity repeats), had coverage >10.0 in both the Illumina (Appendix 1- figure 9) and PacBio libraries (Appendix 1-figure 10), and if the taxonomic phylum was annotated as 'Arthropod' or 'no-hit' by the blobtools pipeline (Appendix 1- figure 11). This approach removed 374 scaffolds (2.1 Mbp), representing 15% of the scaffold number and 0.4% of the nucleotides of Ppyr1.2. Notably, four tenericute scaffolds, likely corresponding to a partially assembled Entomoplasma sp. genome, distinct from the Entomoplasma luminosus var. pyralis assembled from the PacBio library (Appendix 5) were removed. Furthermore, we removed two contigs representing the mitochondrial genome of P. pyralis (complete mtDNA available via Genbank: KY778696). The final filtered assembly, Ppyr1.3, is available at www. fireflybase.org.   In addition to our finalized genome assembly (Ppyr1.3), we sought to better understand the symbiont composition that varied between our P. pyralis PacBio and Illumina libraries. Therefore, we produced a long-read only assembly of our PacBio data to assemble the sequence that might be unique to this library. To achieve this, we first filtered the HDF5 data from the 61 sequence SMRT cells to. FASTQ format subreads using SMRTPortal (v2.3.0.140893) (Pacific Biosciences, 2017) with the 'RS_Subreads.1' protocol with default parameters. These subreads were then input into Canu (Github commit 28ecea5/v1.6)  with parameters 'genomeSize = 450 m corOutCoverage = 200 ovlErrorRate = 0.15 obtErrorRate = 0.15 -pacbio-raw'. The unpolished contigs from this produced genome assembly are dubbed Ppyr0.1-PB.

Mitochondrial genome assembly and annotation
To achieve a full length mitochondrial genome (mtDNA) assembly of P. pyralis, sequences were assembled separately from the nuclear genome. Short insert Illumina reads from a single GSMNP individual (Sample 8369; Appendix 4-table 1) were mapped to the known mtDNA of the closest available relative, Pyrocoelia rufa (NC_003970.1 [Bae et al., 2004]) using bowtie2 v2.3.1 (parameters: -very-sensitive-local). All concordant read pairs were input to SPAdes (v3.8.0) (Nurk et al., 2013) (parameters: -plasmid -only-assembler -k35,55,77,90) for assembly. The resulting contigs were then combined with the P. rufa mitochondrial reference genome for a second round of read mapping and assembly. The longest resulting contig aligned well to the P. rufa mitochondrial genome; however, it was~1 Kbp shorter than expected, with the unresolved region appearing to be the tandem repetitive region (TRU) (Bae et al., 2004), previously described in the P. rufa mitochondrial genome. To resolve this, all PacBio reads were mapped to the draft mitochondrial genome, and a single high-quality PacBio circular-consensus-sequencing (CCS) read that spanned the unresolved region was selected using manual inspection and manually assembled with the contiguous sequence from the Illumina sequencing to produce a complete circular assembly. The full assembly was confirmed by re-mapping the Illumina short-read data using bowtie2 followed by consensus calling with Pilon v1.21 (Walker et al., 2014). Re-mapped PacBio long-read data also confirmed the structure of the mtDNA, and indicated variability in the repeat unit copy number of the TRU amongst the four sequenced P. pyralis individuals (Sample 1611_PpyrPB1; Appendix 4-table 1). The P. pyralis mtDNA was then 'restarted' using seqkit , such that the FASTA record break occurred in the AT-rich region, and annotated using the MITOS2 annotation server (Bernt et al., 2017). Low confidence and duplicate gene predictions were manually removed from the MITOS2 annotation. The final P. pyralis mtDNA with annotations is available on GenBank (KY778696). In order to capture expression from diverse life stages, stranded RNA-Seq libraries were prepared from whole bodies of four life stages/sexes (eggs, 1 st instar larvae, adult male, and adult female; Appendix 1-table 1). Eggs and larvae were derived from a laboratory mating of P. pyralis (Collected MMNJ, July 2016). Briefly, live adult P. pyralis were transported to the lab and allowed to mate in a plastic container over several days. The female, later sequenced, was observed mating with two independent males on two separate nights. The female was then transferred to a plastic container with moss, and allowed to oviposit over several days. Once no more oviposition was observed, the female was removed, flash frozen with liquid N 2 , and stored at À80˚C for RNA extraction. Resulting eggs were washed 3x with dilute bleach/ H 2 O and reared in aggregate in plastic containers on moist Whatman paper.~13 days after the start of egg oviposition, a subset of eggs were flash frozen for RNA extraction. The remaining eggs were allowed to hatch and larvae were flash frozen the day after emergence (first instar). Total RNA was extracted from a single stored adult male (non-paternal to eggs/ larvae), the adult female (maternal to eggs/larvae), seven pooled eggs, and four pooled larvae using the RNeasy Lipid Tissue Mini Kit (QIAGEN) with the optional on-column DNase treatment. Illumina sequencing libraries were prepared by the Whitehead Genome Technology Core (WI-GTC) using the TruSeq Stranded mRNA library prep kit (Illumina) and following the manufacturer's instructions with modification to select for larger insert sizes (~300-350 bp). These samples were multiplexed with unrelated plant RNA-Seq samples and sequenced 150 Â 150 nt on one rapid mode flowcell (two lanes) of a HiSeq2500 (WI-GTC), to a depth of~30M paired reads per library.
To examine gene expression in adult light organs, we generated non-strand specific sequencing of polyA pulldown enriched mRNA from dissected photophore tissue (Appendix 1-table 1). Photophores were dissected from the abdomens of adult P. pyralis males (Collected MMNJ, July 2015) by Dr. Adam South (Harvard School of Public Health), using three individuals per biological replicate. These tissues and libraries were co-prepared and sequenced with other previously published libraries (full library preparation and sequencing details here [Al-Wathiqui et al., 2016]) at a depth of~10M paired reads per library.
To examine gene expression in larval light organs, we performed RNA-seq on dissected larval light organs. We first extracted total RNA from a pool of six dissected larval photophores from three individuals using the RNeasy Lipid Tissue Mini Kit (QIAGEN) with the optional on-column DNase treatment. The larvae were the same larvae described in Appendix 1.3.2. The total RNA was enriched to mRNA via polyA pulldown and prepared into a paired unstranded Illumina sequencing using the Kapa HyperPrep kit (Kapa Biosystems, USA), and sequenced to a depth of 43M 100 Â 100 paired reads on a HiSeq2500 sequencer (Illumina, USA).
All these data were combined with previously published tissue, sex, and stage-specific libraries (Appendix 1-table 1) for reference-guided transcriptome assembly (Appendix 1.9.3). Strand-specific data was used for de novo transcriptome assembly (Appendix 1.9.2).   One strand-specific de novo transcriptome was produced from all available MMNJ strandspecific reads (WholeMale, WholeFemale, eggs, larvae) and strand-specific reads from SRA (SRR3521424) (Appendix 1-table 1). Reads from these five libraries were pooled (158.6M paired-reads) as input for de novo transcriptome assembly. Transcripts were assembled using Trinity (v2.4.0) (Grabherr et al., 2011) with default parameters except the following: (-SS_lib_type RF -trimmomatic -min_glue 2 min_kmer_cov 2 -jaccard_clip -no_normalize_reads). Gene structures were then predicted from alignment of the de novo transcripts to the Ppyr1.3 genome using the PASA pipeline (v2.1.0) (Haas et al., 2008) with the following steps: first, poly-A tails were trimmed from transcripts using the internal seqclean component; next, transcript accessions were extracted using the accession_extractor.pl component; finally, the trimmed transcripts were aligned to the genome with modified parameters (-aligners blat,gmap -ALT_SPLICE -transcribed_is_aligned_orient -tdn tdn.accs). Using both the blat (v. 36 Â 2) (Kent, 2002) and gmap (v2017-09-11) (Wu and Watanabe, 2005) aligners was required, as an appropriate gene model for Luc2 was not correctly produced using only a single aligner. Importantly, it was also necessary to set (-NUM_BP_PERFECT_SPLICE_BOUNDARY = 0) for the validate_alignments_in_db.dbi step, to ensure transcripts with natural variation near the splice sites were not discarded. Post alignment, potentially spurious transcripts were filtered out using a custom script (Fallon, 2017; copy archived at https://github.com/elifesciencespublications/PASA_expression_filter_2017) that removed extremely lowly-expressed transcripts (<1% of the expression of a given PASA assembly cluster). Expression values used for filtering were calculated from the WholeMale library reads using the Trinity align_and_estimate_abundance.pl utility script. The WholeMale library was selected because it was the highest quality library -strand-specific, low contamination, and many readsthereby increasing the reliability of the transcript quantification. Finally, the PASA pipeline was run again with this filtered transcript set to generate reliable transcript structures. Peptides were predicted from the final transcript structures using Transdecoder (v.5.0.2) (Haas, 2018) with default parameters. Direct coding gene models (DCGMs) were then produced with the Transdecoder 'cdna_alignment_orf_to_genome_orf.pl' utility script with the PASA assembly GFF and transdecoder predicted peptide GFF as input. The unaligned de novo transcriptome assembly is dubbed 'PPYR_Trinity_stranded', whereas the aligned direct coding gene models are dubbed 'Ppyr1.3_Trinity-PASA_stranded-DCGM'.

Reference guided transcriptome assembly
Two reference guided transcriptomes, one strand-specific and one non-strand-specific, were produced from all available P. pyralis RNA-Seq reads (Appendix 1-table 1) using HISAT2 (v2.0.5) (Kim et al., 2015) and StringTie (v1.3.3b) (Pertea et al., 2015). For each library, reads were first mapped to the Ppyr1.3 genome assembly with HISAT2 (parameters: -X 2000 -dta -fr) and then assembled using StringTie with default parameters except use of '-rf' for the strand-specific libraries. The resulting library-specific assemblies were then merged into a final assembly using StringTie (-merge), one for the strand-specific and one for the nonstrand specific libraries, producing two final assemblies. For each final assembly, a transcript fasta file was produced and peptides predicted using Transdecoder with default parameters. Then, the StringTie. GTFs were converted to GFF format with the Transdecoder 'gtf_to_alignment_gff3.pl' utility script and direct coding gene models (DCGMs) were produced with the Transdecoder 'cdna_alignment_orf_to_genome_orf.pl' utility script, with the StringTie GFF and transdecoder predicted peptide GFF as input. The final GFFs were validated and sorted with genometools (v1.5.9) with parameters (parameters: gff3 -tidy -sort -retainids), and then sorted again for IGV format with igvtools (parameters: sort). The aligned direct coding gene models for the stranded and unstranded reference guided transcriptomes are dubbed 'Ppyr1.3_Stringtie_stranded-DCGM' and 'Ppyr1.3_Stringtie_unstranded-DCGM'. -l endopterygota_odb9 -long -species tribolium2012). Next, preliminary gene models for prediction training were produced by the alignment of the P. pyralis de novo transcriptome to Ppyr1.2 with the MAKER pipeline (v3.0.0b) (Holt and Yandell, 2011) in 'est2genome' mode. Preliminary gene models were used to train SNAP (v2006-07-28) (Korf, 2004) following the MAKER instructions (MAKER Tutorial for GMOD Online Training, 2014). Augustus and SNAP gene predictions of Ppyr1.3 were then produced through the MAKER pipeline, with hints derived from MAKER blastx/exonerate mediated protein alignments of peptides from Drosophila melanogaster (NCBI GCF_000001215.4_Release_6_plus_ISO1_MT_protein.faa), Tribolium castaneum (NCBI GCF_000002335.3_Tcas5.2_protein), and Aquatica lateralis (AlatOGS1.0; this report), and MAKER blastn/exonerate transcript alignments of the P. pyralis de novo transcriptome. These ab initio coding gene models are dubbed 'Ppyr1.3_abinitio_Augustus-SNAP-MAKER-GMs.gff3' We then integrated the ab initio predictions with our de novo and reference guided direct coding gene models, using EVM. A variety of evidence sources, and EVM evidence weights were empirically tested and evaluated using a combination of inspection of known gene models (e.g. Luc1/Luc2), and the BUSCO score of the geneset. In the final version, six sources of evidence were used for EVM: de novo transcriptome direct coding gene models (Ppyr1.3_Trinity-PASA_stranded-DCGM; weight = 11), protein alignments (D. melanogaster, T. castaneum, A. lateralis; weight = 8), GMAP and BLAT alignments of de novo transcriptome (via PASA; weight = 5), reference guided transcriptome direct coding gene models (Ppyr1.3_Stringtie_stranded-DCGM; weight = 3), Augustus and SNAP ab initio gene models (via MAKER; weight = 2). A custom script (Fallon, 2018a; copy archived at https:// github.com/elifesciences-publications/maker_gff_to_evm_gff_2017) was necessary to convert MAKER GFF format to an EVM compatible GFF format.
Lastly, gene models for luciferase homologs, P450s (Appendix 1.10.1), and de novo methyltransferases (DNMTs) which were fragmented or were incorrect (e.g. fusions of adjacent genes) were manually corrected based on the evidence of the de novo and reference guided direct coding gene models. Manual correction was performed by performing TBLASTN searches with known good genes from these gene families within SequencerServer(v1.10.11) (Priyam et al., 2015), converting the TBLASTN results to gff3 format with a custom script (Fallon, 2018b; copy archived at https://github.com/ elifesciences-publications/firefly_genomes_general_scripts), and viewing these alignments alongside the alternative direct coding gene models (Appendix 1.9.2; 1.9.3) in Integrative Genomics Viewer(v2.4.8) (Thorvaldsdó ttir et al., 2013). The official gene set models gff3 file was manually modified in accordance with the evidence from the direct gene models. Different revision numbers of the official geneset (e.g. PPYR_OGS1.0, PPYR_OGS1.1) represent the improvement of the geneset over time due to these continuing manual gene annotations.

P450 annotation
Translated de novo transcripts were formatted to be BLAST searchable with NCBI's standalone software. The peptides were searched with 58 representative insect P450s in a batch BLAST (evalue = 10). The query set was chosen to cover the diversity of insect P450s. The top 100 hits from each search were retained. The resulting 5837 hit IDs were filtered to remove duplicates, leaving 472 unique hits. To reduce redundancy due to different isoforms, the Trinity transcript IDs (style DNXXX_cX_gX_iX) were filtered down to the 'DN' level, resulting in 136 unique IDs. All peptides with these IDs were retrieved and clustered with CD-Hit (v4.5.4) (Li and Godzik, 2006) to 99% identity to remove short overlapping peptides. These 535 protein sequences were batch BLAST compared to a database of all named insect P450s to identify best hits. False positives were removed and about 30 fungal sequences were removed. These fungal sequences could potentially be from endosymbiotic fungi in the gut. Overlapping sequences were combined and the transcriptome sequences were BLAST searched against the P. pyralis genome assembly to fill gaps and extend the sequences to the ends of the genes were possible. This approach was very helpful with the CYP4G gene cluster, allowing fragments to be assembled into whole sequences. When a new genome assembly and geneset became available, the P450s were compared to the integrated gene models in PPYR_OGS1.0. Some hybrid sequences were corrected. The final set contains 170 named cytochrome P450 sequences (166 genes, two pseudogenes).
The cytochrome P450s in insects belong to four established clans CYP2, CYP3, CYP4 and Mito (Appendix 1- figure 13). P. pyralis has about twice as many P450s as Drosophila melanogaster (86 genes, four pseudogenes) and slightly more than the red flour beetle Tribolium castaneum (137 genes, 10 pseudogenes). Pseudogenes were determined by a lack of conserved sites common to all P450s. The CYP3 clan is the largest, mostly due to three families: CYP9 (40 sequences), CYP6 (36 sequences) and CYP345 (18 sequences). Insects have few conserved sequences across species. These include the halloween genes for 20hydroxyecdysone synthesis and metabolism CYP302A1, CYP306A1, CYP307A2, CYP314A1 and CYP315A1 (Rewitz et al., 2007) in the CYP2 and Mito clans. The CYP4G subfamily makes a hydrocarbon waterproof coating for the exoskeleton (Helvig et al., 2004). Additional conserved P450s are CYP15A1 (juvenile hormone [Helvig et al., 2004]) and CYP18A1 (20-hydroxyecdysone degradation [Guittard et al., 2011]) in the CYP2 clan. Most of the other P450s are limited to a narrower phylogenetic range. Many are unique to a single genus, although this may change as more sampling is done. It is common for P450s to expand into gene blooms (Sezutsu et al., 2013).   2017). RNA sequence reads were first de novo assembled using Trinity v2.4.0 (Grabherr et al., 2011) with default parameters. Resulting transcriptomes were assessed for similarity to known viral sequences by TBLASTN searches (max e-value = 1 Â 10 À5 ) using as a probe the complete predicted non redundant viral Refseq proteins retrieved from NCBI (date accessed: 15th June 2017). Significant hits were explored manually and redundant contigs discarded. False-positives were eliminated by comparing candidate viral contigs to the entire non-redundant nucleotide (nt) and protein (nr) database to remove false-positives. Candidate virus genome segment sequences were curated by iterative mapping of reads using Bowtie 2 (v2.3.2) (Langmead and Salzberg, 2012). Special attention was taken with the segments' terminis -an arbitrary cut off of 10x coverage was used as threshold to support terminal base calls. The complementarity and folded structure of untranslated ends, as would be expected for members of the Orthomyxoviridae, was assessed by Mfold 2.3 (Zuker, 2003). Further, conserved UTR sequences were identified using ClustalW2 (Larkin et al., 2007) (support of >65% required to call a base). To identify/rule out additional segments of no homology to the closely associated viruses we used diverse in silico approaches based on RNA levels including: the sequencing depth of the transcript, predicted gene product structure, or conserved genome termini, and significant coexpression with the remaining viral segments.
After these filtering steps, putative viral sequences were annotated manually. First, potential open-reading frames (ORF) were predicted by ORFfinder (Wheeler et al., 2003) and manually inspected by comparing predicted ORFS to those from the closest-related reference virus genome sequence. Then, translated ORFs were blasted against the nonredundant protein sequences NR database and best hits were retrieved. Predicted ORF protein sequences were also subjected to a domain-based Blast search against the Conserved Domain Database (CDD) (v3.16) (Marchler-Bauer et al., 2017) and integrated with SMART (Letunic and Bork, 2018), Pfam (Finn et al., 2016), and PROSITE (Sigrist et al., 2002) results to characterize the functional domains. Secondary structure was predicted with Garnier as implemented in EMBOSS (v6.6) (Rice et al., 2000), signal and membrane cues were assessed with SignalP (v4.1) (Petersen et al., 2011), and transmembrane topology and signal peptides were predicted by Phobius (Käll et al., 2004). Finally, the potential functions of predicted ORF products were explored using these annotations as well as similarity to viral proteins of known function.
To characterize Orthomyxoviridae viral diversity in P. pyralis in relation to known viruses, predicted P. pyralis viral proteins were used as probes in TBLASTN (max e-value = 1Â10 À5 ) searches of the complete 2754 Transcriptome Shotgun Assembly (TSA) projects on NCBI (date accessed: 15th June 2017). Significant hits were retrieved and the target TSA projects further explored with the complete Orthomyxoviridae refseq collection to assess the presence of additional similar viral segments. Obtained transcripts were extended/curated using the SRA associated libraries for each TSA hit and then the curated virus sequences were characterized and annotated as described above.
To identify P. pyralis viruses to family/genus/species, amino acid sequences of the predicted viral polymerases, specifically the PB1 subunit, were used for phylogenetic analyses with viruses of known taxonomy. To do this, multiple sequence alignment were generated using MAFFT (v7.310) (Katoh and Standley, 2013) and unrooted maximumlikelihood phylogenetic trees were constructed using FastTree (Price et al., 2010) with standard parameters. FastTree accounted for variable rates of evolution across sites by assigning each site to one of 20 categories, with the rates geometrically spaced from 0.05 to 20, and set each site to its most likely rate category using a Bayesian approach with a gamma prior. Support for individual nodes was assessed using an approximate likelihood ratio test with the Shimodaira-Hasegawa-like procedure. Tree topology, support values and substitutions per site were based on 1000 tree resamples.
To facilitate taxonomic identification, we complemented BLASTP data with two levels of phylogenetic insights: (i) Trees based on the complete refseq collection of ssRNA (-) viruses which permitted a conclusive assignment at the virus family level. (ii) Phylogenetic trees based on reported, proposed, and discovered Orthomyxoviridae viruses that allowed tentative species demarcation and genera postulation. PB1-based trees were complemented independently with phylogenetic studies derived from amino acids of predicted nucleoproteins, hemagglutinin protein, PB2 protein, and PA protein which supported species, genera and family demarcation based on solely on PB1, the standard in Orthomyxoviridae. In addition, sequence similarity of concatenated gene products of International Committee on Taxonomy of Viruses (ICTV) allowed demarcation to species and firefly viruses were assessed by Circoletto diagrams (Darzentas, 2010) (e-value = 1e-2). Where definitive identification was not easily assessed, protein Motif signatures were determined by identification of region of high identity between divergent virus species, visualized by Sequence Logo (Crooks et al., 2004), and contrasted with related literature. Heterotrimeric viral polymerase 3D structure prediction was generated with the SWISS-MODEL automated protein structure homology-modeling server (Biasini et al., 2014) with the best fit template 4WSB: the crystal structure of Influenza A virus 4WSB. Predicted structures were visualized in UCSF Chimera (Pettersen et al., 2004) and Needleman-Wunsch sequence alignments from structural superposition of proteins were generated by MatchMaker and the Match->Align Chimera tool. Alternatively, 3D structures were visualized in PyMOL (v1.8.6.0; Schrodinger). Viral RNA levels in the transcriptome sequences were also examined. Virus transcripts RNA levels were obtained by mapping the corresponding raw SRA FASTQ read pairs using either Bowtie2 (Langmead and Salzberg, 2012)  All curation, phylogeny construction, and visualization were conducted in Geneious 8.1.9 (Biomatters, Ltd.). Animal silhouettes in Appendix 5- figure 2 were developed based on non-copyrighted public domain images. Figure compositions were assembled using Photoshop CS5 (Adobe). Bar graphs were generated with Excel 2007 software (Microsoft). RNA levels normalized as mapped transcripts per million per library were visualized using Shinyheatmap (Khomtchouk et al., 2017).
Finally, to identify endogenous viral-like elements, tentative virus detections and the viral refseq collection were contrasted to the P. pyralis genome assembly Ppyr1.2 by BLASTX searches (e-value = 1e-6) and inspected by hand. Then 15 Kbp genome flanking regions were retrieved and annotated. Lastly, transposable elements (TEs) were determined by the presence of characteristic conserved domains (e.g. RNASE_H, RETROTRANSPOSON, INTEGRASE) on predicted gene products and/or significant best BLASTP hits to reported TEs (e-value <1e-10).

Repeat annotation
Repeat prediction for P. pyralis was performed de novo using RepeatModeler (v1.0.9) (Smit and Hubley, 2017) and MITE-Hunter (v11-2011) (Han and Wessler, 2010). RepeatModeler uses RECON (Bao and Eddy, 2002) and RepeatScout (Price et al., 2005) to predict interspersed repeats, and then refines and classifies the consensus repeat models to build a repeat library. MITE-Hunter detects candidate MITEs (miniature inverted-repeat transposable elements) by scanning the assembly for terminal inverted repeats and target site duplications < 2 kb apart. To identify tandem repeats, we also ran Tandem Repeat Finder (v4.09; parameters: 2 7 7 80 10) (Benson, 1999), and added repeats whose repeat block length was >5 kb to the repeat library annotated as 'complex tandem repeat'. The RepeatModeler and MITE-Hunter libraries were combined and classified using RepeatClassifier (RepeatModeler 1.0.9 distribution) (Smit and Hubley, 2017). The complex repeats identified by Tandem Repeat Finder were added to this classified list to create the final library of 3118 repeats. This repeat library is dubbed the P. pyralis Official Repeat Library 1.0 (PPYR_ORL1.0).

Telomere FISH analysis
We synthesized a 5' fluorescein-tagged (TTAGG) 5 oligo probe (FAM; Integrated DNA Technologies) for fluorescence in situ hybridization (FISH). We conducted FISH on squashed larval tissues according to previously published methods (Larracuente and Ferree, 2015), with some modification. Briefly, we dissected larvae in 1X PBS and treated tissues with a hypotonic solution (0.5% Sodium citrate) for 7 min. We transferred treated larval tissues to 45% acetic acid for 30 s, fixed in 2.5% paraformaldehyde in 45% acetic acid for 10 min, squashed, and dehydrated in 100% ethanol. We treated dehydrated slides with detergent (1% SDS), dehydrated again in ethanol, and then stored until hybridization. We hybridized slides with probe overnight at 30˚C, washed in 4X SSCT and 0.1X SSC at 30˚C for 15 min per wash. Slides were mounted in VectaShield with DAPI (Vector Laboratories), visualized on a Leica DM5500 upright fluorescence microscope at 100X, imaged with a Hamamatsu Orca R2 CCD camera. Images were captured and analyzed using Leica's LAX software. Aquatica lateralis (Motschulsky, 1860) (Japanese name, Heike-botaru / ) is one of the most common and popular luminous insects in mainland Japan. This species is a member of the subfamily Luciolinae and had long belonged in the genus Luciola, but was recently moved to the new genus Aquatica with some other Asian aquatic fireflies (Fu et al., 2010).
The life cycle of A. lateralis is usually 1 year. Aquatic larva possesses a pair of outer gills on each abdominal segment and live in still or slow streams near rice paddies, wetlands and ponds. Larvae mainly feed on freshwater snails. They pupate in a mud cocoon under the soil near the water. Adults emerge in early to end of summer. While both males and females are full-winged and can fly, there is sexual dimorphism in adult size: the body length is about 9 mm in males and 12 mm in females (Ohba, 2004).
Like other firefly larvae, A. lateralis larvae are bioluminescent. Larvae possess a pair of lanterns at the dorsal margin of the abdominal segment 8. Adults are also luminescent and possess lanterns at true abdominal segments 6 and 7 in males and at segment six in females (Branham and Wenzel, 2003;Ohba, 2004;Kanda, 1935). The adult is dusk active. Male adults flash yellow-green for about 1.0 s in duration every 0.5-1.0 s while flying~1 m above the ground. Female adults, located on low grass, respond to the male signal with flashes of 1-2 s in duration every 3-6 s. Males immediately approach females and copulate on the grass (Ohba, 2004;Ohba, 1983). Like many other fireflies, A. lateralis is likely toxic: both adults and larvae emit an unpleasant smell when disturbed and both invertebrate (dragonfly) and vertebrate (goby) predators vomit up the larva after biting (Ohba and Hidaka, 2002). A. lateralis larvae have eversible glands on each of the eight abdominal segments (Fu et al., 2010). The contents of the eversible glands is perhaps similar to that reported for A. leii (Fu et al., 2007).

Species distribution
The geographical range of A. lateralis includes Siberia, Northeast China, Kuril Isls, Korea, and Japan (Hokkaido, Honshu, Shikoku, Kyushu, Tsushima Isls.) (Kawashima et al., 2003). Natural habitats of these Japanese fireflies have been gradually destroyed through human activity, and currently these species can be regarded as 'flagship species' for conservation (Higuchi, 1996). For example, in 2017, Japanese Ministry of Environment began efforts to protect the population of A. lateralis in the Imperial Palace, Tokyo, where 3000 larvae cultured in an aquarium were released in the pond beside the Palace (Imperial Palace Outer Garden Management Office, 2017).

Specimen collection
Individuals used for genome sequencing, RNA sequencing, and LC-HRAM-MS were derived from a small population of laboratory-reared fireflies. This population was established from a few individuals collected from rice paddy in Kanagawa Prefecture of Japan in 1989 and 1990 (Ikeya, 2016) by Mr. Haruyoshi Ikeya, a highschool teacher in Yokohama, Japan. Mr. Ikeya collected adult A. lateralis specimens from their natural habitat in Yokohama and has propagated them for over 25 years (~25 generations) in a laboratory aquarium without any addition of wild individuals. This population has since been propagated in the laboratory of YO and JKW, and is dubbed the 'Ikeya-Y90' cultivar. Because of the small number of individuals used to establish the population and the number of generations of propagation, this population likely represents a partially inbred strain. Larvae were kept in aquarium at 19-21˚C and fed using freshwater snails (Physella acuta and Indoplanorbis exustus). Under laboratory rearing conditions, the life cycle is reduced to 7-8 months. The original habitat of this strain has been destroyed and the wild population which led to the laboratory strain is now extinct.

Karyotype and genome size
Unlike P. pyralis, the karyotype of A. lateralis is reported to be 2n = 16 with XY sex determination (male, 14A + XY; female, 14A + XX) (Inoue and Yamamoto, 1987). The Y chromosome is much smaller than X chromosome, and the typical behaviors of XY chromosomes, such as partial conjugation of X/Y at the first meiotic metaphase and a separation delay of X/Y at the first meiotic anaphase, were observed in testis cells (Inoue and Yamamoto, 1987).
We determined the genome size of A. lateralis using flow cytometry-mediated calibratedfluorimetry of DNA content with propidium iodide stained nuclei. First, the head+prothorax of a single pupal female (gender identified by morphological differences in abdominal segment VIII) was homogenized in 100 mL PBS. These tissues were chosen to avoid the ovary tissue. Once homogenized, 900 mL PBS, 1 mL Triton X-100 (Sigma-Aldrich), and 4 mL 100 mg/mL RNase A (QIAGEN) were added. The homogenate was incubated at 4˚C for 15 min, filtered with a 30 mm Cell Tries filter (Sysmex), and further diluted with 1 mL PBS. 20 mL of 0.5 mg/mL propidium iodide was added to the mixture and then average fluorescence of the 2C nuclei determined with a SH-800 flow cytometer (Sony, Japan). Three technical replicates of this sample were performed. Independent runs for extracted Aphid nuclei (Acyrthosiphon pisum; 517 Mbp), and fruit fly nuclei (Drosophila melanogaster; 175 Mbp) were performed as calibration standards. Genome size was estimated at 940 Mbp ±1.4 (S.D.; technical replicates = 3).
Genome size inference via Kmer spectral analysis estimated a genome size of 772 Mbp (Appendix 2- figure 1).

Genomic sequencing and assembly
Genomic DNA was extracted from the whole body of a single laboratory-reared A. lateralis adult female (c.v. Ikeya-Y90) using the QIAamp Kit (Qiagen). Purified DNA was fragmented with a Covaris S2 sonicator (Covaris, Woburn, MA), size-selected with a Pippin Prep (Sage Science, Beverly, MA), and then used to create two paired-end libraries using the TruSeq Nano Sample Preparation Kit (Illumina) with insert sizes of~200 and~800 bp. These libraries were sequenced on an Illumina HiSeq1500 using a 125 Â 125 paired-end sequencing protocol. Mate-pair libraries of 2-20 Kb with a peak at~5 Kb were created from the same genomic DNA using the Nextera Mate Pair Sample Preparation Kit (FC-132-1001, Illumina), and sequenced on HiSeq 1500 using a 100 Â 100 paired-end sequencing protocol at the NIBB Functional Genomics Facility (Aichi, Japan). In total, 133.3 Gb of sequence (159x) was generated.
Reads were assembled using ALLPATHS-LG (build# 48546) (Gnerre et al., 2011), with default parameters and the 'HAPLOIDIFY = True' option. Scaffolds were filtered to remove non-firefly contaminant sequences using blobtools (Laetsch and Blaxter, 2017) , 2017). len = inferred haploid genome length, uniq = percentage nonrepetitive sequence, het = overall rate of genome heterozygosity, kcov = mean kmer coverage for heterozygous bases, err = error rate of the reads, dup: average rate of read duplications. These results are consistent when considering the possible systematic error of kmer spectral analysis and flow cytometry genome size estimates. The heterozygosity is lower than that measured for P. pyralis, possibly reflecting the long-term laboratory rearing in reduced population sizes of A. lateralis strain Ikeya-Y90.

Taxonomic annotation filtering
Potential contaminants in Alat1.2 were identified using the blobtools toolset (v1.0) (Laetsch and Blaxter, 2017). First, scaffolds were compared to known sequences by performing a blastn (v2.5.0+) nucleotide sequence similarity search against the NCBI nt database and a diamond (v0.9.10) (Buchfink et al., 2015) translated nucleotide sequence similarity search against the of Uniprot reference proteomes (July 2017). Using this similarity information, scaffolds were annotated with blobtools (parameters '-x bestsumorder'). We also inspected the read coverage by mapping the paired-end reads (FFGPE_PE200) on the genome using bowtie2. A tab delimited text file containing the results of this blobtools annotation are available on FigShare (DOI: 10.6084/m9.figshare.5688928). The contigs derived from potential contaminants and/or poor quality contigs were then removed: contigs with higher %GC (>50%) with bacterial hits or no database hits and showing low read coverage (<30 x) (see Appendix 2- figure 2). This process removed 1925 scaffolds (1.17

RNA-extraction, library preparation and sequencing
In order to capture transcripts from diverse life-stages and tissues, non-stranded RNA-Seq libraries were prepared from fresh specimens of nine life stages/sexes/tissues (eggs, fifth (the last) instar larvae, both sex of pupae, adult male head, male abdomen (prothorax-to-fifth segment), male lantern, adult female head, and female lantern (Appendix 2-table 1). Live specimens were anesthetized on ice and dissected during the day. The lantern tissue was dissected from the abdomen and contains the cuticle, photocyte layer and reflector layer. For eggs, larvae, and pupae, total RNA was extracted using the RNeasy Mini Kit (QIAGEN) with the optional on-column DNase treatment. For adult specimens, total RNA was extracted using TRIzol reagent (Invitrogen) to avoid contamination of pigments and uric acid. These were then treated with DNase in solution and then cleaned using a RNeasy Mini kit.

Official coding geneset annotation (AQULA_OGS1.0)
A protein-coding gene reference set for A. lateralis was generated by Evidence Modeler (v1.1.1) using both aligned transcripts and aligned proteins. For transcripts, we combined reference guided and de novo transcriptome assembly approaches. Notably, these reference guided and de novo transcriptome assembly approaches differed from the current de novo (Appendix 2.7.1) and reference guided (Appendix 2.7.2) transcriptome assembly approaches.
In the reference-guided approach applied here, RNA-Seq reads were mapped to the genome assembly with TopHat and assembled into transcripts with Cufflinks (parameters:min-intron-length 30) (Trapnell et al., 2010). The Cufflinks transcripts were subjected to the TransDecoder program to extract ORFs. In the de novo transcriptome approach applied here, RNA-seq reads were assembled de novo by Trinity and ORFs were predicted using TransDecoder. We used CD-HIT-EST (Li and Godzik, 2006) to reduce the redundancy of the predicted ORFs. The ORF sequences were mapped to the genome using Exonerate in est2genome mode for splice-aware alignment. We processed homology evidence at the protein level using the reference proteomes of D. melanogaster and T. castaneum. These reference proteins were split-mapped to the A. lateralis genome in two steps: first with BLASTX to find approximate loci, and then with Exonerate in protein2genome mode to obtain more refined alignments. These gene models derived from multiple evidence were merged by the EVM program to obtain the reference annotation for the genomes. We also predicted ab initio gene models using Augustus, but we didn't include Augustus models for the EVM integration because our preliminary analysis showed the ab initio gene models had no positive impact on gene prediction. Lastly, gene models for luciferase homologs, P450s, and de novo methyltransferases (DNMTs) which were fragmented or were incorrect (e.g. fusions of adjacent genes) were manually corrected based on the evidence of the de novo and reference guided direct coding gene models. Manual correction was performed by performing TBLASTN searches with known good genes from these gene families within SequencerServer(v1.10.11) (Priyam et al., 2015), converting the TBLASTN results to gff3 format with a custom script (Fallon, 2018b), and viewing these alignments alongside the alternative direct coding gene models (Appendix 2.7.1; 2.7.2) in Integrative Genomics Viewer(v2.4.8) (Thorvaldsdó ttir et al., 2013). The official gene set. gff3 file was manually modified in accordance with the alternative gene models. Different revision numbers of the official geneset (e.g. AQULA_OGS1.0, AQULA_OGS1.1) represent the improvement of the geneset over time due to these continuing manual gene annotations.

Repeat annotation
A de novo species-specific repeat library for A. lateralis was constructed using RepeatModeler (v1.0.9), and Tandem Repeat Finder (v4.09; settings: 2 7 7 80 10) (Benson, 1999). Only tandem repeats from Tandem Repeat Finder with a repeat block length >5 kb (annotated as 'complex tandem repeat') were added to the RepeatModeler library. This process yielded a final library of 1695 interspersed repeats. We then used this library and RepeatMasker (v4.0.5) (Smit et al., 2015) to identify and mask interspersed and tandem repeats in the genome assembly. This repeat library is dubbed the Aquatica lateralis Official Repeat Library 1.0 (AQULA_ORL1.0).  , which are widespread throughout the globe. Unlike in fireflies, where bioluminescence is universal, only~200 described elaterid species are luminous. These luminous species are recorded only from tropical and subtropical regions of Americas and some small Melanesian islands, such as Fiji and Vanuatu (Costa, 1975;. For instance, the tropical American Pyrophorus noctilucus is considered the largest (~30 mm) and brightest bioluminescent insect (Harvey and Stevens, 1928;Levy, 1998). All luminous species are closely related -luminous click beetles belong to the tribes Pyrophorini and Euplinthini (Costa, 1975;Arias-Bohart, 2015) of the subfamily Agrypninae, with the single exception of Campyloxenus pyrothorax (Chile) in the related subfamily Campyloxeninae (Stibick, 1979). The luminescence of a pair of pronotal 'light organs' of the adult Balgus schnusei (Costa, 1984), a species that has now been assigned to the Thylacosterninae of the Elateridae , has not been confirmed by later observation. This near-monophyly of bioluminescent elaterid taxa is supported by both morphological (Douglas, 2011) and molecular phylogenetic analysis (Sagegami- Oba et al., 2007;Oba and Sagegami-Oba, 2007;Kundrata and Bocak, 2011), although early morphological phylogenies were inconsistent (Stibick, 1979;Hyslop, 1917;Ohira, 1962;Dolin, 1978;Ohira, 2013). This suggests a single origin of bioluminescence in this family. The genus Ignelater was established by Costa in 1975 and I. luminosus was included in this genus (Costa, 1975). Often this species is called Pyrophorus luminosus as an 'auctorum', a name used to describe a variety of taxa (Johnson, 2002). This use of 'Pyrophorus' as an auctorum may be due to the heightened difficulty of classifying Elateridae (Costa, 1975). The genus Ignelater is characterized by the presence of both dorsal and ventral photophores (Costa, 1975;Rosa, 2007). An unreviewed report suggested that the adult I. luminosus has a ventral light organ only in males (Reyes and Lee, 2010). Phylogenetic analyses based on the morphological characters suggested that the genera Ignelater and Photophorus (which contain only two species from Fiji and Vanuatu) are the most closely related genera in the tribe Pyrophorini (Rosa, 2007). The earliest fossil of an Elateridae species was recorded from the Middle Jurassic of Inner Mongolia, China (Chang et al., 2009). McKenna and Farrell suggested that, based on molecular analyses, the family Elateridae originated in the Early Cretaceous (130 Mya) (McKenna and Farrell, 2009). It is expected that many recent genera in Elateroidea were established by the Early Tertiary (<65 Mya) (Grimaldi and Engel, 2005).
The exact function of bioluminescence across different life stages remains unknown for many luminous elaterid species. Bioluminescent elaterid beetles typically have two paired lanterns on the dorsal surface of the prothorax, and a single lantern on the ventral abdomen, which is only exposed during flight. Several bioluminescent Elateridae produce different colored luminescence from their prothorax and abdominal lanterns (Oba et al., 2010a;Feder and Velez, 2009). Harvey reported that there was not a marked difference in the luminescence color of the dorsal and ventral lanterns of Puerto Rican I. luminosus (Harvey, 1952). Like fireflies, elaterid larvae often produce light, with the glowing termite mounds of Brazil that contain the predatory larvae of Pyrearinus termitilluminans being a striking example (Costa and Vanin, 2010). A description of the anatomy of the larval light organ of Pyrophorus is provided by (Harvey, 1952), and a more modern photograph of the larval light organ is provided by (Bechara and Stevani, 2018). Like other bioluminescent elaterid larvae, I. luminosus larvae produce a diffuse light from their prothorax, however they are only luminous when disturbed (Wolcott, 1948). I. luminosus larvae are subterranean predators and are an enthusiastic predator of the white grub (Ancylonycha spp.), reportedly consuming 50 + to reach full size (Wolcott, 1950). Adult I. luminosus are luminescent and a bioluminescent courtship behavior was described in an unreviewed study (Kretsch, 2000). Reportedly, males search during flight with their prothorax lanterns illuminated steadily, while females stay on the ground modulating the intensity of their prothorax lanterns in~2 s intervals. Once a female is observed, the prothorax lanterns of the male go dark, the ventral lantern becomes illuminated, and the male approaches the female via a circular search pattern. Mating is brief, reportedly taking only 5 seconds.
Unlike fireflies, bioluminescent elaterid species are not known to have potent chemical defenses. For example, the Jamaican bioluminescent elaterid beetle Pyrophorus plagiophthalmus, does not appear to be strongly unpalatable, as bats were observed to regularly capture the beetles during their flying bioluminescent displays (Vélez, 2006). A defense role for I. luminosus luminescence to startle predators is possible.

Species distribution
I. luminosus is often considered to be endemic to Puerto Rico (Virkki et al., 1984); however, the genus Ignelater is reported in Florida (USA), Vera Cruz (Mexico), the Bahamas, Cuba, Isla de la Juventud, Hispaniola (Haiti + Dominican Republic), Puerto Rico, and the Lesser Antilles (Costa, 1975). Similarly, I. luminosus itself has been reported on the island of Hispaniola (Kretsch, 2000;Perez-Gelabert, 2008), indicating I. luminosus is not restricted to Puerto Rico. This geographic distribution of Ignelater suggests that Puerto Rico may contain multiple Ignelater species and, given the difficulty of distinguishing different species of bioluminescent Elateridae by morphological characters, a definitive species distribution for I. luminosus cannot be stated, other than this species is seemingly not strictly endemic to Puerto Rico.

Collection
I. luminosus (Illiger, 1807) adult specimens were collected from private land in Mayagü ez, Puerto Rico (18˚13 ' 12.1974' N, 67˚6' 31.6866' W) with permission of the landowner by Dr. David Jenkins (USDA-ARS). Individuals were captured at night on April 20th and April 28th 2015 during flight on the basis of light production. The I. luminosus specimens were frozen in a À80˚C freezer, lyophilized, shipped to the laboratory (MIT) on dry ice, and stored at À80˚C. Full collection metadata is available from the NCBI BioSample records of these specimens (NCBI Bioproject PRJNA418169). Identification to species was performed by comparing antenna and dissected genitalia morphology to published keys (Costa, 1975;Rosa, 2007;Rosa, 2010) (Appendix 3- figure 1). All inspected specimens were male (3/3). Specimens collected at the same time, but not those used for genitalial dissection, were used for sequencing. Although the genitalia morphology of the sequenced specimens was not inspected to confirm their sex, sequenced specimens were inferred to be male, based on the fact that female bioluminescent elaterid beetles are rarely seen in flight (Personal communication: S. Velez) and the dissected specimens collected in the same batch as the sequenced specimens were confirmed to be male.

I. luminosus aedeagus (male genitalia). (A) Dorsal and (B) ventral view of
an Ignelater luminosus aedeagus, dissected from the same batch of specimens used for linkedread sequencing and genome assembly. The species identity of this specimen was confirmed as I. luminosus by comparison of the aedeagus to the keys of Costa and Rosa (Costa, 1975;Rosa, 2007;Rosa, 2010

Genomic sequencing and assembly
HMW DNA (25 mg) was extracted from a single male specimen of I. luminosus using a 100/G Genomic Tip with the Genomic buffers kit (Qiagen, USA). The I. luminosus specimen was first washed with 95% ethanol, and DNA was extracted following the manufacturer's protocol, with the exception of the final precipitation step, where HMW DNA was pelleted with 40 mg RNA grade glycogen (Thermo Scientific, USA) and centrifugation (3000 x g, 30 min, 4˚C) instead of spooling on a glass rod. HMW DNA was sent on dry-ice to the Hudson Alpha Institute of Biotechnology Genomic Services Lab (HAIB-GSL), where pulsed-field-gel-electrophoresis (PFGE) quality control and 10x Genomics Chromium Genome v1 library construction was performed. PFGE quality control indicated the mean size of the input DNA was >35 kbp+. The assembly was exported to FASTA format using Supernova mkoutput (parameters: -style=pseudohap), and modified by taxonomic annotation filtering (Appendix 3.5.2) and polishing (Appendix 3.5.3) to form Ilumi1.1. A Supernova (v2.0.0) assembly was also produced from combined HiSeqX and HiSeq2500 reads, but on a brief inspection the quality was equivalent to Ilumi1.1, so the new assembly was not used for further analyses. Manual longread based scaffolding was then applied to produce a final assembly Ilumi1.2 (Appendix 3.5.4).  , 2017). Before analysis, 10x Chromium barcodes were trimmed off Read1 using cutadapt (v1.8; parameters: -u 23) (Martin, 2011). vlen = inferred haploid genome length, uniq = percentage non-repetitive sequence, het = overall rate of genome heterozygosity, kcov = mean kmer coverage for heterozygous bases, err = error rate of the reads, dup: average rate of read duplications. These results are consistent when considering the possible systematic error of kmer spectral analysis and flow cytometry genome size estimates. The heterozygosity is higher than that measured for P. pyralis and A. lateralis. The read error rate for this library is also significantly higher than the P. pyralis and A. lateralis results, possibly highlighting the difference in raw read error rate between HiSeq2500 and HiSeqX sequencing, or is possibly an artifact of the Chromium library.

Taxonomic annotation filtering
We sought to systematically remove assembled non-elaterid contaminant sequence from Ilumi1.0. Using the blobtools toolset (v1.0.1), (Laetsch and Blaxter, 2017), we taxonomically annotated our scaffolds by performing a blastn (v2.6.0+) nucleotide sequence similarity search against the NCBI nt database, and a diamond (v0.9.10.111) (Buchfink et al., 2015) translated nucleotide sequence similarity search against the of Uniprot reference proteomes (July 2017). Using this similarity information, we taxonomically annotated the scaffolds with blobtools using parameters '-x bestsumorder -rank phylum' (Appendix 3- figure 3). A tab delimited text file containing the results of this blobtools annotation is available on FigShare (DOI: 10. 6084/m9.figshare.5688952). We then generated the final genome assembly by retaining scaffolds that had coverage >10.0 in the 1610_IlumiHiSeqX library, and did not have a high scoring (score >5000) taxonomic assignment for 'Proteobacteria', followed by polishing indels and gap-filling with Pilon (Appendix 3.5.3). This approach removed 235 scaffolds (330 Kbp), representing 0.2% of the scaffold number and 0.03% of the nucleotides of Ilumi1.0. While filtering the Ilumi1.0 assembly, we noted a large contribution of scaffolds taxonomically annotated as Platyhelminthes (1740 scaffolds; 119.56 Mbp). Upon closer inspection, we found conflicting information as to the most likely taxonomic source of these scaffolds. Diamond searches of these scaffolds had hits in Coleoptera, whereas blastn searches showed these scaffold had confident hits (nucleotide identity >90%, evalue = 0) against the Rat Tapeworm Hymenolepis diminuta genome (NCBI BioProject PRJEB507). Removal of these scaffolds decreased the endopterygota BUSCO score, from C:

Ilumi1.1: Indel polishing
Manual inspection of the initial gene-models for Ilumi1.0 revealed a key luciferase homolog had an unlikely frameshift occurring after a polynucleotide run. Mapping of the 1610_IlumiHiSeqX and 1706_IlumiHiSeq2500 reads (Appendix 4-table 1) with Bowtie2 using parameters (-local), revealed that this indel was not supported by the majority of the data, and that indels were present at a notable frequency after polynucleotide runs. As a greatly increased indel rate after polynucleotide runs (~10% error) is a known systematic error of Illumina sequencing, and has been noted as the major error type in Supernova assemblies (Weisenfeld et al., 2017), we therefore sought to correct these errors globally through the use of Pilon (v1.2.2) (Walker et al., 2014). In order to run Pilon efficiently, we split the taxonomically filtered Ilumi1.0 reference (dubbed Ilumi1.0b; Appendix 3.5.2) using Kirill Kryukov's fasta_splitter.pl script (v0.2.6) (Kryukov, 2017), partitioned the previously mapped 1610_IlumiHiSeqX paired-end reads to these references using samtools, and ran Pilon in parallel on the partitioned reads and records with parameters (-fix gaps,indels -changes -vcfdiploid). The final consensus FASTAs produced by Pilon were merged to produce the polished assembly (Ilumi1.1). Ilumi1.1 (842,900,589 nt; 91,325 scaffolds) was slightly smaller than Ilumi1.0b (845,332,796 nt; 91,325 scaffolds), indicating the gaps filled by Pilon were smaller than their predicted size. The BUSCO score increased modestly after polishing (C:93.3% to C:94.8%), suggesting that indel polishing and gap filling had a net positive effect.

Ilumi1.2: Manual long-read scaffolding
We determined via manual gene-model annotation of Ilumi1.1 (Appendix 3.8), that the second through seventh exon of IlumPACS4 (ILUMI_06433 PA) were present on Ilumi1.1_Scaffold13255, but that the first exon was missing from this scaffold. Targeted tblastn using PangPACS (AB479114.1) (Oba et al., 2010a), the most closely related gene sequence to IlumPACS4, indicated that the most similar region in the I. luminosus genome to the predicted PangPACS first exon was a right-pointing region on Ilumi1.1_Scaffold11560, not captured in any gene model, but downstream of the existing luciferase homolog genes IlumPACS1 and IlumPACS2. We surmised that this region was the correct first exon for IlumPACS4, and that the IlumPACS4 gene model spanned Ilumi1.1_Scaffold13255 and Ilumi1.1_Scaffold11560, and thus that the right edge of Ilumi1.1_Scaffold13255 and the left edge of the reverse complement of Ilumi1.1_Scaffold11560 should be joined. To substantiate this, we performed long-read Oxford Nanopore MinION sequencing at the MIT BioMicroCenter. The HMW DNA used was the same DNA used for Chromium library prep, and had been stored at À80˚C since extraction. Thawing of DNA and size distribution QC on a FEMTO Pulse capillary electrophoresis instrument (Advanced Analytical Technologies Inc, USA) indicated the DNA had a mean size distribution peak of~17 kbp. A 1D Nanopore library was prepared from this DNA using the standard kit and protocol (Part #: SQK-LSK108 kbp read with seven kbp antiparallel alignment to the right edge of Scaffold13255. Inspection of the extension of this read off Scaffold13255 revealed it contained 10 Kbp+ of a nonpalindromic complex tandem repeat DNA with an~100 bp repeat unit (Appendix 3- figure  4). The repeat unit of this complex tandem repeat DNA (Appendix 3-table 1) is annotated in our de novo repeat library construction as 'Ilumi.complex.repeat.1' (Appendix 3.9), and via blastn is clearly interspersed at low copy numbers throughout the Ilumi1.1 genome assembly. Notably, this repeat unit was present the right edge of Ilumi1.1_Scaffold13255, while the reverse complement of this repeat unit was present on the right edge of Ilumi1.1_Scaffold11560, supporting that these scaffolds were adjacent to one another, but the assembly had been broken by this large stretch of tandem repetitive DNA. Although our Nanopore sequencing did not unambiguously span this repetitive element and bridge the two scaffolds, we surmised that this information was sufficient to manually merge these scaffolds (Appendix 3- figure 5). The long Ilumi1.1_Scaffold13255 extending read was adaptor trimmed with porechop (v0.2.3) (Wick, 2018), removing 35 bp from the start of the read. Next, the 3' end of the read which aligned up to the last nucleotide of Ilumi1.1_Scaffold13255 was trimmed. Finally, the remaining read was reverse complemented, and concatenated to the right edge of Ilumi1.1_Scaffold13255. 1337 Ns were concatenated to the right edge of the extended Ilumi1.1_Scaffold13255 to indicate an uncertainty in the repeat copy number, and Ilumi1.1_Scaffold11560 was reverse complemented and concatenated to Ilumi1.1_Scaffold13255 to produce the final version of Ilumi1.2_Scaffold13255 (Appendix 3- figure 5). Further whole genome scaffolding using this Nanopore data and the LINKS pipeline (v1.8.5) (Warren et al., 2015) with parameters (-d 4000,8000,10000,14000,16000,20000 t 2,3,5,9 l 2 -a 0.75) was attempted, but only a single additional pair of scaffolds was merged, so this whole-genome scaffolding was not used further.
repair, an 'A' base was added at the 3'-end of each strand. The Ad153-2B adapters with barcode was ligated to both ends of the end repaired/dA tailed DNA fragments, then amplification by ligation-mediated PCR. Following this, a single strand DNA was separated at a high temperature and then a Splint oligo sequence was used as bridge for DNA cyclization to obtain the final library. Then rolling circle amplification (RCA) was performed to produce DNA Nanoballs (DNBs). The qualified DNBs were loaded into the patterned nanoarrays and the libraries were sequenced as 50 Â 50 bp (PE-50) read through on the BGISEQ-500 platform. Sequencing-derived raw image files were processed by BGISEQ-500 base-calling software with the default parameters, generating the 'raw data' for each sample stored in FASTQ format. This library preparation and sequencing was provided free of charge as an evaluation of the BGISEQ-500 platform.

Transcriptome analysis
Both de novo (Appendix 3.7.1) and reference guided (Appendix 3.7.2) transcriptome assembly approaches using Trinity and Stringtie were used, respectively.

Reference guided transcriptome alignment and assembly
A reference guided transcriptome was produced from all available I. luminosus RNA-seq reads (head + prothorax, mesothorax + metathorax, abdomen -both Illumina and BGISEQ-500) using HISAT2 (v2.0.5) (Kim et al., 2015) and StringTie (v1.3.3b) (Pertea et al., 2015). Reads were first mapped to the I. luminosus draft genome with HISAT2 (parameters: -X 2000 -dta -fr). Then StringTie assemblies were performed on each separate bam file corresponding to the original libraries using default parameters. Finally, the produced GTF files were merged using StringTie (-merge). A transcript fasta file was produced from the StringTie GTF file with the transdecoder 'gtf_genome_to_cdna_fasta.pl' utility script, and peptides were predicted for these transcripts using Transdecoder (v5.0.2) with default parameters. The StringTie GTF was converted to GFF format with the Transdecoder 'gtf_to_alignment_gff3.pl' utility script, and direct coding gene models (DCGMs) were then produced with the Transdecoder 'cdna_alignment_orf_to_genome_orf.pl' utility script, with the StringTie-provided GFF and transdecoder predicted peptide GFF as input. The resulting DCGM GFF3 file was manually lifted over to the Ilumi1.2 assembly. The reference guided transcriptome assembled was dubbed 'ILUMI_Stringtie_unstranded', whereas the aligned direct coding gene models were dubbed 'Ilumi1.2_Stringtie_unstranded-DCGM'
Lastly, gene models for luciferase homologs, P450s, and de novo methyltransferases (DNMTs) which were fragmented or were incorrectly assembled (e.g. adjacent gene fusions) were manually corrected based on the evidence of the de novo and reference guided direct coding gene models (Appendix 3.7.1; 3.7.2). Manual correction was performed by performing TBLASTN searches with known good genes from these gene families within SequencerServer(v1.10.11) (Priyam et al., 2015), converting the TBLASTN results to gff3 format with a custom script (Fallon, 2018b), and viewing these TBLASTN alignments alongside the alternative direct coding gene models and the official geneset in Integrative Genomics Viewer (v2.4.8) (Thorvaldsdó ttir et al., 2013). The official gene set models gff3 file was then manually modified based on the observed evidence. Different revision numbers of the official geneset (e.g. ILUMI_OGS1.0, ILUMI_OGS1.1) represent the improvement of the geneset over time due to these continuing manual gene annotations.

Repeat annotation
A de novo species-specific repeat library for I. luminosus was constructed using RepeatModeler (v1.0.9), and Tandem Repeat Finder (v4.09; settings: 2 7 7 80 10) (Benson, 1999). Only tandem repeats from Tandem Repeat Finder with a repeat block length >5 kb (annotated as 'complex tandem repeat') were added to the RepeatModeler library. This process yielded a final library of 2259 interspersed repeats. We then used this library and RepeatMasker (v4.0.5) (Smit et al., 2015) to identify and mask interspersed and tandem repeats in the genome assembly. This repeat library is dubbed the Ignelater luminosus Official Repeat Library 1.0 (ILUMI_ORL1.0).

Mitochondrial genome assembly and annotation
The mitochondrial genome sequence of I. luminosus was assembled by a targeted subassembly approach. First, Chromium linked-reads were mapped to the previously sequenced mitochondrial genome of the Brazilian elaterid beetle Pyrophorus divergens (NCBI ID: NC_009964.1) (Arnoldi et al., 2007), using Bowtie2 (v2.3.1; parameters: -very-sensitivelocal ) (Langmead et al., 2009). Although these reads still contain the 16 bp Chromium library barcode on read 1 (R1), Bowtie2 in local mapping mode can accurately map these reads. Mitochondrial mapping R1 reads with a mapping read 2 (R2) pair were extracted with 'samtools view -bh -F 4 f 8', whereas mapping R2 reads with a mapping R1 pair were extracted with 'samtools view -bh -F 8 f 4'. R1 and R2 singleton mapping reads were extracted with 'samtools view -bh -F 12' for diagnostic purposes, but were not used further in the assembly. The R1, R2, and singleton reads in. BAM format were merged, sorted, and converted to FASTQ format with samtools and 'bedtools bamtofastq', respectively. The resultant R1 and R2 FASTQ files containing only the paired mapped reads (995523 pairs, 298 Mbp) were assembled with SPAdes (Nurk et al., 2013) without error correction and with the plasmidSPAdes module (Antipov et al., 2016) enabled (parameters: -t 16 -plasmid -k55,127 -cov-cutoff 1000 -only-assembler). The resulting 'assembly_graph.fastg' file was viewed in Bandage (Wick et al., 2015), revealing a 16,088 bp node with 1119x average coverage that circularized through two possible paths: a 246 bp node with 252x average coverage, or a 245 bp node with 1690x coverage. The lower coverage path was observed to differ only in a 'T' insertion after a 10-nucleotide poly-T stretch when compared to the higher coverage path. Given that increased levels of insertions after polynucleotide stretches are a known systematic error of Illumina sequencing, it was concluded that the lower coverage path represented technical error rather than an authentic genetic variant and was deleted. This produced a single 16,070 bp circular contig. This contig was 'restarted' with seqkit (v0.7.0)  to place the FASTA record break in the AT-rich region, and was submitted to the MITOSv2 mitochondrial genome annotation web server. Small misannotations (e.g. low scoring additional predictions of already annotated mitochondrial genes) were manually inspected and removed. This annotation indicated that all expected features were present on the contig, including subunits of the NAD + dehydrogenase complex (NAD1, NAD2, NAD3, NAD4, NAD4l, NAD5, NAD6), the large and small ribosomal RNAs (rrnL, rrnS), subunits of the cytochrome c oxidase complex (COX1, COX2, COX3), cytochrome b oxidase (COB), ATP synthase (atp6, atp8), and tRNAs. BLASTN of the Ignelater luminosus mitochondrial genome against published complete mitochondrial genomes from beetles indicated 96-89% alignment with 86-73% nucleotide identity, with poor or no sequence level alignment in the A-T rich region. Like other reported elaterid beetle genomes, the I. luminosus mitochondrial genome does not contain the tandem repeat unit (TRU) previously reported in Lampyridae (Bae et al., 2004).   (Emms and Kelly, 2015). See Appendix 4.2.1 for description of clustering method. OGs = Orthogroups, OGS = Official gene set, *=Not completely filtered to single peptide per gene. Figure  produced with InteractiVenn (Heberle et al., 2015). Intermediate scripts and species specific overlaps are available as Figure 2-source data 1. For differential expression testing, Kallisto transcript expression results for P. pyralis (Appendix 1.9.4) and A. lateralis (Appendix 2.7.3) were independently between-sample normalized using Sleuth (v0.30.0) (Pimentel et al., 2017) with default parameters, producing between-sample-normalized transcripts-per-million reads (BSN-TPM). Differential expression (DE) tests for P. pyralis (adult male dissected fatbody vs. adult male dissected lantern -three biological replicates per condition), and for A. lateralis (adult male thorax + abdominal segments 1-5 vs. adult male dissected lantern -three biological replicates per condition), were performed using the Wald test within Sleuth. Genes whose mean BSN-TPM across bioreplicates was above the 90th percentile were annotated as 'highly expressed' (HE). Genes with a Sleuth DE q-value <0.05 were annotated as 'differentially expressed.' (DE). Enzyme encoding (E/NotE) genes were predicted from the InterProScan functional annotations using a custom script (Fallon, 2018d; copy archived at https://github.com/elifesciences-publications/ interproscan_to_enzyme_go) and GOAtools (Tang et al., 2018), with the modification that the enzymatic activity GO term was manually added to select InterPro annotations: IPR029058, IPR036291, and IPR001279. These enzyme lists are available as supporting files associated with the official geneset filesets. Orthogroup membership was determined from the OrthoFinder analysis (Appendix 4.2.1). The enzyme HE/DE/E + NotE gene filtering and overlaps ( Figure 5) were performed using custom scripts. These custom scripts and results of the differential expression testing are available on FigShare (10.6084/m9.figshare.5715151).  . DNA and tRNA methyltransferase gene phylogeny. Levels and patterns of mCG in P. pyralis are corroborated by the presence of de novo and maintenance DNMTs (DNMT3 and DNMT1, respectively). Notably, P. pyralis possesses two copies of DNMT1, and 3 copies of DNMT3, in contrast to a single copy of DNMT1 and DNMT3 in the firefly Aquatica lateralis. The evolutionary history was inferred by using the Maximum Likelihood method with the LG + G (five gamma categories) (Le and Gascuel, 2008). Evolutionary analyses were conducted in MEGA7   (Suzuki et al., 2007), and confidently recovers the presence/absence of DNA methylation in insects (Bewick et al., 2017). In a mixture of loci that are DNA methylated and low to un-methylated, a bimodal distribution of CpG [O/E] values is expected. Conversely, a unimodal distribution is suggestive of a set of loci that are mostly low to un-methylated.
The modality of CpG [O/E] distributions was tested using Gaussian mixture modeling in R (https://www.r-project.org/: mclust v5.4 and mixtools v1.0.4). Two modes were modeled for each CpG [O/E] distribution, and the subsequent means and 95% confidence interval (CI) of the means were compared with overlapping or nonoverlapping CI's signifying unimodality or bimodality, respectively.

CYP303 evolutionary analysis (Figure 6C)
Candidate P450s were identified using BLASTP (e-value: 1 Â 10 À20 ) of a P. pyralis CYP303 family member (PPYR_OGS1.0: PPYR_14345-PA) against the P. pyralis, A. lateralis, and I. luminosus reference set of peptides, and the D. melanogaster (NCBI GCF_000001215.4) and T. castaneum (NCBI GCF_000002335.3) geneset peptides. Resulting hits were merged, aligned with MAFFT E-INS-i (v7.243) (Katoh and Standley, 2013), and a preliminary neighborjoining (NJ) tree was generated using MEGA7 . Genes descending from the common ancestor of the CYP303 and CYP304 genes were selected from this NJ tree, and the peptides within this subset re-aligned with MAFFT using the L-INS-i algorithm. Then the maximum likelihood evolutionary history of these genes was inferred within MEGA7 using the LG + G model (five gamma categories (+G, parameter = 2.4805). Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with the best log likelihood value. The resulting tree was rooted using D. melanogaster Cyp6a17 (NP_652018.1). The tree shown in Figure 6C was truncated in Dendroscope (v3.5.9) (Huson and Scornavacca, 2012) to display only the CYP303 clade. The multiple sequence alignment FASTA files and newick files of the full and truncated tree are available in Figure 6-source data 1.

Luciferase genetics overview
The gene for firefly luciferase was first isolated from the North American firefly P. pyralis (de Wet et al., 1985;Wood et al., 1984;de Wet et al., 1987) and then identified from the Japanese fireflies Luciola cruciata (Masuda et al., 1989) and Aquatica lateralis (Tatsumi et al., 1992). To date, firefly luciferase genes have been isolated from more than 30 lampyrid species in the world. Two different types of luciferase genes, Luc1 and Luc2, have been reported from Photuris pennsylvanica (Ye et al., 1997) (Photurinae), L. cruciata (Oba et al., 2010b) (Luciolinae), A. lateralis (Oba et al., 2013a)  Luciferase genes have also been isolated from members of the other luminous beetles families: Phengodidae, Rhagophthalmidae, and Elateridae (Wood et al., 1989;Viviani et al., 1999a;Viviani et al., 1999b;Ohmiya et al., 2000) with amino acid identities to firefly luciferases at >48% (Oba, 2014). The chemical structures of the substrates for these enzymes are identical to firefly luciferin. These results that the bioluminescence systems of luminous beetles are essentially the same, supports a single origin of the bioluminescence in elateroid beetles. Recent molecular analyses based on the mitochondrial genome sequences strongly support a sister relationship between the three luminous families: Lampyridae, Phengodidae, and Rhagophthalmidae (Timmermans et al., 2010;Timmermans and Vogler, 2012), suggesting the monophyly of Elateroidea and a single origin of the luminescence in the ancestor of these three lineages (Oba, 2014). However, ambiguity in the evolutionary relationships among luminous beetles, including luminous Elaterids, does not yet exclude multiple origins.
Molecular analyses have suggested that the origin of Lampyridae was dated back to late Jurassic (McKenna and Farrell, 2009) or mid-Cretaceous periods (Mckenna et al., 2015). Luciolinae and Lampyrinae was diverged at the basal position of the Lampyridae (Martin et al., 2017) and the fossil of the Luciolinae firefly dated at Cretaceous period was discovered in Burmese amber (Shi et al., 2012;Kazantsev, 2015). Taken  Between fireflies and click-beetles, the structure of the luciferase genes are globally similar, with seven exons, similar intron lengths, and identical splice junction locations (Appendix 4- figure 5). The intron-exon structure of IlumLuc is consistent with the reported intron-exon structure of Pyrophorus plagiophthalamus luciferase (Velez and Feder, 2006).

Alignment and Gene Phylogeny
The 20 merged CDS sequences were multiple-sequenced-aligned with MUSCLE (Edgar, 2004) in 'codon' mode within MEGA7 , using parameters (Gap Open = À0.2.9; Gap Extend = 0; Hydrophobicity Multiplier 1.2, Clustering Method = UPGMB, Min Diag Length (lambda) = 24, Genetic Code = Standard), producing a nucleotide multiplesequence-alignment (MSA). A maximum likelihood gene tree was produced from the nucleotide MSA within MEGA7 using the General Time Reversible model (Nei and Kumar, 2000), with five gamma categories (+G, parameter = 0.8692). The analysis involved 20 nucleotide sequences. Codon positions included were 1st + 2nd + 3rd + Noncoding. There were a total of 1659 positions in the final dataset. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with the superior log likelihood value. The tree with the highest log likelihood (À16392.22) was selected. 1000 bootstrap replicates were performed to evaluate the topology, and the percentage of trees in which the associated taxa clustered together is shown next to the branches in Figure 4B.

Tests of selection: aBSREL
An adaptive branch-site REL test for episodic diversification was performed on the previously mentioned gene-tree and nucleotide MSA using the adaptive branch-site REL test for episodic diversification (aBSREL) method (Smith et al., 2015) within the HyPhy program (v2.3.11) (Pond et al., 2005). The input MSA contained 20 sequences with 553 sites (codons). All 37 branches of the gene phylogeny were formally tested for diversifying selection. The aBSREL analysis found evidence of episodic diversifying selection on 3 out of 37 branches in the phylogeny. Significance was assessed using the Likelihood Ratio Test at a threshold of p 0.01, after the Holm-Bonferroni correction for multiple hypothesis testing. The intermediate files and results of this analysis, including the nucleotide MSA, GTR based gene-tree, and aBSREL produced adaptive rate class model gene tree are available as Figure 4-source data 2.

Tests of selection: MEME
After identification of the selected branch via the aBSREL method, we turned to the MEME method within the HyPhy program (v2.3.11) (Pond et al., 2005), to identify those sites which may have adaptively evolved. We tested the branch leading to EAncLuc, which was previously identified as under selection in the aBSREL analysis. A single partition was recovered with 28 sites under episodic diversifying positive selection at p<=0.1 (Appendix 4-table 5). Input files and full results are available on FigShare (10.6084/m9.figshare.6626651).

Tests of selection: PAML-BEB
To validate our findings from aBSREL and MEME using a different method, we applied Phylogenetic Analysis by Maximum Likelihood (PAML) branch by site analysis to the luciferase sequences. We tested the alternative hypothesis, that there is a class of sites under selection (w >1) on the EAncLuc ancestral branch identified as under selection in the aBSREL analysis, against the null hypotheses, that all classes of sites on all branches are evolving either under constraint (w <1) or neutrality (w = 1). A likelihood ratio test supported the alternative hypothesis, that 13% of sites in luciferase were in a positively selected class (w = 3.25). Subsequent Bayes Empirical Bayes (BEB) estimation identified 31 sites with evidence of selection on these branches, 5 of which were significant. Full results are available on FigShare (10.6084/m9.figshare.6725081).

Tests of selection: Overlap
Nineteen of the overall sites were shared between the MEME analysis, and are shown in Appendix 4-table 5. The frequency of extant amino acids at these sites are shown in Appendix 4-figure 7.
Appendix 4-table 4. Results of PAML branch x sites analysis. Proportion indicates the proportion of sites in each site class (0, 1, 2a, 2b). Site classes 0 and 1 are those in the constrained and neutral classes, respectively. 2a are sites that were constrained on the background branches, but are either neutral (H0) or in the selective class (HA) on the foreground branches. 2b are sites that were neutral on the background branches, but are either neutral (H0) or in the selective class (HA) on the foreground branches.  4.4 Non-enzyme highly and differentially expressed genes of the firefly lantern PPYR_04589, a predicted fatty acid binding protein is almost certainly orthologous to the light organ fatty acid binding protein reported from Luciola cerata (Goh and Li, 2011). This fatty acid binding protein was previously reported to bind strongly to fatty acids, and weakly to luciferin. Notably, PPYR_04589 is the most highly expressed gene in the P. pyralis adult lantern, ahead of firefly luciferase. Three G-coupled protein receptors (GCPRs) with similarity to annotated octopamine/tyramine receptors were also detected to be highly and differentially expressed in the P. pyralis light organ (PPYR_11673-PA, PPYR_11364-PA, PPYR_12266-PA). Octopamine is known to be the key effector neurotransmitter of the adult and larval firefly lantern and this identified GPCR likely serves as the upstream receptor of octopamine activated adenylate cyclase, previously reported as abundant in P. pyralis lanterns (Nathanson et al., 1989). The neurobiology of flash control, including regulation of flash pattern and intensity, is a fascinating area of behavioral research. Our data generate new hypotheses regarding the molecular players in flash control. A particularly interesting highly and differentially expressed gene in both P. pyralis and A. lateralis is the full length 'octopamine binding secreted hemocyanin'(PPYR_14966; AQULA_008529; Appendix 4-table 6) previously identified from P. pyralis light organ extracts via photoaffinity labeling with an octopamine analog and partial N-terminal Edman degradation (Nathanson et al., 1989). This protein is intriguing as hemocyanins are typically thought to be oxygen binding. We speculate that this octopamine binding secreted hemocyanin, previous demonstrated to be abundant, species is also HE, DE, NotE. BSN-TPM = between sample normalized TPM.  1). As this is a genome-wide analysis where bootstrap replicates would be computationally prohibitive, no bootstrap replicates were performed to evaluate the support of the tree topology. PTS1 sequences were predicted from the peptide sequence using the PTS1 predictor server (Neuberger et al., 2017).

Opsin analysis
Opsins are G-protein-coupled receptors that, together with a bound chromophore, form visual pigments that detect light, reviewed here (Briscoe and Chittka, 2001). While opsin genes are known for their expression in photoreceptors and function in vision, they have also been found to be expressed in other tissues, suggesting non-visual functions in some cases. Insects generally use rhabdomeric opsins (r-opsins) for vision, while mammals generally use ciliary opsins (c-opsins) for vision, products of an ancient gene duplication (Briscoe and Chittka, 2001;Porter et al., 2012). Both insects and mammals may retain the alternate opsin type, generally in a non-visual capacity. The ancestral insect is hypothesized to have three visual opsins -one sensitive to long-wavelengths of light (LW), one to blue-wavelengths (B), and one to ultraviolet light (UV). Previously, two opsins, one with sequence similarity to other insect LW opsins and one with similarity to other insect UV opsins, were identified as highly expressed in firefly heads (Sander and Hall, 2015;Martin et al., 2015). A likely nonvisual c-opsin was also detected, although not highly expressed (Sander and Hall, 2015;Martin et al., 2015).
To confirm the previously documented opsin presence and expression patterns, we collected candidate opsin genes via BLASTP searches (e-value threshold: 1 Â 10 À20 ) of the PPYR_OGS1.0, AQULA_OGS1.0 and ILUMI_OGS1.0 reference genesets against UV opsin of P. pyralis (Genbank Accession: ALB48839.1), as well as collected non-firefly opsin sequences via literature searches, followed by maximum likelihood phylogenetic reconstruction (Appendix 4-figure 9A), and expression analyses of the opsins (Appendix 4- figure 9B.B). The amino acid sequences of opsin were multiple aligned using MAFFT and trimmed using trimAL (parameters: -gt 0.5). The amino acid substitution model for ML analysis was estimated using Aminosan (v1.0.2016.11.07) (Tanabe, 2011). In P. pyralis, A. lateralis, and I. luminosus, we detected three r-opsins, including LW, UV, and an r-opsin homologous to Drosophila Rh7 opsin, and one c-opsin. While LW and UV opsins were highly and differentially expressed in heads of both fireflies, c-opsin was lowly expressed, in P. pyralis head tissue only (Appendix 4- figure 9B). In contrast, Rh7 was not expressed in the P. pyralis light organ, but was differentially expressed in the light organ of A. lateralis (Appendix 4- figure 9B). The detection of Rh7 in our genomes is unusual in beetles (Feuda et al., 2016), although emerging genomic resources across the order have detected it in two taxa: Anoplophora glabripennis (McKenna et al., 2016) and Leptinotarsa decemlineata (Schoville et al., 2017). Rh7 has an enigmatic function -a recent study in Drosophila melanogaster showed that Rh7 is expressed in the brain, functions in circadian photoentrainment, and has broad UV-to-visible spectrum sensitivity (Ni et al., 2017;Sakai et al., 2017). Extraocular opsin expression has been detected in other eukaryotes: a photosensory organ is located in the genitalia at the posterior abdominal segments in butterfly (Lepidoptera) (Arikawa and Aoki, 1982). In the bioluminescent Ctenophore Mnemiopsis leidyi, three c-opsins are co-expressed with the luminous photoprotein in the photophores (Schnitzler et al., 2012). In the bobtail squid, Euprymna scolopes, one of the c-opsin isoforms is expressed in the bacterial symbiotic light organ (Tong et al., 2009;Pankey et al., 2014). Thus, it is possible that Rh7 has a photo sensory function in the lantern of A. lateralis, although this putative function is seemingly not conserved in P. pyralis. Future study will confirm and further explore the biological, physiological, and evolutionary significance of Rh7 expression in the light organ across firefly taxa. 4.6 LC-HRAM-MS of lucibufagin content in P. pyralis, A. lateralis, and I. luminosus We assayed the hemolymph of adult P. pyralis and A. lateralis, as well as body extracts from P. pyralis and A. lateralis larvae, and I. luminosus adult male thorax, for lucibufagin content using liquid-chromatography high-resolution accurate-mass mass-spectrometry (LC-HRAM-MS) and MS 2 spectral similarity networking approaches. We chose to analyze extracted hemolymph from both P. pyralis, and A. lateralis for lucibufagin content, as lucibufagins are known to accumulate in the adult hemolymph and hemolymph samples give less complex extracts than tissue extracts. For P. pyralis and A. lateralis larvae, and I. luminosus thorax, tissue extracts were sampled as we do not have a reliable hemolymph extraction protocol for these life stages and species. Specific tissues were chosen for extracts to enable a smaller quantity of tissue to go into the metabolite extraction, and to explore possible difference in compound abundance across tissues, but we expected that defense compounds like lucibufagins would be roughly equally abundant present in all tissues. Adult male P. pyralis and A. lateralis hemolymph was extracted by the following methods: A single live adult P. pyralis male was placed in a 1.5 mL microcentrifuge tube with a 5-mmglass bead underneath the specimen, and centrifuged at maximum speed (~20,000 xg) for 30 s in a benchtop centrifuge. This centrifugation crushed the specimen on top of the bead, and allowed the hemolymph to collect at the bottom of the tube. Approximately 5 mL was obtained. The extracted hemolymph was diluted with 50 mL methanol to precipitate proteins and other macromolecules. For A. lateralis adult hemolymph, three adult male individuals were placed in individual 1.5 mL microcentrifuge tubes with 5-mm-glass beads, and spun at 5000 RPM for 1 min in a benchtop centrifuge. The pooled extracted hemolymph (~5 mL), was diluted with 50 mL MeOH, and air dried. The P. pyralis extracted hemolymph was filtered through a 0.2 mm PFTE filter (Filter Vial, P/No. 15530-100, Thomson Instrument Company), whereas the A. lateralis hemolymph residue was redissolved in 100 mL 50% MeOH, and then filtered through the filter vial.
For extraction of P. pyralis larval partial body, the posterior two abdominal segments were first cut off from a single laboratory reared larvae (Appendix 1.3.2), and the remaining partial body was placed in 180 mL 50% acetonitrile, and macerated with a pipette tip. The extract was sonicated in a water bath sonicator for~10 min, not letting the temperature of the bath go above 50˚C. The extract was then centrifuged (20,000 x g for 10 min), and filtered through a 0.2 mm PFTE filter (Filter Vial, P/No. 15530-100, Thomson Instrument Company).
For extraction of A. lateralis larval whole body, laboratory reared A. lateralis larvae were flash frozen in liquid N 2 , lyophilized, and the whole body (dry weight: 29.1 mg) was placed in 200 mL 50% methanol, and macerated with a pipette tip. The extract was sonicated in a water bath sonicator for 30 min, centrifuged (20,000xg for 10 min), and filtered through a 0.2 mm PFTE filter (Filter Vial, P/No. 15530-100, Thomson Instrument Company).
For extraction of I. luminosus adult thorax, the mesothorax through the two most anterior abdominal segments (ventral lantern containing segment +1 segment) of a lyophilized I. luminosus adult male (Appendix 3.3), was separated from the prothorax plus head and posterior three abdominal segments. This mesothorax + abdomen fragment was then placed in 0.5 mL 50% methanol, and macerated with a pipette tip. The extract was then sonicated in a water bath sonicator for~10 min, not letting the temperature of the bath go above 50˚C, centrifuged (20,000xg for 10 min), and filtered through a 0.2 mm PFTE filter (Filter Vial, P/No. 15530-100, Thomson Instrument Company).
Two different instrument methods were used, a slow~44 min method, and an optimized~28 min method. Chromatographically both methods are identical up to 20 min. P. pyralis hemolymph compounds were separated by the optimized method (28 min), with separation via reversed-phase chromatography on a C18 column using a gradient of Solvent A (0.1% formic acid in H 2 O) and Solvent B (0.1% formic acid in acetonitrile); 5% B for 2 min, 5-40% B until 20 min, 40-95% B until 22 min, 95% B for 4 min, and 5% B for 5 min; flow rate 0.8 mL/min. All other sample extracts were separated by the slow (44 min) reversed-phase chromatography method, using a C18 column with a gradient of Solvent A (0.1% formic acid in H 2 O) and Solvent B (0.1% formic acid in acetonitrile); 5% B for 2 min, 5-80% B until 40 min, 95% B for 4 min, and 5% B for 5 min; flow rate 0.8 mL/min.
The mass spectrometer was configured to perform one MS 1 scan from m/z 120-1250 followed by 1-3 data-dependent MS 2 scans using HCD fragmentation with a stepped collision energy of 10, 15, 25 normalized collision energy (NCE). Positive mode and negative mode MS 1 and MS 2 data were obtained in a single run via polarity switching for the optimized method, or in separate runs for the slow method. Data was collected as profile data. The instrument was always used within 7 days of the last mass accuracy calibration. The ion source parameters were as follows: spray voltage (+) at 3000 V, spray voltage (-) at 2000 V, capillary temperature at 275˚C, sheath gas at 40 arb units, aux gas at 15 arb units, spare gas at one arb unit, max spray current at 100 (mA), probe heater temp at 350˚C, ion source: HESI-II. The raw data in Thermo format was converted to mzML format using ProteoWizard MSConvert (Chambers et al., 2012). Data analysis was performed with Xcalibur (Thermo Scientific) and MZmine2 (v2.30) (Pluskal et al., 2010). Raw LC-MS data is available on MetaboLights (Accession: MTBLS698).
Coding regions, tRNAs, and rRNAs were predicted via the MITOSv2 mitochondrial genome annotation web server (Bernt et al., 2017). Small mis-annotations (e.g. low scoring additional predictions of already annotated mitochondrial genes) were manually inspected and removed. Tandem repetitive regions were manually annotated. The complete A. antennatus genome annotation plus assembly is available on NCBI Genbank (Accession: MG546669).

Taxonomic identification of Phorid mitochondrial genome origin
After the successful metagenomic assembly of the mitochondrial genome of an unknown Phorid fly species from the P. pyralis PacBio library (Appendix 5.2), we sought to characterize the species of origin for this mitochondrial genome. We planned to achieve this by collecting the Phorid flies which emerged from adult P. pyralis, taxonomically identifying them, and performing targeted mitochondrial PCR and sequencing experiments to correlate their mitochondrial genome sequence to our mtDNA assembly. We successfully obtained phorid fly larvae emerging from P. pyralis adult males collected from MMNJ (identical field site to PacBio collection), and Rochester, NY (RCNY), in the summer of 2017. The MMNJ phorid larvae did not successfully pupate, however we obtained five adult specimens from successful pupations of the RCNY larvae. Two adults from this batch were identified as A. antennatus (Malloch), by Brian V. Brown, Entomology Curator of the Natural History Museum of Los Angeles County. DNA was extracted from one of the remaining three specimens and a COI fragment was PCR-amplified and Sanger sequenced. The forward primer was 5'-TTTGATTCTTCGGCCACCCA-3', the reverse primer 5'-AGCATCGGGGTAGTCTGAGT-3'. This COI fragment from had 99% identity (558/563 nt) to the COI gene of our mitochondrial assembly. This sequenced COI fragment has been submitted to GenBank (GenBank Accession: MG517481). We conclude that this is sufficient evidence to denote that our assembled Phorid mitochondrial genome is the mitochondrial genome of A. antennatus. Notably, A. antennatus was previously reported by Lloyd (1973) to be a parasite of several firefly species in genera Photuris, Photinus, and Pyractomena, from collection sites ranging from Florida to New York. To our knowledge, this is the first report of a mitochondrial genome which was first assembled via an untargeted metagenomic approach and then later correlated to its species of origin.

Photinus pyralis orthomyxo-like viruses
We identified the first two viruses associated to P. pyralis and the Lampyridae family. The proposed Photinus pyralis orthomyxo-like virus 1 and 2 (PpyrOMLV1 and 2) present a multipartite genome conformed by five RNA segments encoding a putative nucleoprotein (NP), hemagglutinin-like glycoprotein (HA) and a heterotrimeric viral RNA polymerase (PB1, PB2 and PA). The viral genomes for Photinus pyralis orthomyxo-like virus 1 and 2 are available on NCBI Genbank with accessions MG972985-MG972994. Expression analyses on 24 RNA libraries of diverse individuals/developmental stages/tissues and geographic origins of P. pyralis indicate a dynamic presence, widespread prevalence, a pervasive tissue tropism, a low isolate variability, and a persistent life cycle through transovarial transmission of PpyrOMLV1 and 2. Genomic and phylogenetic studies suggest that the detected viruses correspond to a new lineage within the Orthomyxoviridae family (ssRNA(-)) (Appendix 5-figure 2A-I). The concomitant occurrence in the P. pyralis genome of species-specific signatures of Endogenous viral-like elements (EVEs) associated to retrotransposons linked to the identified Orthomyxoviruses, suggest a past evolutionary history of host-virus interaction (Appendix 5.5, Appendix 5- figure 2J). This tentative interface is correlated to low viral RNA levels, persistence and no apparent phenotypes associated with infection. We suggest that the identified viruses are potential endophytes of high prevalence as a result of potential evolutionary modulation of viral levels associated to EVEs. Photinus pyralis orthomyxo-like virus 1 and 2 (PpyrOMLV1 and PpyrOMLV2) share their genomic architecture and evolutionary clustering (Appendix 5- figure 2A-H, Appendix 5- figure 3).  figure 2E). Despite sharing genome architecture and structural and functional domains of their predicted proteins, PpyrOMLV1 and PpyrOMLV2 pairwise identity of ortholog gene products range between 21.4% (HA) to 49.8% (PB1), suggesting although a common evolutionary history, a strong divergence indicating separated species, borderline to be considered even members of different virus genera (Appendix 5- figure 3). The conserved 3' sequence termini of the viral genomic RNAs are (vgRNA ssRNA(-) 3'-end) 5'-GUUCUUACU-3' for PpyrOMLV1, and and 5'-(G/A)U(U/G)(G/U/C)(A/C/U)UACU-3'. for PpyrOMLV2. The 5' termini of the vgRNAs are partially complementary to the 3' termini, supporting a panhandle structure and a hook like structure of the 5' end by a terminal short stem loop. PpyrOMLV1 and PpyrOMLV2 genome segments present an overall high identity in their respective RNA segments ends (Appendix 5- figure 2F). These primary and secondary sequence cues are associated to polymerase binding and promotion of both replication and transcription. In influenza viruses, and probably every OMV, the first 10 nucleotides of the 3 0 end form a stem-loop or 'hook' with four base-pairs (two canonical base-pairs flanked by an A-A base-pair). This compact RNA structure conforms the promoter, which activates polymerase initiation of RNA synthesis . The presence of eventual orthologs of OMV additional genome segments and proteins, such as Neuraminidase (NA), Matrix (M) and Non-structural proteins (NS1, NS2) was assessed retrieving no results by TBLASTN relaxed searches, nor with in silico approaches involving co-expression, expression levels, or conserved terminis. Given that the presence of those additional segments varies among diverse OMV genera, and that 35 related tentative new virus species identified in TSA did not present any additional segments, we believe that these lineages of viruses are conformed by five genome segments. Further experiments based on specific virus particle purification and target sequencing could corroborate our results. Based on sequence homology to best BLASTP hits, amino acid sequence alignments, predicted proteins and domains, and phylogenetic comparisons to reported species we assigned PpyrOMLV1 and PpyrOMLV2 to the OMV virus family. These are the first viruses that have been associated with the Lampyridae beetle family, which includes over 2000 species. The OMV virus members share diverse structural, functional and biological characters that define and restrict the family. OMV virions are 80-120 nm in diameter, of spherical or pleomorphic morphology. The virion envelope is derived from the host cell membrane,  14.6 Kbp (King et al., 2011). As described previously, every OMV RNA segment possess conserved and partially complementary 5 0 -and 3 0 -end sequences with promoter activity (Hsu et al., 1987). OMV structural proteins are tentatively common to all genera involving the three polypeptides subunits that form the viral RdRP (PA, PB1, PB2) (Pflug et al., 2017); a nucleoprotein (NP), which binds with each genome ssRNA segment to form RNPs; and the hemagglutinin protein (HA, HE or GP), which is a type I membrane integral glycoprotein involved in virus attachment, envelope fusion and neutralization. In addition, a nonglycosylated matrix protein (M) is present in most species. There are some species-specific divergence in some structural OMVs proteins. For instance, HA of FLUAV is acylated at the membrane-spanning region and has widespread N-linked glycans (Eisfeld et al., 2015). The HA protein of FLUCV, besides its hemagglutinating and envelope fusion function, has an esterase activity that induces host receptor enzymatic destruction (King et al., 2011). In contrast, the HA of THOV is divergent to influenzavirus HA proteins, and presents high sequence similarity to a baculovirus surface glycoprotein (Leahy et al., 1997). The HA protein has been described to have an important role in determining OMV host specificity. For instance, human infecting Influenza viruses selectively bind to glycolipids that contain terminal sialyl-galactosyl residues with a 2-6 linkage, in contrast, avian influenza viruses bind to sialylgalactosyl residues with a 2-3 linkage (King et al., 2011). Furthermore, FLUAV and FLUBV share a neuraminidase protein (NA), which is an integral, type II envelope glycoprotein containing sialidase activity. Some OMVs possess additional small integral membrane proteins (M2, NB, BM2, or CM2) that may be glycosylated and have diverse functions. As an illustration, M2 and BM2 function during un-coating and fusion by equilibrating the intralumenal pH of the trans-Golgi apparatus and the cytoplasm. In addition, some viruses encode two nonstructural proteins (NS1, NS2) (King et al., 2011). OMV share replication properties, which have been studied mostly in Influenza viruses. It is important to note that gene reassortment has been described to occur during mixed OMV infections, involving viruses of the same genus, but not between viruses of different genera (Kimble, 2013). This is used also as a criteria for OMV genus demarcation. Influenza virus replication and transcription occurs in the cell nucleus and comprises the production of the three types of RNA species (i) genomic RNA (vRNA) which are found in virions; (ii) cRNA molecules which are complementary RNA in sequence and identical in length to vRNA; and also (iii) virus mRNA molecules which are 5' capped by cap snatching of host RNAs and 3' polyadenylated by polymerase stuttering on U rich stretches. These remarkable dynamic multifunction characters of OMV polymerases are associated with its complex tertiary structure, of this modular heterotrimeric replicase (Te Velthuis and Fodor, 2016). We explored in detail the putative polymerase subunits of the identified firefly viruses. The PB1 subunit catalyzes RNA synthesis in its internal active site opening, which is formed by the highly conserved polymerase motifs I-III. Motifs I and III (Appendix 5- figure 2H) present three conserved aspartates (PpyrOMLV1: Asp 346, Asp 491 and Asp 492; PpyrOMLV2: Asp 348, Asp 495 and Asp 496) which coordinate and promote nucleophilic attack of the terminal 3' OH from the growing transcript on the alpha-phosphate of the inbound NTP (Pflug et al., 2017). Besides presenting, with high confidence, the putative functional domains associated with their potential replicase/transcriptase function, we assessed whether the potential spatial and functional architecture was conserved at least in part in FOML viruses. In this direction we employed the SWISS-MODEL automated protein structure homology-modelling server to generate a 3D structure of PpyrOMLV1 heterotrimeric polymerase. The SWISS server selected as best-fit template the trimeric structure of Influenza A virus polymerase, generating a structure for each polymerase subunit of PpyrOMLV1. The generated structure shared structural cues related to its multiple role of RNA nucleotide binding, endonuclease, cap binding, and nucleotidyl transferase (Appendix 5- figure 2G-H).
The engendered subunit structures suggest a probable conservation of PpyrOMLV1 POL, that could allow the predicted functional enzymatic activity of this multiple gene product. The overall polymerase rendered structure presents a typical U shape with two upper protrusions corresponding to the PA endonuclease and the PB2 cap-binding domain. The PB1 subunit appears to plug into the interior of the U and has the distinctive fold of related viral RNA polymerases with fingers, palm and thumb adjacent to a tentative central active site opening where RNA synthesis may occur Hengrung et al., 2015). OMV Pol activity is central in the virus cycle of OMVs, which have been extensively studied. The life cycle of OMVs starts with virus entry involving the HA by receptor-mediated endocytosis. For Influenza, sialic acid bound to glycoproteins or glycolipids function as receptor determinants of endocytosis. Fusion between viral and cell membranes occurs in endosomes. The infectivity and fusion of influenza is associated to the post-translational cleavage of the virion HA. Cleavability depends on the number of basic amino acids at the target cleavage site (King et al., 2011). In thogotoviruses, no requirement for HA glycoprotein cleavage have been demonstrated (Leahy et al., 1997). Integral membrane proteins migrate through the Golgi apparatus to localized regions of the plasma membrane. New virions form by budding, incorporating matrix proteins and viral RNPs. Viral RNPs are transported to the cell nucleus where the virion polymerase complex synthesizes mRNA species (Hara et al., 2017). Another tentative function of the NP could be associated to the potential interference of the host immune response in the nucleus mediated by capsid proteins of some RNA virus, which could inhibit host transcription and thus liberate and direct it to viral RNA synthesis (Wulan et al., 2015). mRNA synthesis is primed by capped RNA fragments 10-13 nt in length that are generated by cap snatching from host nuclear RNAs which are sequestered after cap recognition by PB2 and incorporated to vRNA by PB1 and PA proteins which present viral endonuclease activity (Sikora et al., 2017). In contrast, thogotoviruses have capped viral mRNA without host-derived sequences at the 5 0 end. Virus mRNAs are polyadenylated at the 3 0 termini through iterative copying by the viral polymerase stuttering on a poly U track in the vRNA template. Some OMV mRNAs are spliced generating alternative gene products with defined functions. Protein synthesis of influenza viruses occurs in the cytoplasm. Partially complementary vRNA molecules act as templates for new viral RNA synthesis and are neither capped nor polyadenylated. These RNAs exist as RNPs in infected cells. Given the diverse hosts of OMV, biological properties of virus infection diverge between species. Influenzaviruses A infect humans and cause respiratory disease, and they have been found to infect a variety of bird species and some mammalian species. Interspecies transmission, although rare, is well documented. Influenza B virus infect humans and cause epidemics, and have been rarely found in seals. Influenzaviruses C cause limited outbreaks in humans and have been occasionally found on dogs. Influenza spreads globaly in a yearly outbreak, resulting in about three to five million cases of severe illness and about 250,000 to 500,000 human deaths (Thompson et al., 2009). Influenzavirus D has been recently reported and accepted and infects cows and swine (Hause et al., 2013). Natural transmission of influenzaviruses is by aerosol (human and non-aquatic hosts) or is water-borne (avians). In contrast, Thogoto and Dhori viruses which also infect humans, are transmitted by, and able to replicate in ticks. Thogoto virus was identified in Rhipicephalus sp. ticks collected from cattle in the Thogoto forest in Kenya, and Dhori virus was first isolated in India from Hyalomma dromedarri, a species of camel ticks (Anderson and Casals, 1973;Haig et al., 1965). Dhori virus infection in humans causes a febrile illness and encephalitis. Serological evidence suggests that cattle, camel, goats, and ducks might be also susceptible to this virus. Experimental hamster infection with THOV may be lethal. Unlike influenzaviruses, these viruses do not cause respiratory disease. The transmission of fish infecting isaviruses (ISAV) is via water, and virus infection induces the agglutination of erythrocytes of many fish species, but not avian or mammalian erythrocytes (Mjaaland et al., 1997). Quaranfil and Johnston Atoll are transmitted by ticks and infect avian species (Presti et al., 2009). We have limited biological data of the firefly detected viruses. Nevertheless, a significant consistency in the genomic landscape and predicted gene products of the detected viruses in comparison with accepted OMV species sufficed to suggest for PpyrOMLV1 and PpyrOMLV2 a tentative taxonomic assignment within the OMV family. Besides relying on the OMV structural and functional signatures determined by virus genome annotation, we explored the evolutionary clustering of the detected viruses by phylogenetic insights. We generated MAFFT alignments and phylogenetic trees of the predicted viral polymerase of firefly viruses and the corresponding replicases of all 493 proposed and accepted species of ssRNA(-) virus. The generated trees consistently clustered the diverse sequences to their corresponding taxonomical niche, at the level of genera. Interestingly, PpyrOMLV1 and PpyrOMLV2 replicases were placed unequivocally within the OMV family (Appendix 5- figure 2B). When the genetic distances of firefly viruses proteins and ICTV accepted OMV species were computed, a strong similarity was evident (Appendix 5- figure 2B-D). Overall similarity levels of PpyrOMLV polymerase subunits ranged between 11.03% to as high as 37.30% among recognized species, while for the more divergent accepted OMV (ISAV -Isavirus genus) these levels ranged only from 8.54% to 20.74%, illustrating that PpyrOMLV are within the OMV by genetic standards. Phylogenetic trees based on aa alignments of structural gene products of recognized species and PpyrOMLV supported this assignment, placing ISAV and issavirus as the most distant species and genus within the family, and clustering PpyrOMLV1 and PpyrOMLV2 in a distinctive lineage within OMV, more closely related to the Quaranjavirus and Thogotovirus genera than the Influenza A-D or Isavirus genera (Appendix 5- figure 3).
Furthermore, it appears that virus genomic sequence data, while it has been paramount to separate species, in the case of genera, there are some contrasting data that should be taken into consideration. For instance, DHOV and THOV are both members of the Thogotovirus genus, sharing a 61.9% and a 34.9% identity at PB1 and PB2, respectively. However, FLUCV and FLUDV are assigned members of two different genus, Influenzavirus C and Influenzavirus D, while sharing a higher 72.2% and a 52.2% pairwise identity at PB1 and PB2, respectively (Appendix 5- figure 3). In addition, FLUAV and FLUBV, assigned members of two different genus, Influenzavirus A and Influenzavirus D present a comparable identity to that of DHOV and THOV thogotoviruses, sharing a 61% and a 37.9% identity at PB1 and PB2, respectively. It is worth noting that similarity thresholds and phylogenetic clustering based in genomic data have been used differently to demarcate OMV genera, hence there is a need to eventually reevaluate a series of consensus values, which in addition to biological data, would be useful to redefine the OMV family. Perhaps, these criteria discrepancies are more related to a historical evolution of the OMV taxonomy than to pure biological or genetic standards. In contrast to FLUDV, JOV and QUAV, the other virus members of OMV have been described, proposed and assigned at least 34 years ago.
The potential prevalence, tissue/organ tropism, geographic dispersion and lifestyle of PpyrOMLV1 and 2 were assessed by the generation and analyses of 29 specific RNA-Seq libraries of P. pyralis (Appendix 1-table 1). As RNA was isolated from independent P. pyralis individuals of diverse origin, wild caught or lab reared, the fact that we found at least one of the PpyrOMLV present in 82% of the libraries reflects a widespread presence and potentially a high prevalence of these viruses in P. pyralis (Appendix 5-figure 2J, Appendix 5-table 3, S5.4.6). Wild caught individuals were collected in period spanning six years, and locations separated as much as 900 miles (New Jersey -Georgia, USA). Interestingly PpyrOMLV1 and 2 were found in individuals of both location, and the corresponding assembled isolate virus sequences presented negligible differences, with an inter-individual variability equivalent to that of isolates (0.012%). A similar result was observed for virus sequences identified in RNA libraries generated from samples collected in different years. We were not able to identified fixed mutations associated to geographical or chronological cues. Further experiments should explore the mutational landscape of PpyrOMLV1 and 2, which appears to be significantly lower than of Influenzaviruses, specifically Influenza A virus, which are characterized by high mutational rate (ca. one mutation per genome replication) associated to the absence of RNA proofreading enzymes (Pauly et al., 2017). In addition we evaluated the presence of PpyrOMLV1 and 2 on diverse tissues and organs of P. pyralis. Overall virus RNA levels were generally low, with an average of 9.47 FPKM on positive samples. However, PpyrOMLV1 levels appear to be consistently higher than PpyrOMLV2, with an average of 20.50 FPKM for PpyrOMLV1 versus 4.22 FPKM for PpyrOMLV2 on positive samples. When the expression levels are scrutinized by genome segment, HA and NP encoding segments appear to be, for both viruses, at higher levels, which would be in agreement with other OMV such as Influenzaviruses, in which HA and NP proteins are the most expressed proteins, and thus viral mRNAs are consistently more expressed (King et al., 2011). Nevertheless, these preliminary findings related to expression levels should be taken cautiously, given the small sample size. Perhaps, the more remarkable allusion derived from the analyses of virus presence is related to tissue and organ deduced virus tropism. Strikingly, we found virus transcripts in samples exclusively obtained from light organs, complete heads, male or female thorax, female spermatheca, female spermatophore digesting glands and bursa, abdominal fat bodies, male reproductive spiral gland, and other male reproductive accessory glands (Appendix 5-table 3, S5.4.6), indicating a widespread tissue/organ tropism of PpyrOMLV1 and 2. This tentatively pervasive tropism of PpyrOMLV1 and 2 emerges as a differentiation character of these viruses and accepted OMV. For instance, influenza viruses present a epithelial cell-specific tropism, restricted typically to the nose, throat, and lungs of mammals, and intestines of birds. Tropism has consequences on host restriction. Human influenza viruses mainly infect ciliated cells, because attachment of all influenza A virus strains to cells requires sialic acids. Differential expression of sialic acid residues in diverse tissues may prevent cross-species or zoonotic transmission events of avian influenza strains to man (Zeng et al., 2013). Tropism has also influence in disease associated effects of OMV. Some influenza A virus strains are more present in tracheal and bronchial tissue which is associated with the primary lesion of tracheobronchitis observed in typical epidemic influenza. Other influenza A virus strains are more prevalent in type II pneumocytes and alveolar macrophages in the lower respiratory tract, which is correlated to diffuse alveolar damage with avian influenza (Mansfield, 2007). The presence of PpyrOMLV1 and 2 virus RNA in reproductive glands raises some potential of the involvement of sex in terms of prospective horizontal transmission. Given that most libraries corresponded to 3-6 pooled individuals samples of specific organs/tissue, direct comparisons of virus RNA levels were not always possible. However, this valuable data gives important insights into the widespread potential presence of the viruses in every analyzed organ/tissue. Importantly, RNA levels of the putative virus segments shared co-expression levels and a systematic pattern of presence/absence, supporting the suggested multipartite nature of the viruses. We observed the presence of virus RNA of both PpyrOMLV1 and 2 in eight of the RNA-Seq libraries, thus mixed infections appear to be common. Interestingly, we did not observe in any of the 24 virus positive samples evidence of reassortment. Reassortment is a common event in OMV, a process by which influenza viruses swap gene segments. Genetic exchange is possible due to the segmented nature of the OMV viral genome and may occur during mixed infections. Reassortment generates viral diversity and has been associated to host gain of Influenzavirus (Steel and Lowen, 2014). Reassorted Influenzavirus have been reported to occasionally cross the species barrier, into birds and some mammalian species like swine and eventually humans. These infections are usually dead ends, but sporadically, a stable lineage becomes established and may spread in an animal population (Kimble, 2013). Besides its evolutionary role, reassortment has been used as a criterion for species/genus demarcation, thus the lack of observed gene swap in our data supports the phylogenetic and sequence similarity insights that indicates species separation of PpyrOMLV1 and 2.
In light of the presence of virus RNA in reproductive glands, we further explored the potential life style of PpyrOMLV1 and 2 related to eventual vertical transmission. Vertical transmission is extremely exceptional for OMV, and has only been conclusively described for the Infectious salmon anemia virus (Isavirus) (Marshall et al., 2014). In this direction, we were able to generate a strand-specific RNA-Seq library of one P. pyralis adult female PpyrOMLV1 virus positive (parent), another library from seven eggs of this female at~13 days post fertilization, and lastly an RNA-Seq library of four 1 st instar larvae (offspring). When we analyzed the resulting RNA reads, we found as expected virus RNA transcripts of every genome segment of PpyrOMLV1 in the adult female library. Remarkably, we also found PpyrOMLV1 sequence reads of every genome segment of PpyrOMLV1 in both the eggs and larvae samples. Moreover, virus RNA levels fluctuated among the different developmental stages of the samples. The average RNA levels of the adult female were 41.10 FPKM, in contrast, the fertilized eggs sample had higher levels of virus related RNA, averaging at 61.61 FPKM and peaking at the genome segment encoding NP (104.49 FPKM). Interestingly, virus RNA levels appear to drop in first instar larvae, in the sequenced library average virus RNA levels were of 10.42 FPKM. Future experiments should focus on PpyrOMLV1 and 2 virus titers at extended developmental stages to complement these preliminary results. However, it is interesting to note that the tissue specific library corresponding to female spermatheca, where male sperm are stored prior to fertilization, presented relatively high levels of both PpyrOMLV1 and 2 virus RNAs, suggesting that perhaps during early reproductive process and during egg development virus RNAs tend to raise. This tentatively differential and variable virus RNA titers observed during development could be associated to an unknown mechanism of modulation of latent antiviral response that could be repressed in specific life cycle stages. Further studies may validate these results and unravel a mechanistic explanation of this phenomenon. Nevertheless, besides the preliminary developmental data, the consistent presence of PpyrOMLV1 in lab-reared, isolated offspring of an infected P. pyralis female is robust evidence demonstrating mother-to-offspring vertical transmission for this newly identified OMV.
One of many questions that remains elusive here is whether PpyrOMLV1 and 2 are associated with any potential alteration of phenotype of the infected host. We failed to unveil any specific effect of the presence of PpyrOMLV1 and 2 on fireflies. It is worth noting that subtle alterations or symptoms would be difficult to pinpoint in these insects. Future studies should enquire whether PpyrOMLV1 and 2 may have any influence in biological attributes of fireflies such as fecundity, life span or life cycle. Nevertheless, we observed in our data some hints that could be indicative of a chronic state status, cryptic or latent infection of firefly individuals: (i) virus positive individuals presented in general relatively low virus RNA levels. (ii) virus RNA was found in every assessed tissue/organ. (iii) vertical transmission of the identified viruses. The first hint is hardly conclusive, it is difficult to define what a relatively low RNA level is, and high virus RNA loads are not directly associated with disease on reported OMV. The correlation of high prevalence, prolonged host infection, and vertical transmission observed in several new mosquito viruses has resulted in their classification as 'commensal' microbes. A shared evolutionary history of viruses and host, based in strategies of immune evasion of the viruses and counter antiviral strategies of the host could occasionally result in a modulation of viral loads and a chronic but latent state of virus infection (Hall et al., 2016).  identify FEVEs that covered as high as ca. 60% of the corresponding gene product, and in addition, although at specific short protein regions of the putative related FOLMV, similarity values were as high as 89% pairwise identity. In addition, most of the detected FEVEs were flanked by Transposable Elements (TE) (Appendix 5-figure 2J) suggesting that integration followed ectopic recombination between viral RNA and transposons. We found several conserved domains associated to reverse transcriptases and integrases adjacent to the corresponding FEVEs, which supports the hypothesis that these virus-like elements could be reminiscent of an OMV-like ancestral virus that could have been integrated into the genome by occasional sequestering of viral RNAs by the TE machinery. The finding of EVEs in the P. pyralis genome is not trivial, OMV EVEs are extremely rare. There has been only one report of OMV like sequences integrated into animal host genomes, which is the case of Ixodes scapularis, the putative vector of Quaranfil virus and Johnston Atoll virus corresponding to genus Quaranjavirus (Katzourakis and Gifford, 2010). The fact that besides FEVEs, the only other OMV EVE corresponded to an Arthropod genome, given the ample studies of bird and mammal genomes, is suggestive that perhaps OMV EVEs are restricted to Arthropod hosts. Sequence similarity of FEVEs and firefly viruses suggest that these viral 'molecular fossils' could have been tightly associated to PpyrOLMV1 and 2 ancestors. Moreover, we found potential NP and PB1 EVEs in our genome of light emitting click beetle Ignelater luminosus (Elateridae), an evolutionary distant coleoptera. Sequence similarity levels of the corresponding EVEs averaging 52%, could not be related with evolutionary distances of the hosts. We were not able to generate conclusive phylogenetic insights of the detected EVEs, given their partial, truncated and altered nature of the virus like sequences. In specific cases such as PB1-like EVEs there appears to be a trend suggesting an indirect relation between sequence identity and evolutionary status of the firefly host, but this preceding findings should be taken cautiously until more gathered data is available. The widespread presence of DNA sequences significantly similar to OMV in the explored firefly and related genomes are an interesting and intriguing result. At this stage is prudently not to venture to suggest more likely one of the two plausible explanations of the presence of these sequences in related beetles genomes: (i) Ancestral OMV like virus sequences were retrotranscribed and incorporated to an ancient beetle, followed by speciation and eventual stabilization or lost of EVEs in diverse species. (ii) Recent and recursive integration of OMV like virus sequences in fireflies and horizontal transmission between hosts. These propositions are not mutually exclusive, and may be indistinctly applied to specific cases. Future studies should enquire in this genome dark matter to better understand this interesting phenomenon. When more data is available EVE sequences may be combined with phylogenetic data of host species to expose eventual patterns of inter-class virus transmission. Either way, more studies are needed to explore these proposals, Katzourakis and Gifford (Katzourakis and Gifford, 2010) suggested that EVEs could reveal novel virus diversity and indicate the likely host range of virus clades. After identification and confirmation that firefly related EVEs are present in the host DNA genome, an obvious question follows: Are these EVEs just signatures of an evolutionary vestige of stochastic past infections; or could they be associated with an intrinsic function? It has been suggested that intensity and prevalence of infection may be a determinant of EVEs integration, and that exposure to environmental viruses may not (Olson and Bonizzoni, 2017). Previous reports have suggested that EVEs may firstly function as restriction factors in their hosts by conferring resistance to infection by exogenous viruses, and the eventual counter-adaptation of virus populations of EVE positive hosts, could reduce the EVE restriction mechanism to a non-functional status (Aiewsakun and Katzourakis, 2015). Recently, in mosquitoes, a new mechanism of antiviral immunity against RNA viruses has been proposed, relying in the production and expression of EVEs DNA (Goic et al., 2016). Alternatively, eventual EVE expression could lend to the production viral like truncated proteins that may compete in trans with virus proteins from infecting viruses and limit viral replication, transcription or virion assembly (Aaskov et al., 2006). In addition, integration and eventual modulation in the host genome may be associated with an interaction between viral RNA and the mosquito RNAi machinery (Goic et al., 2013). The piRNA pathway mediates through small RNAs and Piwi-Argonaut proteins the repression of TE-derived nucleic acids based on sequence complementarity, and has also been associated to regulation of arbovirus viral-related RNA, suggesting a functional connection among resistance mechanisms against RNA viruses and TEs (Palatini et al., 2017;Miesen et al., 2016). Furthermore, arbovirus EVEs have been linked to the production of viral-derived piRNAs and virus-specific siRNA, inducing host cell immunity without limiting viral replication, supporting persistent and chronic infection (Goic et al., 2016). Perhaps, an EVEdependent mechanism of modulation of virus infection could have some level of reminiscence to the paradigmatic CRISPR/Cas system which mediates bacteriophage resistance in prokaryotic hosts.
In sum, genomic studies are a great resource for the understanding of virus and host evolution. Here, we glimpsed an unexpected hidden evolutionary tale of firefly viruses and related FEVEs. Animal genomes appear to reflect as a book, with many dispersed sentences, an antique history of ancestral interaction with microbes, and EVEs functioning as virus related bookmarks. The exponential growth of genomic data would help to further understand this complex and intriguing interface, in order to advance not only in the apprehension of the phylogenomic insights of the host, but also explore a multifaceted and dynamic virome that has accompanied and even might have shifted the evolution of the host.