Firefly genomes illuminate parallel origins of bioluminescence in beetles

1Fireflies (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, while non-adult bioluminescence is likely a warning signal of their unpalatable chemical defenses, such as the cardiotoxic lucibufagins of Photinus fireflies. The biochemical understanding of firefly luminescence: an ATP, Mg2+, and O2-dependent luciferase-mediated oxidation of the substrate luciferin, along with the cloning of the luciferase gene, led to the widespread use of luciferase as a reporter with unique applications in biomedical research and industry. 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). These four closely related families (superfamily Elateroidea) have homologous luciferases and structurally identical luciferins, implying a single origin of beetle bioluminescence. However, as Darwin recognized in his ‘Difficulties on Theory’, the 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.

2To 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 and the cloning of luciferase, while A. lateralis, a species with specialized aquatic larvae, is one of the few fireflies that can be reliably cultured in the laboratory. These two fireflies represent the two major firefly subfamilies, Lampyrinae and Luciolinae, which diverged from a common ancestor over 100 Mya (Figure 1B). 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 (1849-1929) to first establish the enzymatic basis of bioluminescence in the late 1800s. Comparative analyses of the genomes of these three species allowed us to reconstruct the origin(s) and evolution of beetle bioluminescence.

3Photinus 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 short-insert, mate-pair, Hi-C, and long-read Pacific Biosciences (PacBio) approaches (Appendix 4—table 1). These datasets were combined in a MaSuRCA hybrid genome assembly (Appendix 1.5). The Aquatica lateralis genome was derived from an ALL-PATHs 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 linked-read draft genome (Appendix 3.5), which was further manually scaffolded using low-coverage long-read Oxford Nanopore MinION sequencing (Appendix 3.5.4).

4The 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), 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; Appendix 4—table 2) 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.

5To 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). 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). We hypothesized that the P. pyralis orthologs of known bioluminescence genes, including the canonical luciferase Luc1 and the specialized luciferin sulfotransferase LST, would be located on the same linkage group to facilitate chromosomal looping and enhancer assisted co-expression within the light organ. We, however, found these genes on separate linkage groups (Figure 2A).

6In addition to nuclear genome assembly and coding gene annotation, we also assembled the complete mitochondrial genomes (mtDNA) of P. pyralis (Figure 2C; Appendix 1.8) and I. luminosus (Appendix 3.10), while the mtDNA sequence of A. lateralis was recently published. These mtDNA assemblies show high conservation of gene content and synteny, with the exception of the variable ~1 Kbp tandem repeat unit (TRU) found in the firefly mtDNAs.

7As repetitive elements are common participants and drivers of genome evolution, we next sought to characterize the repeat content of our genome assemblies. Overall, 42.6%, 19.8%, and 34.1% of the P. pyralis, A. lateralis, and I. luminosus assemblies were found to be repetitive, respectively (Appendix 1.11; 2.9; 3.9). Of these repeats 66.7%, 39.4%, and 55% could not be classified as any known repetitive sequence, respectively. Helitrons, DNA transposons that transpose through rolling circle replication, are among the most abundant individual repeat elements in the P. pyralis assembly. Via in situ hybridization, we identified that P. pyralis chromosomes have canonical telomeres with telomeric repeats (TTAGG) (Figure 2B; Appendix 1.13).

8DNA methylation is common in eukaryotes, but varies in degree across insects, especially within Coleoptera. Furthermore, the functions of DNA methylation across insects remain obscure. To examine firefly cytosine methylation, we characterized the methylation status of P. pyralis DNA with whole genome bisulfite sequencing (WGBS). Methylation at CpGs (mCG) was unambiguously detected at ~20% within the genic regions of P. pyralis and its methylation levels were at least twice those reported from other holometabolous insects (Figure 2D; Appendix 1.12). Molecular evolution analyses of the DNA methyltransferases (DNMTs) show that direct orthologs of both DNMT1 and DNMT3 were conserved in P. pyralis, A. lateralis, and I. luminosus (Appendix 4—figure 2; Appendix 4.2.3), implying that our three study species, and inferentially likely most firefly lineages, possess mCG. Corroborating this claim, CpG[O/E] analysis of methylation indicated our three study species had DNA methylation (Appendix 4—figure 3).

9Two luciferase paralogs have been previously described in fireflies. P. pyralis Luc1 was the first firefly luciferase cloned, and its direct orthologs have been widely identified from other fireflies. The luciferase paralog Luc2 was previously known only from a handful of Asian taxa, including A. lateralis. 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.

10Firefly luciferase is hypothesized to be derived from an ancestral peroxisomal fatty acyl-CoA synthetase (PACS) (Figure 3A). We found that, in both firefly species, Luc1 is genomically clustered with its closely related homologs, including PACSs and non-peroxisomal 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 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).

11Luc2 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).

12To 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, 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.

13We 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 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 dN/dS > 1 (ω or max dN/dS = 3.98) as compared to the evolution of its paralogs using the aBSREL branch-site selection test (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. 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.

14Beyond 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. 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.

15Moreover, 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 and recycling. 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 4—figure 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.

16We 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 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.

17Firefly bioluminescence is postulated to have first evolved as an aposematic warning of larval chemical defenses. Lucibufagins are abundant unpalatable defense steroids described from certain North American firefly species, most notably in the genera Photinus, Lucidota, and Ellychnia, 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-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, but to date there are no reports of lucibufagins in the Luciolinae.

18The 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). We hypothesized that cytochrome P450s, an enzyme family widely involved in metabolic diversification of organic substrates, 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, 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. 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. 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.

19Given the increasingly recognized contributions of symbionts to host metabolism, 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 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.

20Within our unfiltered A. lateralis genomic assembly (Alat1.2), we also found 43 scaffolds (2.3 Mbp; GC:29.8%, ~64x coverage), whose taxonomic annotation corresponded to the Tenericutes (Appendix 2.5.2), suggesting that A. lateralis may also harbor a mollicute symbiont. Alat1.2 also contains 2119 scaffolds (13.0 Mbp, GC:63.7%, ~25x coverage) annotated as of Proteobacterial origin. Limited Proteobacterial symbionts were detected in the I. luminosus assembly (0.4 Mbp; GC:30–65% ~10x coverage) (Appendix 3.5.2), suggesting no stable symbiont is present in adult I. luminosus. 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, and the second report of Orthomyxoviridae derived EVEs. 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.

21Here, 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. 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. Toward this goal, the recently published Pyrocoelia pectoralis Lampyrinae firefly genome is an important advance which will contribute to future phylogenetic and evolutionary studies.

22The 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. 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, and that a compound resulting from the spontaneous coupling of benzoquinone and cysteine acts as a luciferin biosynthetic intermediate in A. lateralis. Benzoquinone is known to be a defense compound of distantly related beetles and other arthropods (e.g. millipedes). 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.

23The 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 startle-response mediated increase in hemolymph oxygenation through spiracle opening and escape locomotion caused a concomitant increase in luminescence. Alternatively, an early role for firefly luminescence in mate attraction has not been ruled out. 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 MS2 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.

24Overall, 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.

25Detailed 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.

26Genomic assemblies (Ppyr1.3, Alat1.3, and Ilumi1.2), associated official geneset data, a SequenceServer BLAST server, and a JBrowse genome browser are available at www.fireflybase.org. Raw genomic and RNA-Seq reads for P. pyralis, A. lateralis, and I. luminosus, are available under the NCBI/EBI/DDBJ BioProjects PRJNA378805, PRJDB6460, and PRJNA418169 respectively. Raw WGBS reads can be found on the NCBI Gene Expression Omnibus (GSE107177). Mitochondrial genomes for P. pyralis and I. luminosus and A. antennatus are available on NCBI GenBank with accessions KY778696, MG242621, and MG546669. The complete genome of Entomoplasma luminosum subsp. pyralis is available on NCBI GenBank with accession CP027019. The viral genomes for Photinus pyralis orthomyxo-like virus 1 and 2 are available on NCBI Genbank with accessions MG972985-MG972994. LC-MS data is available on MetaboLights (Accession MTBLS698). Other supporting datasets are available on FigShare (Appendix 6.1).